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.

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).

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.

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.

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.

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.

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.

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.

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.

The CC domain structure from the wheat stem rust resistance protein Sr33 challenges paradigms for dimerization in plant NLR proteins.

Plants use intracellular immunity receptors, known as nucleotide-binding oligomerization domain-like receptors (NLRs), to recognize specific pathogen effector proteins and induce immune responses. These proteins provide resistance to many of the world's most destructive plant pathogens, yet we have a limited understanding of the molecular mechanisms that lead to defense signaling. We examined the wheat NLR protein, Sr33, which is responsible for strain-specific resistance to the wheat stem rust pathogen, Puccinia graminis f. sp. tritici We present the solution structure of a coiled-coil (CC) fragment from Sr33, which adopts a four-helix bundle conformation. Unexpectedly, this structure differs from the published dimeric crystal structure of the equivalent region from the orthologous barley powdery mildew resistance protein, MLA10, but is similar to the structure of the distantly related potato NLR protein, Rx. We demonstrate that these regions are, in fact, largely monomeric and adopt similar folds in solution in all three proteins, suggesting that the CC domains from plant NLRs adopt a conserved fold. However, larger C-terminal fragments of Sr33 and MLA10 can self-associate both in vitro and in planta, and this self-association correlates with their cell death signaling activity. The minimal region of the CC domain required for both cell death signaling and self-association extends to amino acid 142, thus including 22 residues absent from previous biochemical and structural protein studies. These data suggest that self-association of the minimal CC domain is necessary for signaling but is likely to involve a different structural basis than previously suggested by the MLA10 crystallographic dimer.

A potential new, stable state of the E-cadherin strand-swapped dimer in solution.

E-cadherin is a transmembrane glycoprotein that facilitates inter-cellular adhesion in the epithelium. The ectodomain of the native structure is comprised of five repeated immunoglobulin-like domains. All E-cadherin crystal structures show the protein in one of three alternative conformations: a monomer, a strand-swapped trans homodimer and the so-called X-dimer, which is proposed to be a kinetic intermediate to forming the strand-swapped trans homodimer. However, previous studies have indicated that even once the trans strand-swapped dimer is formed, the complex is highly dynamic and the E-cadherin monomers may reorient relative to each other. Here, molecular dynamics simulations have been used to investigate the stability and conformational flexibility of the human E-cadherin trans strand-swapped dimer. In four independent, 100 ns simulations, the dimer moved away from the starting structure and converged to a previously unreported structure, which we call the Y-dimer. The Y-dimer was present for over 90% of the combined simulation time, suggesting that it represents a stable conformation of the E-cadherin dimer in solution. The Y-dimer conformation is stabilised by interactions present in both the trans strand-swapped dimer and X-dimer crystal structures, as well as additional interactions not found in any E-cadherin dimer crystal structures. The Y-dimer represents a previously unreported, stable conformation of the human E-cadherin trans strand-swapped dimer and suggests that the available crystal structures do not fully capture the conformations that the human E-cadherin trans homodimer adopts in solution.

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.

The Molecular Origin of Anisotropic Emission in an Organic Light-Emitting Diode.

Atomistic nonequilibrium molecular dynamics simulations have been used to model the induction of molecular orientation anisotropy within the emission layer of an organic light-emitting diode (OLED) formed by vapor deposition. Two emitter species were compared: racemic fac-tris(2-phenylpyridine)iridium(III) (Ir(ppy)3) and trans-bis(2-phenylpyridine)(acetylacetonate)iridium(III) (Ir(ppy)2(acac)). The simulations show that the molecular symmetry axes of both emitters preferentially align perpendicular to the surface during deposition. The molecular arrangement formed on deposition combined with consideration of the transition dipole moments provides insight into experimental reports that Ir(ppy)3 shows isotropic emission, while Ir(ppy)2(acac) displays improved efficiency due to an apparent preferential alignment of the transition dipole vectors parallel to the substrate. The simulations indicate that this difference is not due to differences in the extent of emitter alignment, but rather differences in the direction of the transition dipoles within the two complexes.

Validation of Molecular Simulation: An Overview of Issues.

Computer simulation of molecular systems enables structure-energy-function relationships of molecular processes to be described at the sub-atomic, atomic, supra-atomic, or supra-molecular level. To interpret results of such simulations appropriately, the quality of the calculated properties must be evaluated. This depends on the way the simulations are performed and on the way they are validated by comparison to values Qexp of experimentally observable quantities Q. One must consider 1) the accuracy of Qexp , 2) the accuracy of the function Q(rN ) used to calculate a Q-value based on a molecular configuration rN of N particles, 3) the sensitivity of the function Q(rN ) to the configuration rN , 4) the relative time scales of the simulation and experiment, 5) the degree to which the calculated and experimental properties are equivalent, and 6) the degree to which the system simulated matches the experimental conditions. Experimental data is limited in scope and generally corresponds to averages over both time and space. A critical analysis of the various factors influencing the apparent degree of (dis)agreement between simulations and experiment is presented and illustrated using examples from the literature. What can be done to enhance the validation of molecular simulation is also discussed.

Do All X-ray Structures of Protein-Ligand Complexes Represent Functional States? EPOR, a Case Study.

Based on differences between the x-ray crystal structures of ligand-bound and unbound forms, the activation of the erythropoietin receptor (EPOR) was initially proposed to involve a cross-action scissorlike motion. However, the validity of the motions involved in the scissorlike model has been recently challenged. Here, atomistic molecular dynamics simulations are used to examine the structure of the extracellular domain of the EPOR dimer in the presence and absence of erythropoietin and a series of agonistic or antagonistic mimetic peptides free in solution. The simulations suggest that in the absence of crystal packing effects, the EPOR chains in the different dimers adopt very similar conformations with no clear distinction between the agonist and antagonist-bound complexes. This questions whether the available x-ray crystal structures of EPOR truly represent active or inactive conformations. The study demonstrates the difficulty in using such structures to infer a mechanism of action, especially in the case of membrane receptors where just part of the structure has been considered in addition to potential confounding effects that arise from the comparison of structures in a crystal as opposed to a membrane environment. The work highlights the danger of assigning functional significance to small differences between structures of proteins bound to different ligands in a crystal environment without consideration of the effects of the crystal lattice and thermal motion.

Interaction of Tarantula Venom Peptide ProTx-II with Lipid Membranes Is a Prerequisite for Its Inhibition of Human Voltage-gated Sodium Channel NaV1.7.

ProTx-II is a disulfide-rich peptide toxin from tarantula venom able to inhibit the human voltage-gated sodium channel 1.7 (hNaV1.7), a channel reported to be involved in nociception, and thus it might have potential as a pain therapeutic. ProTx-II acts by binding to the membrane-embedded voltage sensor domain of hNaV1.7, but the precise peptide channel-binding site and the importance of membrane binding on the inhibitory activity of ProTx-II remain unknown. In this study, we examined the structure and membrane-binding properties of ProTx-II and several analogues using NMR spectroscopy, surface plasmon resonance, fluorescence spectroscopy, and molecular dynamics simulations. Our results show a direct correlation between ProTx-II membrane binding affinity and its potency as an hNaV1.7 channel inhibitor. The data support a model whereby a hydrophobic patch on the ProTx-II surface anchors the molecule at the cell surface in a position that optimizes interaction of the peptide with the binding site on the voltage sensor domain. This is the first study to demonstrate that binding of ProTx-II to the lipid membrane is directly linked to its potency as an hNaV1.7 channel inhibitor.

Validating lipid force fields against experimental data: Progress, challenges and perspectives.

Biological membranes display a great diversity in lipid composition and lateral structure that is crucial in a variety of cellular functions. Simulations of membranes have contributed significantly to the understanding of the properties, functions and behaviour of membranes and membrane-protein assemblies. This success relies on the ability of the force field used to describe lipid-lipid and lipid-environment interactions accurately, reproducibly and realistically. In this review, we present some recent progress in lipid force-field development and validation strategies. In particular, we highlight how a range of properties obtained from various experimental techniques on lipid bilayers and membranes, can be used to assess the quality of a force field. We discuss the limitations and assumptions that are inherent to both computational and experimental approaches and how these can influence the comparison between simulations and experimental data. This article is part of a Special Issue entitled: Membrane Proteins edited by J.C. Gumbart and Sergei Noskov.

Understanding the accumulation of P-glycoprotein substrates within cells: The effect of cholesterol on membrane partitioning.

The apparent activity of the multidrug transporter P-glycoprotein (P-gp) is enhanced by the presence of cholesterol. Whether this is due to the direct effect of cholesterol on the activity of P-gp, its effect on the local concentration of substrate in the membrane, or its effect on the rate of entry of the drug into the cell, is unknown. In this study, molecular dynamics simulation techniques coupled with potential of mean force calculations have been used to investigate the role of cholesterol in the movement of four P-gp substrates across a POPC bilayer in the presence or absence of 10% cholesterol. The simulations suggest that the presence of cholesterol lowers the free energy associated with entering the middle of the bilayer in a substrate-specific manner. These findings suggest that P-gp substrates may preferentially accumulate in cholesterol-rich regions of the membrane, which may explain its enhanced transport activity.

Membrane-binding properties of gating modifier and pore-blocking toxins: Membrane interaction is not a prerequisite for modification of channel gating.

Many venom peptides are potent and selective inhibitors of voltage-gated ion channels, including channels that are validated therapeutic targets for treatment of a wide range of human diseases. However, the development of novel venom-peptide-based therapeutics requires an understanding of their mechanism of action. In the case of voltage-gated ion channels, venom peptides act either as pore blockers that bind to the extracellular side of the channel pore or gating modifiers that bind to one or more of the membrane-embedded voltage sensor domains. In the case of gating modifiers, it has been debated whether the peptide must partition into the membrane to reach its binding site. In this study, we used surface plasmon resonance, fluorescence spectroscopy and molecular dynamics to directly compare the lipid-binding properties of two gating modifiers (µ-TRTX-Hd1a and ProTx-I) and two pore blockers (ShK and KIIIA). Only ProTx-I was found to bind to model membranes. Our results provide further evidence that the ability to insert into the lipid bilayer is not a requirement to be a gating modifier. In addition, we characterised the surface of ProTx-I that mediates its interaction with neutral and anionic phospholipid membranes and show that it preferentially interacts with anionic lipids.

Elucidating the Spatial Arrangement of Emitter Molecules in Organic Light-Emitting Diode Films.

The effect of varying the emitter concentration on the structural properties of an archetypal phosphorescent blend consisting of 4,4'-bis(N-carbazolyl)biphenyl and tris(2-phenylpyridyl)iridium(III) has been investigated using non-equilibrium molecular dynamics (MD) simulations that mimic the process of vacuum deposition. By comparison with reflectometry measurements, we show that the simulations provide an accurate model of the average density of such films. The emitter molecules were found not to be evenly distributed throughout film, but rather they can form networks that provide charge and/or energy migration pathways, even at emitter concentrations as low as ≈5 weight percent. At slightly higher concentrations, percolated networks form that span the entire system. While such networks would give improved charge transport, they could also lead to more non-radiative pathways for the emissive state and a resultant loss of efficiency.