On the Validation of Protein Force Fields Based on Structural Criteria.

Molecular dynamics simulations depend critically on the quality of the force field used to describe the interatomic interactions and the extent to which it has been validated for use in a specific application. Using a curated test set of 52 high-resolution structures, 39 derived from X-ray diffraction and 13 solved using NMR, we consider the extent to which different parameter sets of the GROMOS protein force field can be distinguished based on comparing a range of structural criteria, including the number of backbone hydrogen bonds, the number of native hydrogen bonds, polar and nonpolar solvent-accessible surface area, radius of gyration, the prevalence of secondary structure elements, J-coupling constants, nuclear Overhauser effect (NOE) intensities, positional root-mean-square deviations (RMSD), and the distribution of backbone ϕ and ψ dihedral angles. It is shown that while statistically significant differences between the average values of individual metrics could be detected, these were in general small. Furthermore, improvements in agreement in one metric were often offset by loss of agreement in another. The work establishes a framework and test set against which protein force fields can be validated. It also highlights the danger of inferring the relative quality of a given force field based on a small range of structural properties or small number of proteins.

Molecular Insights into the Dynamics of Amyloid Fibril Growth: Elongation and Lateral Assembly of GNNQQNY Protofibrils.

The self-assembly of peptides and proteins into ß-sheet rich amyloid fibrils is linked to both functional and pathological states. In this study, the growth of fibrillar structures of the short peptide GNNQQNY, a fragment from the yeast prion Sup35 protein, was examined. Molecular dynamics simulations were used to study alternative mechanisms of fibril growth, including elongation through binding of monomers as well as fibril self-assembly into larger, more mature structures. It was found that after binding, monomers diffused along preformed fibrils toward the ends, supporting the mechanism of fibril growth via elongation. Lateral assembly of protofibrils was found to occur readily, suggesting that this could be the key to transitioning from isolated fibrils to mature multilayer structures. Overall, the work provides mechanistic insights into the competitive pathways that govern amyloid fibril growth.

OFraMP: a fragment-based tool to facilitate the parametrization of large molecules.

An Online tool for Fragment-based Molecule Parametrization (OFraMP) is described. OFraMP is a web application for assigning atomic interaction parameters to large molecules by matching sub-fragments within the target molecule to equivalent sub-fragments within the Automated Topology Builder (ATB, atb.uq.edu.au) database. OFraMP identifies and compares alternative molecular fragments from the ATB database, which contains over 890,000 pre-parameterized molecules, using a novel hierarchical matching procedure. Atoms are considered within the context of an extended local environment (buffer region) with the degree of similarity between an atom in the target molecule and that in the proposed match controlled by varying the size of the buffer region. Adjacent matching atoms are combined into progressively larger matched sub-structures. The user then selects the most appropriate match. OFraMP also allows users to manually alter interaction parameters and automates the submission of missing substructures to the ATB in order to generate parameters for atoms in environments not represented in the existing database. The utility of OFraMP is illustrated using the anti-cancer agent paclitaxel and a dendrimer used in organic semiconductor devices. OFraMP applied to paclitaxel (ATB ID 35922).

Engineering Transferable Atomic Force Fields: Empirical Optimization of Hydrocarbon Lennard-Jones Interactions by Direct Mapping of Parameter Space.

The utility of atomistic simulations depends on the accuracy of the force field used to represent the potential energy landscape, the consistency with which interaction parameters can be assigned, and the extent to which parameters can be transferred between chemical entities. Here, parameter space mapping, a simple and robust procedure for atom typing (parameter assignment) and parameter optimization, is used to identify a minimal set of parameters capable of simultaneously reproducing the density, heat of vaporization, and solvation free energies for a targeted set of simple hydrocarbons. Using an atom-centered fixed charge model and a 6-12 Lennard-Jones potential, the experimental densities and the heats of vaporization for 22 hydrocarbons (linear, cyclic, and aromatic) could be predicted with high precision: average unsigned error (AUE) of 6.1 kg/m3 and 0.5 kJ/mol, respectively, and R2 values of 0.991 and 0.999, respectively. For the 17 compounds with experimental solvation free energy values in water, the AUE was 1.3 kJ/mol, and the slope and R2 for the line of best fit were 0.968 and 0.991, respectively. A key element in ensuring transferability in this work was minimizing confounding variables by ensuring that the calculation of observables was independent of the precise choice of simulation settings (cutoff, bond constraints, etc.) and the explicit consideration of correlations between parameters.

PyThinFilm: Automated Molecular Dynamics Simulation Protocols for the Generation of Thin Film Morphologies.

The performance of organic optoelectronic devices, such as organic light-emitting diodes (OLEDs) and organic solar cells (OSCs), is intrinsically related to the molecular-scale morphology of the thin films from which they are composed. However, the experimental characterization of morphology at the molecular level is challenging due to the often amorphous or at best semicrystalline nature of these films. Classical molecular modeling techniques, such as molecular dynamics (MD) simulation, are increasingly used to understand the relationship between morphology and the properties of thin-film devices. PyThinFilm (github.com/ATB-UQ/PyThinFilm) is an open-source Python package which allows fully automated MD simulations of thin film growth to be performed using vacuum and/or solution deposition processes. PyThinFilm utilizes the GROMACS simulation package in combination with interaction parameters from the Automated Topology Builder (atb.uq.edu.au). Here, PyThinFilm is described along with an overview of applications in which PyThinFilm has been used to study the thin films of organic semiconductor materials typically used in OLEDs and OSCs.

Understanding the performance differences between solution and vacuum deposited OLEDs: A computational approach.

Solution-processing of organic light-emitting diode films has potential advantages in terms of cost and scalability over vacuum-deposition for large area applications. However, solution processed small molecule films can have lower overall device performance. Here, novel molecular dynamics techniques are developed to enable faster simulation of solvent evaporation that occurs during solution processing and give films of thicknesses relevant to real devices. All-atom molecular dynamics simulations are then used in combination with kinetic Monte Carlo transport modeling to examine how differences in morphology stemming from solution or vacuum film deposition affect charge transport and exciton dynamics in films consisting of light-emitting bis(2-phenylpyridine)(acetylacetonate)iridium(III) [Ir(ppy)2(acac)] guest molecules in a 4,4'-bis(N-carbazolyl)biphenyl host. While the structures of the films deposited from vacuum and solution were found to differ, critically, only minor variations in the transport properties were predicted by the simulations even if trapped solvent was present.

Understanding the Effect of pH on the Solubility and Aggregation Extent of Humic Acid in Solution by Combining Simulation and the Experiment.

Molecular dynamics (MD) simulations were performed to investigate the dynamics of humic acid (HA) in an aqueous solution and the influence of pH, temperature, and HA concentration. The HA model employed in MD simulations was chosen and validated using experimental chemical composition data and Fourier transform infrared (FTIR) spectra. The simulations showed that the HA molecule has a strong propensity to adopt a compact conformation in water independent of pH, while the aggregation of HA was found to be pH-dependent. At high pH, the ionized HAs assembled into a thread-like structure, maximizing contact with water. At low pH, the neutral HAs formed a droplet-like aggregate, minimizing contact with the solvent. The simulation results are consistent with experimental data from dynamic light scattering (DLS) measurements and transmission electron microscopy (TEM) imaging. This work provides new insight into the folding and aggregation of HA as a function of pH and a molecular-level understanding of the relationship between the acidity and the structure, solubility, and aggregation of HA, with direct implications for HA-based remediation strategies of contaminated sites.

Unraveling exciton processes in Ir(ppy)3:CBP OLED films upon photoexcitation.

Emissive layers in phosphorescent organic light-emitting diodes commonly make use of guest-host blends such as Ir(ppy)3:CBP to achieve high external quantum efficiencies. However, while the Ir(ppy)3:CBP blend has been studied experimentally, crucial questions remain regarding how exciton diffusion is dependent on the distribution of the guest in the host, which can currently only be addressed at the atomic level via computational modeling. In this work, kinetic Monte Carlo simulations are utilized to gain insight into exciton diffusion in Ir(ppy)3:CBP blend films. The effects of both guest concentration and exciton density on various system properties are analyzed, including the probability of singlet excitons being converted to triplets, and the probability of those triplets decaying radiatively. Significantly, these simulations suggest that triplet diffusion occurs almost exclusively via guest-guest Dexter transfer and that concentration quenching of triplets induced by guest-guest intermolecular dipole-dipole interactions has a negligible effect at high exciton densities due to the prevalence of triplet-triplet annihilation. Furthermore, results for vacuum deposited morphologies derived from molecular dynamics simulations are compared to the results obtained using a simple cubic lattice approximation with randomly distributed guest molecules. We show that while differences in host-based processes such as singlet diffusion are observed, overall, the results on the fate of the excitons are in good agreement for the two morphology types, particularly for guest-based processes at low guest concentrations where guest clustering is limited.

On the Effect of the Various Assumptions and Approximations used in Molecular Simulations on the Properties of Bio-Molecular Systems: Overview and Perspective on Issues.

Computer simulations of molecular systems enable structure-energy-function relationships of molecular processes to be described at the sub-atomic, atomic, supra-atomic or supra-molecular level and plays an increasingly important role in chemistry, biology and physics. To interpret the results of such simulations appropriately, the degree of uncertainty and potential errors affecting the calculated properties must be considered. Uncertainty and errors arise from (1) assumptions underlying the molecular model, force field and simulation algorithms, (2) approximations implicit in the interatomic interaction function (force field), or when integrating the equations of motion, (3) the chosen values of the parameters that determine the accuracy of the approximations used, and (4) the nature of the system and the property of interest. In this overview, advantages and shortcomings of assumptions and approximations commonly used when simulating bio-molecular systems are considered. What the developers of bio-molecular force fields and simulation software can do to facilitate and broaden research involving bio-molecular simulations is also discussed.

Understanding the Activated Form of a Class-I Fusion Protein: Modeling the Interaction of the Ebola Virus Glycoprotein 2 with a Lipid Bilayer.

The fusion of the viral and target cell membranes is a key step in the life cycle of all enveloped viruses. Here, a range of structural data is used to generate an evidence-based model of the active conformation of an archetypical type-I fusion protein, the Ebola glycoprotein 2 (GP2). The stability of the trimeric complex is demonstrated using molecular dynamics and validated by simulating the interaction of the complex with a lipid bilayer. In this model, the fusion peptides project away from the central helix bundle parallel to the target membrane. This maximizes contact with the host membrane, enhances lateral stability, and would explain why, when activated, viral fusion proteins are trimeric.

Evolution and Morphology of Thin Films Formed by Solvent Evaporation: An Organic Semiconductor Case Study.

The crucial role played by the solution-vapor interface in determining the growth and morphology of an organic semiconductor thin film formed by solvent evaporation has been examined in atomic detail. Specifically, how the loss of individual solvent molecules from the surface of the solution induces solute assembly has been studied using molecular dynamics simulations. The system consisted of bis(2-phenylpyridine) (acetylacetonate)iridium(III) [Ir(ppy)2(acac)] and 4,4'-bis(N-carbazolyl)-1,1'-biphenyl (CBP) in chloroform at 310 K. The simulations clearly indicate that (a) the system does not undergo uniform phase separation (spinodal decomposition), (b) solute aggregation initiates at the solution-vapor interface, (c) the distribution of solvent in the film is nonhomogeneous, (d) this nonhomogeneous distribution can induce preferential alignment of host molecules, and (e) a portion of the solvent likely remains trapped within the film. The work not only demonstrates the ability to directly model evaporation in atomic detail on the relevant length scales but also shows that atomistic simulations have the potential to shed new light on morphological properties of a wide range of organic semiconductor devices manufactured using solution-processing methods.

Curved or linear? Predicting the 3-dimensional structure of α-helical antimicrobial peptides in an amphipathic environment.

α-Helical membrane-active antimicrobial peptides (AMPs) are known to act via a range of mechanisms, including the formation of barrel-stave and toroidal pores and the micellisation of the membrane (carpet mechanism). Different mechanisms imply that the peptides adopt different 3D structures when bound at the water-membrane interface, a highly amphipathic environment. Here, an evolutionary algorithm is used to predict the 3D structure of a range of α-helical membrane-active AMPs at the water-membrane interface by optimising amphipathicity. This amphipathic structure prediction (ASP) is capable of distinguishing between curved and linear peptides solved experimentally, potentially allowing the activity and mechanism of action of different membrane-active AMPs to be predicted. The ASP algorithm is accessible via a web interface at http://atb.uq.edu.au/asp/.

Effect of Triclosan and Chloroxylenol on Bacterial Membranes.

Triclosan and chloroxylenol are broad-spectrum biocides used extensively in healthcare and consumer products. They have been suggested to perturb the structure of bacterial membranes, but studies so far have not considered that most bacterial membranes contain large amounts of branched-chain lipids. Here, molecular dynamics simulation is used to examine the effect of the two biocides on membranes consisting of lipids with methyl-branched chains, cyclopropanated chains, and nonbranched chains. It is shown that triclosan and chloroxylenol induced a phase transition in membranes from a liquid-crystalline to a liquid-ordered phase irrespective of the presence and nature of branching groups. At high concentration, chloroxylenol promoted chain interdigitation. Our results suggest that triclosan and chloroxylenol decrease the degree of fluidity of membranes and that this effect is more pronounced in bacterial membranes. As a result, their biocidal activity could be associated with a change in the function of membrane proteins.

Response of microbial membranes to butanol: interdigitation vs. disorder.

Biobutanol production by fermentation is potentially a sustainable alternative to butanol production from fossil fuels. However, the toxicity of butanol to fermentative bacteria, resulting largely from cell membrane fluidization, limits production titers and is a major factor limiting the uptake of the technology. Here, studies were undertaken, in vitro and in silico, on the butanol effects on a representative bacterial (i.e. Escherichia coli) inner cell membrane. A critical butanol : lipid ratio for stability of 2 : 1 was observed, computationally, consistent with complete interdigitation. However, at this ratio the bilayer was ∼20% thicker than for full interdigitation. Furthermore, butanol intercalation induced acyl chain bending and increased disorder, measured as a 27% lateral diffusivity increase experimentally in a supported lipid bilayer. There was also a monophasic Tm reduction in butanol-treated large unilamellar vesicles. Both behaviours are inconsistent with an interdigitated gel. Butanol thus causes only partial interdigitation at physiological temperatures, due to butanol accumulating at the phospholipid headgroups. Acyl tail disordering (i.e. splaying and bending) fills the subsequent voids. Finally, butanol short-circuits the bilayer and creates a coupled system where interdigitated and splayed phospholipids coexist. These findings will inform the design of strategies targeting bilayer stability for increasing biobutanol production titers.

Automated partial atomic charge assignment for drug-like molecules: a fast knapsack approach.

A key factor in computational drug design is the consistency and reliability with which intermolecular interactions between a wide variety of molecules can be described. Here we present a procedure to efficiently, reliably and automatically assign partial atomic charges to atoms based on known distributions. We formally introduce the molecular charge assignment problem, where the task is to select a charge from a set of candidate charges for every atom of a given query molecule. Charges are accompanied by a score that depends on their observed frequency in similar neighbourhoods (chemical environments) in a database of previously parameterised molecules. The aim is to assign the charges such that the total charge equals a known target charge within a margin of error while maximizing the sum of the charge scores. We show that the problem is a variant of the well-studied multiple-choice knapsack problem and thus weakly NP -complete. We propose solutions based on Integer Linear Programming and a pseudo-polynomial time Dynamic Programming algorithm. We demonstrate that the results obtained for novel molecules not included in the database are comparable to the ones obtained performing explicit charge calculations while decreasing the time to determine partial charges for a molecule from hours or even days to below a second. Our software is openly available.

Probing the Pharmacological Binding Sites of P-Glycoprotein Using Umbrella Sampling Simulations.

The human multidrug transporter P-glycoprotein (P-gp) transports over 200 chemically diverse substrates, influencing their bioavailability and tissue distribution. Pharmacological studies have identified both competitive and noncompetitive P-gp substrates, but neither the precise location of the substrate binding sites, nor the basis of competitive and noncompetitive interactions has been fully characterized. Here, potential of mean force (PMF) calculations are used to identify the transport-competent minimum free energy binding locations of five compounds, Hoechst 33342, Rhodamine 123, paclitaxel, tariquidar, and verapamil to P-gp. Unrestrained molecular dynamics simulations were also performed to confirm the substrates were stable in the energy wells determined using the PMF calculations. All compounds had energy minima within the P-gp transmembrane (TM) pore. For Hoechst 33342 and Rhodamine 123, a second minimum outside the TM pore was also identified. Based on this and previous studies of nicardipine and morphine [ Subramanian et al. J. Chem. Inf. Model. 2015 , 55 , 1202 ], a general scheme that accounts for the observed noncompetitive and competitive substrate interactions with P-gp is proposed.

Automated Topology Builder Version 3.0: Prediction of Solvation Free Enthalpies in Water and Hexane.

The ability of atomic interaction parameters generated using the Automated Topology Builder and Repository version 3.0 (ATB3.0) to predict experimental hydration free enthalpies (Δ Gwater) and solvation free enthalpies in the apolar solvent hexane (Δ Ghexane) is presented. For a validation set of 685 molecules the average unsigned error (AUE) between Δ Gwater values calculated using the ATB3.0 and experiment is 3.8 kJ·mol-1. The slope of the line of best fit is 1.00, the intercept -1.0 kJ·mol-1, and the R2 0.90. For the more restricted set of 239 molecules used to validate OPLS3 ( J. Chem. Theory Comput. 2016 , 12 , 281 - 296 , DOI: 10.1021/acs.jctc.5b00864 ) the AUE using the ATB3.0 is just 2.7 kJ·mol-1 and the R2 0.93. A roadmap for further improvement of the ATB parameters is presented together with a discussion of the challenges of validating force fields against the available experimental data.

Morphology of a Bulk Heterojunction Photovoltaic Cell with Low Donor Concentration.

Atomistic nonequilibrium molecular dynamics simulations have been used to model the morphology of small-molecule bulk heterojunction films formed by vapor deposition as used in organic photovoltaics. Films comprising C60 and 1, 5, 10, and 50 wt % of 1,1-bis[4-bis(4-methylphenyl)aminophenyl]cyclohexane (TAPC) were compared to films of neat C60. The simulations suggest that if holes can hop between donor molecules separated by as little as 1.2-1.5 nm, then a TAPC concentration of 5 wt % is sufficient to form a percolating donor network and facilitate charge extraction. The results provide an explanation for why low donor content organic photovoltaics can still have high efficiencies. In addition, the roughness, porosity, and crystallinity of the films were found to decrease with increasing TAPC content.

Predicting the Prevalence of Alternative Warfarin Tautomers in Solution.

Warfarin, a widely used oral anticoagulant, is prescribed as a racemic mixture. Each enantiomer of neutral Warfarin can exist in 20 possible tautomeric states leading to complex pharmacokinetics and uncertainty as to the relevant species under different conditions. Here, the ability of alternative computational approaches to predict the preferred tautomeric form(s) of neutral Warfarin in different solvents is examined. It is shown that varying the method used to estimate the heat of formation in vacuum (direct or via homodesmic reactions), whether entropic corrections were included, and the method used to estimate the free enthalpy of solvation (i.e., PCM, COSMO, or SMD implicit models or explicit solvent) lead to large differences in the predicted rank and relative populations of the tautomers. In this case, only a combination of the enthalpy of formation using homodesmic reactions and explicit solvent to estimate the free enthalpy of solvation yielded results compatible with the available experimental data. The work also suggests that a small but significant subset of the possible Warfarin tautomers are likely to be physiologically relevant.

Effect of Binding on Enantioselectivity of Epoxide Hydrolase.

Molecular dynamics simulations and free energy calculations have been used to investigate the effect of ligand binding on the enantioselectivity of an epoxide hydrolase (EH) from Aspergillus niger. Despite sharing a common mechanism, a wide range of alternative mechanisms have been proposed to explain the origin of enantiomeric selectivity in EHs. By comparing the interactions of ( R)- and ( S)-glycidyl phenyl ether (GPE) with both the wild type (WT, E = 3) and a mutant showing enhanced enantioselectivity to GPE (LW202, E = 193), we have examined whether enantioselectivity is due to differences in the binding pose, the affinity for the ( R)- or ( S)- enantiomers, or a kinetic effect. The two enantiomers were easily accommodated within the binding pockets of the WT enzyme and LW202. Free energy calculations suggested that neither enzyme had a preference for a given enantiomer. The two substrates sampled a wide variety of conformations in the simulations with the sterically hindered and unhindered carbon atoms of the GPE epoxide ring both coming in close proximity to the nucleophilic aspartic acid residue. This suggests that alternative pathways could lead to the formation of a ( S)- and ( R)-diol product. Together, the calculations suggest that the enantioselectivity is due to kinetic rather than thermodynamic effects and that the assumption that one substrate results in one product when interpreting the available experimental data and deriving E-values may be inappropriate in the case of EHs.