Data-Driven Mapping of Gas-Phase Quantum Calculations to General Force Field Lennard-Jones Parameters.

Molecular dynamics simulations are helpful tools for a range of applications, ranging from drug discovery to protein structure determination. The successful use of this technology largely depends on the potential function, or force field, used to determine the potential energy at each configuration of the system. Most force fields encode all of the relevant parameters to be used in distinct atom types, each associated with parameters for all parts of the force field, typically bond stretches, angle bends, torsions, and nonbonded terms accounting for van der Waals and electrostatic interactions. Much attention has been paid to the nonbonded parameters and their derivation, which are important in particular due to their governance of noncovalent interactions, such as protein-ligand binding. Parametrization involves adjusting the nonbonded parameters to minimize the error between simulation results and experimental properties, such as heats of vaporization and densities of neat liquids. In this setting, determining the best set of atom types is far from trivial, and the large number of parameters to be fit for the atom types in a typical force field can make it difficult to approach a true optimum. Here, we utilize a previously described Minimal Basis Iterative Stockholder (MBIS) method to carry out an atoms-in-molecules partitioning of electron densities. Information from these atomic densities is then mapped to Lennard-Jones parameters using a set of mapping parameters much smaller than the typical number of atom types in a force field. This approach is advantageous in two ways: it eliminates atom types by allowing each atom to have unique Lennard-Jones parameters, and it greatly reduces the number of parameters to be optimized. We show that this approach yields results comparable to those obtained with the typed GAFF 1.7 force field, even when trained on a relatively small amount of experimental data.

Bind3P: Optimization of a Water Model Based on Host-Guest Binding Data.

We report a water model, Bind3P (Version 0.1), which was obtained by using sensitivity analysis to readjust the Lennard-Jones parameters of the TIP3P model against experimental binding free energies for six host-guest systems, along with pure liquid properties. Tests of Bind3P against >100 experimental binding free energies and enthalpies for host-guest systems distinct from the training set show a consistent drop in the mean signed error, relative to matched calculations with TIP3P. Importantly, Bind3P also yields some improvement in the hydration free energies of small organic molecules and preserves the accuracy of bulk water properties, such as density and the heat of vaporization. The same approach can be applied to more sophisticated water models that can better represent pure water properties. These results lend further support to the concept of integrating host-guest binding data into force field parametrization.

Selectively manipulable acoustic-powered microswimmers.

Selective actuation of a single microswimmer from within a diverse group would be a first step toward collaborative guided action by a group of swimmers. Here we describe a new class of microswimmer that accomplishes this goal. Our swimmer design overcomes the commonly-held design paradigm that microswimmers must use non-reciprocal motion to achieve propulsion; instead, the swimmer is propelled by oscillatory motion of an air bubble trapped within the swimmer's polymer body. This oscillatory motion is driven by the application of a low-power acoustic field, which is biocompatible with biological samples and with the ambient liquid. This acoustically-powered microswimmer accomplishes controllable and rapid translational and rotational motion, even in highly viscous liquids (with viscosity 6,000 times higher than that of water). And by using a group of swimmers each with a unique bubble size (and resulting unique resonance frequencies), selective actuation of a single swimmer from among the group can be readily achieved.

Calculation and visualization of atomistic mechanical stresses in nanomaterials and biomolecules.

Many biomolecules have machine-like functions, and accordingly are discussed in terms of mechanical properties like force and motion. However, the concept of stress, a mechanical property that is of fundamental importance in the study of macroscopic mechanics, is not commonly applied in the biomolecular context. We anticipate that microscopical stress analyses of biomolecules and nanomaterials will provide useful mechanistic insights and help guide molecular design. To enable such applications, we have developed Calculator of Atomistic Mechanical Stress (CAMS), an open-source software package for computing atomic resolution stresses from molecular dynamics (MD) simulations. The software also enables decomposition of stress into contributions from bonded, nonbonded and Generalized Born potential terms. CAMS reads GROMACS topology and trajectory files, which are easily generated from AMBER files as well; and time-varying stresses may be animated and visualized in the VMD viewer. Here, we review relevant theory and present illustrative applications.

Acoustofluidic chemical waveform generator and switch.

Eliciting a cellular response to a changing chemical microenvironment is central to many biological processes including gene expression, cell migration, differentiation, apoptosis, and intercellular signaling. The nature and scope of the response is highly dependent upon the spatiotemporal characteristics of the stimulus. To date, studies that investigate this phenomenon have been limited to digital (or step) chemical stimulation with little control over the temporal counterparts. Here, we demonstrate an acoustofluidic (i.e., fusion of acoustics and microfluidics) approach for generating programmable chemical waveforms that permits continuous modulation of the signal characteristics including the amplitude (i.e., sample concentration), shape, frequency, and duty cycle, with frequencies reaching up to 30 Hz. Furthermore, we show fast switching between multiple distinct stimuli, wherein the waveform of each stimulus is independently controlled. Using our device, we characterized the frequency-dependent activation and internalization of the ß2-adrenergic receptor (ß2-AR), a prototypic G-protein coupled receptor (GPCR), using epinephrine. The acoustofluidic-based programmable chemical waveform generation and switching method presented herein is expected to be a powerful tool for the investigation and characterization of the kinetics and other dynamic properties of many biological and biochemical processes.

Bridging Calorimetry and Simulation through Precise Calculations of Cucurbituril-Guest Binding Enthalpies.

We used microsecond time scale molecular dynamics simulations to compute, at high precision, binding enthalpies for cucurbit[7]uril (CB7) with eight guests in aqueous solution. The results correlate well with experimental data from previously published isothermal titration calorimetry studies, and decomposition of the computed binding enthalpies by interaction type provides plausible mechanistic insights. Thus, dispersion interactions appear to play a key role in stabilizing these complexes, due at least in part to the fact that their packing density is greater than that of water. On the other hand, strongly favorable Coulombic interactions between the host and guests are compensated by unfavorable solvent contributions, leaving relatively modest electrostatic contributions to the binding enthalpies. The better steric fit of the aliphatic guests into the circular host appears to explain why their binding enthalpies tend to be more favorable than those of the more planar aromatic guests. The present calculations also bear on the validity of the simulation force field. Somewhat unexpectedly, the TIP3P water yields better agreement with experiment than the TIP4P-Ew water model, although the latter is known to replicate the properties of pure water more accurately. More broadly, the present results demonstrate the potential for computational calorimetry to provide atomistic explanations for thermodynamic observations.

Quantum Mechanical Calculation of Noncovalent Interactions: A Large-Scale Evaluation of PMx, DFT, and SAPT Approaches.

Quantum mechanical (QM) calculations of noncovalent interactions are uniquely useful as tools to test and improve molecular mechanics force fields and to model the forces involved in biomolecular binding and folding. Because the more computationally tractable QM methods necessarily include approximations, which risk degrading accuracy, it is essential to evaluate such methods by comparison with high-level reference calculations. Here, we use the extensive Benchmark Energy and Geometry Database (BEGDB) of CCSD(T)/CBS reference results to evaluate the accuracy and speed of widely used QM methods for over 1200 chemically varied gas-phase dimers. In particular, we study the semiempirical PM6 and PM7 methods; density functional theory (DFT) approaches B3LYP, B97-D, M062X, and ωB97X-D; and symmetry-adapted perturbation theory (SAPT) approach. For the PM6 and DFT methods, we also examine the effects of post hoc corrections for hydrogen bonding (PM6-DH+, PM6-DH2), halogen atoms (PM6-DH2X), and dispersion (DFT-D3 with zero and Becke-Johnson damping). Several orders of the SAPT expansion are also compared, ranging from SAPT0 up to SAPT2+3, where computationally feasible. We find that all DFT methods with dispersion corrections, as well as SAPT at orders above SAPT2, consistently provide dimer interaction energies within 1.0 kcal/mol RMSE across all systems. We also show that a linear scaling of the perturbative energy terms provided by the fast SAPT0 method yields similar high accuracy, at particularly low computational cost. The energies of all the dimer systems from the various QM approaches are included in the Supporting Information, as are the full SAPT2+(3) energy decomposition for a subset of over 1000 systems. The latter can be used to guide the parametrization of molecular mechanics force fields on a term-by-term basis.

The SAMPL4 host-guest blind prediction challenge: an overview.

Prospective validation of methods for computing binding affinities can help assess their predictive power and thus set reasonable expectations for their performance in drug design applications. Supramolecular host-guest systems are excellent model systems for testing such affinity prediction methods, because their small size and limited conformational flexibility, relative to proteins, allows higher throughput and better numerical convergence. The SAMPL4 prediction challenge therefore included a series of host-guest systems, based on two hosts, cucurbit[7]uril and octa-acid. Binding affinities in aqueous solution were measured experimentally for a total of 23 guest molecules. Participants submitted 35 sets of computational predictions for these host-guest systems, based on methods ranging from simple docking, to extensive free energy simulations, to quantum mechanical calculations. Over half of the predictions provided better correlations with experiment than two simple null models, but most methods underperformed the null models in terms of root mean squared error and linear regression slope. Interestingly, the overall performance across all SAMPL4 submissions was similar to that for the prior SAMPL3 host-guest challenge, although the experimentalists took steps to simplify the current challenge. While some methods performed fairly consistently across both hosts, no single approach emerged as consistent top performer, and the nonsystematic nature of the various submissions made it impossible to draw definitive conclusions regarding the best choices of energy models or sampling algorithms. Salt effects emerged as an issue in the calculation of absolute binding affinities of cucurbit[7]uril-guest systems, but were not expected to affect the relative affinities significantly. Useful directions for future rounds of the challenge might involve encouraging participants to carry out some calculations that replicate each others' studies, and to systematically explore parameter options.

Blind prediction of SAMPL4 cucurbit[7]uril binding affinities with the mining minima method.

Accurate methods for predicting protein-ligand binding affinities are of central interest to computer-aided drug design for hit identification and lead optimization. Here, we used the mining minima (M2) method to predict cucurbit[7]uril binding affinities from the SAMPL4 blind prediction challenge. We tested two different energy models, an empirical classical force field, CHARMm with VCharge charges, and the Poisson-Boltzmann surface area solvation model; and a semiempirical quantum mechanical (QM) Hamiltonian, PM6-DH+, coupled with the COSMO solvation model and a surface area term for nonpolar solvation free energy. Binding affinities based on the classical force field correlated strongly with the experiments with a correlation coefficient (R(2)) of 0.74. On the other hand, binding affinities based on the QM energy model correlated poorly with experiments (R(2) = 0.24), due largely to two major outliers. As we used extensive conformational search methods, these results point to possible inaccuracies in the PM6-DH+ energy model or the COSMO solvation model. Furthermore, the different binding free energy components, solute energy, solvation free energy, and configurational entropy showed significant deviations between the classical M2 and quantum M2 calculations. Comparison of different classical M2 free energy components to experiments show that the change in the total energy, i.e. the solute energy plus the solvation free energy, is the key driving force for binding, with a reasonable correlation to experiment (R(2) = 0.56); however, accounting for configurational entropy further improves the correlation.

The SAMPL4 hydration challenge: evaluation of partial charge sets with explicit-water molecular dynamics simulations.

We used blind predictions of the 47 hydration free energies in the SAMPL4 challenge to test multiple partial charge models in the context of explicit solvent free energy simulations with the general AMBER force field. One of the partial charge models, IPolQ-Mod, is a fast continuum solvent-based implementation of the IPolQ approach. The AM1-BCC, restrained electrostatic potential (RESP) and IpolQ-Mod approaches all perform reasonably well (R(2) > 0.8), while VCharge, though faster, gives less accurate results (R(2) of 0.5). The AM1-BCC results are more accurate than those of RESP for tertiary amines and nitrates, but the overall difference in accuracy between these methods is not statistically significant. Interestingly, the IPolQ-Mod method is found to yield partial charges in very close agreement with RESP. This observation suggests that the success of RESP may be attributed to its fortuitously approximating the arguably more rigorous IPolQ approach.

The electrostatic response of water to neutral polar solutes: implications for continuum solvent modeling.

Continuum solvation models are widely used to estimate the hydration free energies of small molecules and proteins, in applications ranging from drug design to protein engineering, and most such models are based on the approximation of a linear dielectric response by the solvent. We used explicit-water molecular dynamics simulations with the TIP3P water model to probe this linear response approximation in the case of neutral polar molecules, using miniature cucurbituril and cyclodextrin receptors and protein side-chain analogs as model systems. We observe supralinear electrostatic solvent responses, and this nonlinearity is found to result primarily from waters' being drawn closer and closer to the solutes with increased solute-solvent electrostatic interactions; i.e., from solute electrostriction. Dielectric saturation and changes in the water-water hydrogen bonding network, on the other hand, play little role. Thus, accounting for solute electrostriction may be a productive approach to improving the accuracy of continuum solvation models.

Substrate stiffness regulates cellular uptake of nanoparticles.

Nanoparticle (NP)-bioconjugates hold great promise for more sensitive disease diagnosis and more effective anticancer drug delivery compared with existing approaches. A critical aspect in both applications is cellular internalization of NPs, which is influenced by NP properties and cell surface mechanics. Despite considerable progress in optimization of the NP-bioconjugates for improved targeting, the role of substrate stiffness on cellular uptake has not been investigated. Using polyacrylamide (PA) hydrogels as model substrates with tunable stiffness, we quantified the relationship between substrate stiffness and cellular uptake of fluorescent NPs by bovine aortic endothelial cells (BAECs). We found that a stiffer substrate results in a higher total cellular uptake on a per cell basis, but a lower uptake per unit membrane area. To obtain a mechanistic understanding of the cellular uptake behavior, we developed a thermodynamic model that predicts that membrane spreading area and cell membrane tension are two key factors controlling cellular uptake of NPs, both of which are modulated by substrate stiffness. Our experimental and modeling results not only open up new avenues for engineering NP-based cancer cell targets for more effective in vivo delivery but also contribute an example of how the physical environment dictates cellular behavior and function.

Enzyme molecules as nanomotors.

Using fluorescence correlation spectroscopy, we show that the diffusive movements of catalase enzyme molecules increase in the presence of the substrate, hydrogen peroxide, in a concentration-dependent manner. Employing a microfluidic device to generate a substrate concentration gradient, we show that both catalase and urease enzyme molecules spread toward areas of higher substrate concentration, a form of chemotaxis at the molecular scale. Using glucose oxidase and glucose to generate a hydrogen peroxide gradient, we induce the migration of catalase toward glucose oxidase, thereby showing that chemically interconnected enzymes can be drawn together.

Tunable, pulsatile chemical gradient generation via acoustically driven oscillating bubbles.

We present a novel concept of generating both static and pulsatile chemical gradients using acoustically activated bubbles designed in a ladder-like arrangement. Furthermore, by regulating the amplitude of the bubble oscillation, we demonstrate that the chemical gradient profiles can be effectively tuned.

Entropy-enthalpy transduction caused by conformational shifts can obscure the forces driving protein-ligand binding.

Molecular dynamics simulations of unprecedented duration now can provide new insights into biomolecular mechanisms. Analysis of a 1-ms molecular dynamics simulation of the small protein bovine pancreatic trypsin inhibitor reveals that its main conformations have different thermodynamic profiles and that perturbation of a single geometric variable, such as a torsion angle or interresidue distance, can select for occupancy of one or another conformational state. These results establish the basis for a mechanism that we term entropy-enthalpy transduction (EET), in which the thermodynamic character of a local perturbation, such as enthalpic binding of a small molecule, is camouflaged by the thermodynamics of a global conformational change induced by the perturbation, such as a switch into a high-entropy conformational state. It is noted that EET could occur in many systems, making measured entropies and enthalpies of folding and binding unreliable indicators of actual thermodynamic driving forces. The same mechanism might also account for the high experimental variance of measured enthalpies and entropies relative to free energies in some calorimetric studies. Finally, EET may be the physical mechanism underlying many cases of entropy-enthalpy compensation.

Calculation of Host-Guest Binding Affinities Using a Quantum-Mechanical Energy Model.

The prediction of protein-ligand binding affinities is of central interest in computer-aided drug discovery, but it is still difficult to achieve a high degree of accuracy. Recent studies suggesting that available force fields may be a key source of error motivate the present study, which reports the first mining minima (M2) binding affinity calculations based on a quantum mechanical energy model, rather than an empirical force field. We apply a semi-empirical quantum-mechanical energy function, PM6-DH+, coupled with the COSMO solvation model, to 29 host-guest systems with a wide range of measured binding affinities. After correction for a systematic error, which appears to derive from the treatment of polar solvation, the computed absolute binding affinities agree well with experimental measurements, with a mean error 1.6 kcal/mol and a correlation coefficient of 0.91. These calculations also delineate the contributions of various energy components, including solute energy, configurational entropy, and solvation free energy, to the binding free energies of these host-guest complexes. Comparison with our previous calculations, which used empirical force fields, point to significant differences in both the energetic and entropic components of the binding free energy. The present study demonstrates successful combination of a quantum mechanical Hamiltonian with the M2 affinity method.

Blind prediction of host-guest binding affinities: a new SAMPL3 challenge.

The computational prediction of protein-ligand binding affinities is of central interest in early-stage drug-discovery, and there is a widely recognized need for improved methods. Low molecular weight receptors and their ligands--i.e., host-guest systems--represent valuable test-beds for such affinity prediction methods, because their small size makes for fast calculations and relatively facile numerical convergence. The SAMPL3 community exercise included the first ever blind prediction challenge for host-guest binding affinities, through the incorporation of 11 new host-guest complexes. Ten participating research groups addressed this challenge with a variety of approaches. Statistical assessment indicates that, although most methods performed well at predicting some general trends in binding affinity, overall accuracy was not high, as all the methods suffered from either poor correlation or high RMS errors or both. There was no clear advantage in using explicit versus implicit solvent models, any particular force field, or any particular approach to conformational sampling. In a few cases, predictions using very similar energy models but different sampling and/or free-energy methods resulted in significantly different results. The protonation states of one host and some guest molecules emerged as key uncertainties beyond the choice of computational approach. The present results have implications for methods development and future blind prediction exercises.

Tuning membrane phase separation using nonlipid amphiphiles.

Lipid phase separation may be a mechanism by which lipids participate in sorting membrane proteins and facilitate membrane-mediated biochemical signaling in cells. To provide new tools for membrane lipid phase manipulation that avoid direct effects on protein activity and lipid composition, we studied phase separation in binary and ternary lipid mixtures under the influence of three nonlipid amphiphiles, vitamin E (VE), Triton-X (TX)-100, and benzyl alcohol (BA). Mechanisms of additive-induced phase separation were elucidated using coarse-grained molecular dynamics simulations of these additives in a liquid bilayer made from 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) and 1,2-dilinoleoyl-sn-glycero-3-phosphocholine [corrected]. From simulations, the additive's partitioning preference, changes in membrane thickness, and alterations in lipid order were quantified. Simulations showed that VE favored the DPPC phase but partitioned predominantly to the domain boundaries and lowered the tendency for domain formation, and therefore acted as a linactant. This simulated behavior was consistent with experimental observations in which VE promoted lipid mixing and dispersed domains in both gel/liquid and liquid-ordered/liquid-disordered systems. From simulation, BA partitioned predominantly to the DUPC phase, decreased lipid order there, and thinned the membrane. These actions explain why, experimentally, BA promoted phase separation in both binary and ternary lipid mixtures. In contrast, TX, a popular detergent used to isolate raft membranes in cells, exhibited equal preference for both phases, as demonstrated by simulations, but nonetheless, was a strong domain promoter in all lipid mixtures. Further analysis showed that TX increased membrane thickness of the DPPC phase to a greater extent than the DUPC phase and thus increased hydrophobic mismatch, which may explain experimental observation of phase separation in the presence of TX. In summary, these nonlipid amphiphiles provide new tools to tune domain formation in model vesicle systems and could provide the means to form or disperse membrane lipid domains in cells, in addition to the well-known methods involving cholesterol enrichment and sequestration.

Prediction of SAMPL3 host-guest binding affinities: evaluating the accuracy of generalized force-fields.

We used the second-generation mining minima method (M2) to compute the binding affinities of the novel host-guest complexes in the SAMPL3 blind prediction challenge. The predictions were in poor agreement with experiment, and we conjectured that much of the error might derive from the force field, CHARMm with Vcharge charges. Repeating the calculations with other generalized force-fields led to no significant improvement, and we observed that the predicted affinities were highly sensitive to the choice of force-field. We therefore embarked on a systematic evaluation of a set of generalized force fields, based upon comparisons with PM6-DH2, a fast yet accurate semi-empirical quantum mechanics method. In particular, we compared gas-phase interaction energies and entropies for the host-guest complexes themselves, as well as for smaller chemical fragments derived from the same molecules. The mean deviations of the force field interaction energies from the quantum results were greater than 3 kcal/mol and 9 kcal/mol, for the fragments and host-guest systems respectively. We further evaluated the accuracy of force-fields for computing the vibrational entropies and found the mean errors to be greater than 4 kcal/mol. Given these errors in energy and entropy, it is not surprising in retrospect that the predicted binding affinities deviated from the experiment by several kcal/mol. These results emphasize the need for improvements in generalized force-fields and also highlight the importance of systematic evaluation of force-field parameters prior to evaluating different free-energy methods.