The successful application of high throughput molecular simulations to determine biochemical properties would be of great importance to the biomedical community if such simulations could be turned around in a clinically relevant timescale. An important example is the determination of antiretroviral inhibitor efficacy against varying strains of HIV through calculation of drug−protein binding affinities. We describe the Binding Affinity Calculator (BAC), a tool for the automated calculation of HIV-1 protease−ligand binding affinities. The tool employs fully atomistic molecular simulations alongside the well established molecular mechanics Poisson−Boltzmann solvent accessible surface area (MMPBSA) free energy methodology to enable the calculation of the binding free energy of several ligand−protease complexes, including all nine FDA approved inhibitors of HIV-1 protease and seven of the natural substrates cleaved by the protease. This enables the efficacy of these inhibitors to be ranked across several mutant strains of the protease relative to the wildtype. BAC is a tool that utilizes the power provided by a computational grid to automate all of the stages required to compute free energies of binding: model preparation, equilibration, simulation, postprocessing, and data-marshaling around the generally widely distributed compute resources utilized. Such automation enables the molecular dynamics methodology to be used in a high throughput manner not achievable by manual methods. This paper describes the architecture and workflow management of BAC and the function of each of its components. Given adequate compute resources, BAC can yield quantitative information regarding drug resistance at the molecular level within 96 h. Such a timescale is of direct clinical relevance and can assist in decision support for the assessment of patient-specific optimal drug treatment and the subsequent response to therapy for any given genotype.

Accurate calculation of important thermodynamic properties, such as macromolecular binding free energies, is one of the principal goals of molecular dynamics simulations. However, single long simulation frequently produces incorrectly converged quantitative results due to inadequate sampling of conformational space in a feasible wall-clock time. Multiple short (ensemble) simulations have been shown to explore conformational space more effectively than single long simulations, but the two methods have not yet been thermodynamically compared. Here we show that, for end-state binding free energy determination methods, ensemble simulations exhibit significantly enhanced thermodynamic sampling over single long simulations and result in accurate and converged relative binding free energies that are reproducible to within 0.5 kcal/mol. Completely correct ranking is obtained for six HIV-1 protease variants bound to lopinavir with a correlation coefficient of 0.89 and a mean relative deviation from experiment of 0.9 kcal/mol. Multidrug resistance to lopinavir is enthalpically driven and increases through a decrease in the protein−ligand van der Waals interaction, principally due to the V82A/I84V mutation, and an increase in net electrostatic repulsion due to water-mediated disruption of protein−ligand interactions in the catalytic region. Furthermore, we correctly rank, to within 1 kcal/mol of experiment, the substantially increased chemical potency of lopinavir binding to the wild-type protease compared to saquinavir and show that lopinavir takes advantage of a decreased net electrostatic repulsion to confer enhanced binding. Our approach is dependent on the combined use of petascale computing resources and on an automated simulation workflow to attain the required level of sampling and turn around time to obtain the results, which can be as little as three days. This level of performance promotes integration of such methodology with clinical decision support systems for the optimization of patient-specific therapy.

Optimization of ligand binding affinity to the target protein of interest is a primary objective in small-molecule drug discovery. Until now, the prediction of binding affinities by computational methods has not been widely applied in the drug discovery process, mainly because of its lack of accuracy and reproducibility as well as the long turnaround times required to obtain results. Herein we report on a collaborative study that compares tropomyosin receptor kinase A (TrkA) binding affinity predictions using two recently formulated fast computational approaches, namely, Enhanced Sampling of Molecular dynamics with Approximation of Continuum Solvent (ESMACS) and Thermodynamic Integration with Enhanced Sampling (TIES), to experimentally derived TrkA binding affinities for a set of Pfizer pan-Trk compounds. ESMACS gives precise and reproducible results and is applicable to highly diverse sets of compounds. It also provides detailed chemical insight into the nature of ligand–protein binding. TIES can predict and thus optimize more subtle changes in binding affinities between compounds of similar structure. Individual binding affinities were calculated in a few hours, exhibiting good correlations with the experimental data of 0.79 and 0.88 from the ESMACS and TIES approaches, respectively. The speed, level of accuracy, and precision of the calculations are such that the affinity predictions can be used to rapidly explain the effects of compound modifications on TrkA binding affinity. The methods could therefore be used as tools to guide lead optimization efforts across multiple prospective structurally enabled programs in the drug discovery setting for a wide range of compounds and targets.

The accurate prediction of the binding affinities of ligands to proteins is a major goal in drug discovery and personalized medicine. The time taken to make such predictions is of similar importance to their accuracy, precision, and reliability. In the past few years, an ensemble based molecular dynamics approach has been proposed that provides a route to reliable predictions of free energies based on the molecular mechanics Poisson–Boltzmann surface area method which meets the requirements of speed, accuracy, precision, and reliability. Here, we describe an equivalent methodology based on thermodynamic integration to substantially improve the speed, accuracy, precision, and reliability of calculated relative binding free energies. We report the performance of the method when applied to a diverse set of protein targets and ligands. The results are in very good agreement with experimental data (90% of calculations agree to within 1 kcal/mol), while the method is reproducible by construction. Statistical uncertainties of the order of 0.5 kcal/mol or less are achieved. We present a systematic account of how the uncertainty in the predictions may be estimated.

Binding free energies of bromodomain inhibitors are calculated with recently formulated approaches, namely ESMACS (enhanced sampling of molecular dynamics with approximation of continuum solvent) and TIES (thermodynamic integration with enhanced sampling). A set of compounds is provided by GlaxoSmithKline, which represents a range of chemical functionality and binding affinities. The predicted binding free energies exhibit a good Spearman correlation of 0.78 with the experimental data from the 3-trajectory ESMACS, and an excellent correlation of 0.92 from the TIES approach where applicable. Given access to suitable high end computing resources and a high degree of automation, we can compute individual binding affinities in a few hours with precisions no greater than 0.2 kcal/mol for TIES, and no larger than 0.34 and 1.71 kcal/mol for the 1- and 3-trajectory ESMACS approaches.

We introduce the lattice-Boltzmann code LB3D, version 7.1. Building on a parallel program and supporting tools which have enabled research utilising high performance computing resources for nearly two decades, LB3D version 7 provides a subset of the research code functionality as an open source project. Here, we describe the theoretical basis of the algorithm as well as computational aspects of the implementation. The software package is validated against simulations of meso-phases resulting from self-assembly in ternary fluid mixtures comprising immiscible and amphiphilic components such as water–oil–surfactant systems. The impact of the surfactant species on the dynamics of spinodal decomposition are tested and quantitative measurement of the permeability of a body centred cubic (BCC) model porous medium for a simple binary mixture is described. Single-core performance and scaling behaviour of the code are reported for simulations on current supercomputer architectures.Program summaryProgram Title: LB3DProgram Files doi: http://dx.doi.org/10.17632/9g9x2wr8z8.1 Licensing provisions: BSD 3-clauseProgramming language: FORTRAN90, Python, CNature of problem: Solution of the hydrodynamics of single phase, binary immiscible and ternary amphiphilic fluids. Simulation of fluid mixtures comprising miscible and immiscible fluid components as well as amphiphilic species on the mesoscopic scale. Observable phenomena include self-organisation of mesoscopic complex fluid phases and fluid transport in porous media.Solution method: Lattice-Boltzmann (lattice-Bhatnagar–Gross–Krook, LBGK) [1, 2, 3] method describing fluid dynamics in terms of the single particle velocity distribution function in a 3-dimensional discrete phase space (D3Q19) [4, 5, 6]. Multiphase interactions are modelled using a phenomenological pseudo-potential approach [7, 8] with amphiphilic interactions utilising an additional dipole field [9, 10]. Solid boundaries are modelled using simple bounce-back boundary conditions and additional pseudo-potential wetting interactions [11].Additional comments including Restrictions and Unusual features: The purpose of the release is the provision of a refactored minimal version of LB3D suitable as a starting point for the integration of additional features building on the parallel computation and IO functionality. [1]S. Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond, Oxford University Press, 2001.[2]B. Dünweg, A. Ladd, Lattice Boltzmann simulations of soft matter systems, Adv. Poly. Sci. 221 (2009) 89–166[3]C. K. Aidun, J. R. Clausen, Lattice-Boltzmann Method for Complex Flows, Annual Review of Fluid Mechanics 42 (2010) 439.[4]X. He, L.-S. Luo, A priori derivation of the lattice-Boltzmann equation, Phys. Rev. E 55 (1997) R6333.[5]X. He, L.-S. Luo, Theory of the lattice Boltzmann method: from the Boltzmann equation to the lattice Boltzmann equation, Phys. Rev. E 56.[6]Y. H. Qian, D. D’Humiéres, P. Lallemand, Lattice BGK Models for Navier–Stokes Equation, Europhysics Letters 17 (1992) 479.[7]X. Shan, H. Chen, Lattice-Boltzmann model for simulating flows with multiple phases and components, Physical Review E 47 (1993) 1815.[8]X. Shan, G. Doolen, Multicomponent lattice-Boltzmann model with interparticle interaction, Journal of Statistical Physics 81 (1995) 379.[9]H. Chen, B. Boghosian, P.V. Coveney, M. Nekovee, A ternary lattice-Boltzmann model for amphiphilic fluids, Proceedings of the Royal Society of London A 456 (2000) 2043.[10]M. Nekovee, P. V. Coveney, H. Chen, B. M. Boghosian, Lattice-Boltzmann model for interacting amphiphilic fluids, Phys. Rev. E 62 (2000) 8282.[11]N. S. Martys, H. Chen, Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice-Boltzmann method, Phys. Rev. E 53 (1996) 743.

We present FabSim, a toolkit developed to simplify a range of computational tasks for researchers in diverse disciplines. FabSim is flexible, adaptable, and allows users to perform a wide range of tasks with ease. It also provides a systematic way to automate the use of resources, including HPC and distributed machines, and to make tasks easier to repeat by recording contextual information. To demonstrate this, we present three use cases where FabSim has enhanced our research productivity. These include simulating cerebrovascular bloodflow, modelling clay-polymer nanocomposites across multiple scales, and calculating ligand–protein binding affinities.Program summaryProgram title: FabSimCatalogue identifier: AFAO_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AFAO_v1_0.htmlProgram obtainable from: CPC Programme Library, Queen’s University, Belfast, N. IrelandLicensing provisions: BSD 3-ClauseNo. of lines in distributed program, including test data, etc.: 268282No. of bytes in distributed program, including test data, etc.: 2791792Distribution format: tar.gzProgramming language: Python.Computer: PC or Mac.Operating system: Unix, OSX.RAM: 1 GbytesClassification: 3, 4, 6.5.External routines: NumPy, SciPy, Fabric (1.5 or newer), PyYamlNature of problem:Running advanced computations using remote resources is an activity that requires considerable time and human attention. These activities, such as organizing data, configuring software and setting up individual runs often vary slightly each time they are performed. To lighten this burden, we required an approach that introduced little burden of its own to set up and adapt, beyond which very substantial productivity ensues.Solution method:We present a toolkit which helps to simplify and automate the activities which surround computational science research. FabSim is aimed squarely at the experienced computational scientist, who can use the command line interface and a system of modifiable content to quickly automate sets of research tasks.Restrictions:FabSim relies on a command-line interface, and assumes some level of scripting knowledge from the user.Unusual features:FabSim has a proven track record of being easy to adapt. It has already been extensively adapted to facilitate leading research in the modelling of bloodflow, nanomaterials, and ligand–protein binding.Running time:FabSim can be used interactively, typically requiring a few seconds to perform a basic task.

We describe the mechanism that leads to full exfoliation and dispersion of organophilic clays when mixed with molten hydrophilic polymers. This process is of fundamental importance for the production of clay–polymer nanocomposites with enhanced materials properties. The chemically specific nature of our multiscale approach allows us to probe how chemistry, in combination with processing conditions, produces such materials properties at the mesoscale and beyond. In general agreement with experimental observations, we find that a higher grafting density of charged quaternary ammonium surfactant ions promotes exfoliation, by a mechanism whereby the clay sheets slide transversally over one another. We can determine the elastic properties of these nanocomposites; exfoliated and partially exfoliated morphologies lead to substantial enhancement of the Young’s modulus, as found experimentally.

The presentation of potentially pathogenic peptides by major histocompatibility complex (MHC) molecules is one of the most important processes in adaptive immune defense. Prediction of peptide–MHC (pMHC) binding affinities is therefore a principal objective of theoretical immunology. Machine learning techniques achieve good results if substantial experimental training data are available. Approaches based on structural information become necessary if sufficiently similar training data are unavailable for a specific MHC allele, although they have often been deemed to lack accuracy. In this study, we use a free energy method to rank the binding affinities of 12 diverse peptides bound by a class I MHC molecule HLA-A*02:01. The method is based on enhanced sampling of molecular dynamics calculations in combination with a continuum solvent approximation and includes estimates of the configurational entropy based on either a one or a three trajectory protocol. It produces precise and reproducible free energy estimates which correlate well with experimental measurements. If the results are combined with an amino acid hydrophobicity scale, then an extremely good ranking of peptide binding affinities emerges. Our approach is rapid, robust, and applicable to a wide range of ligand–receptor interactions without further adjustment.

Near attack conformations (NACs) are conformations extending from the ground state (GS) that lie on the transition path of a chemical reaction. Here, we develop a method for computing the thermodynamic contribution to catalysis due to NAC formation in bimolecular reactions, within the limit of a classical molecular dynamics force field. We make use of the Bürgi–Dunitz theory applied to large-scale unbiased all-atom ensemble molecular dynamics simulations. We apply this to HIV-1 protease peptide hydrolysis, known to achieve a rate enhancement of ∼1011 (ΔGcat⧧ ∼ 15 kcal/mol) over the uncatalyzed bimolecular reaction (ΔGnon⧧ ∼ 30 kcal/mol). The ground state consists of a nucleophilic water molecule bound to an octapeptide substrate in the active site. We first observe multiple and reversible binding of a nucleophilic water molecule into the active site giving a free energy of binding of ΔG = −1 kcal/mol to form the GS. The free energy barriers for catalyzed and uncatalyzed NAC formation are both equivalent: ΔGNAC⧧ = 4.6 kcal/mol, constituting ∼30% and ∼15% of the overall barriers, respectively. Therefore, not only does adoption of NACs only account for minor progress along the transition path in both catalyzed and uncatalyzed reactions, but there is no preferential formation of them in the catalyzed reaction. Analysis of the catalytic hydrogen bond network reveals interactions that stabilize the GS; however, subsequent NAC formation does not preferentially favor any of the possible hydrogen bond configurations. This supports the view that the catalytic power of HIV-1 protease is not due to NAC formation.

Compared with proteins, the relationship between structure, dynamics, and function of RNA enzymes (known as ribozymes) is far less well understood, despite the fact that ribozymes are found in many organisms and are often conceived as “molecular fossils” of the first self-replicating molecules to have arisen on Earth. To investigate how ribozymal function is governed by structure and dynamics, we study the full hammerhead ribozyme in bulk water and in an aqueous clay mineral environment by computer simulation using replica-exchange molecular dynamics. Through extensive sampling of the major conformational states of the hammerhead ribozyme, we are able to show that the hammerhead manifests a free-energy landscape reminiscent of that which is well known in proteins, exhibiting a “funnel” topology that guides the ribozyme into its globally most stable conformation. The active-site geometry is found to be closely correlated to the tertiary structure of the ribozyme, thereby reconciling conflicts between previously proposed mechanisms for the self-scission of the hammerhead. The conformational analysis also accounts for the differences reported experimentally in the catalytic activity of the hammerhead ribozyme, which is reduced when interacting with clay minerals as compared with bulk water.

The use of molecular simulation to estimate the strength of macromolecular binding free energies is becoming increasingly widespread, with goals ranging from lead optimization and enrichment in drug discovery to personalizing or stratifying treatment regimes. In order to realize the potential of such approaches to predict new results, not merely to explain previous experimental findings, it is necessary that the methods used are reliable and accurate, and that their limitations are thoroughly understood. However, the computational cost of atomistic simulation techniques such as molecular dynamics (MD) has meant that until recently little work has focused on validating and verifying the available free energy methodologies, with the consequence that many of the results published in the literature are not reproducible. Here, we present a detailed analysis of two of the most popular approximate methods for calculating binding free energies from molecular simulations, molecular mechanics Poisson–Boltzmann surface area (MMPBSA) and molecular mechanics generalized Born surface area (MMGBSA), applied to the nine FDA-approved HIV-1 protease inhibitors. Our results show that the values obtained from replica simulations of the same protease–drug complex, differing only in initially assigned atom velocities, can vary by as much as 10 kcal mol–1, which is greater than the difference between the best and worst binding inhibitors under investigation. Despite this, analysis of ensembles of simulations producing 50 trajectories of 4 ns duration leads to well converged free energy estimates. For seven inhibitors, we find that with correctly converged normal mode estimates of the configurational entropy, we can correctly distinguish inhibitors in agreement with experimental data for both the MMPBSA and MMGBSA methods and thus have the ability to rank the efficacy of binding of this selection of drugs to the protease (no account is made for free energy penalties associated with protein distortion leading to the over estimation of the binding strength of the two largest inhibitors ritonavir and atazanavir). We obtain improved rankings and estimates of the relative binding strengths of the drugs by using a novel combination of MMPBSA/MMGBSA with normal mode entropy estimates and the free energy of association calculated directly from simulation trajectories. Our work provides a thorough assessment of what is required to produce converged and hence reliable free energies for protein–ligand binding.

The human histone deacetylase 8 (HDAC8) is a key hydrolase in gene regulation and has been identified as a drug target for the treatment of several cancers. Previously the HDAC8 enzyme has been extensively studied using biochemical techniques, X-ray crystallography, and computational methods. Those investigations have yielded detailed information about the active site and have demonstrated that the substrate entrance surface is highly dynamic. Yet it has remained unclear how the dynamics of the entrance surface tune and influence the catalytic activity of HDAC8. Using long time scale all atom molecular dynamics simulations we have found a mechanism whereby the interactions and dynamics of two loops tune the configuration of functionally important residues of HDAC8 and could therefore influence the activity of the enzyme. We subsequently investigated this hypothesis using a well-established fluorescence activity assay and a noninvasive real-time progression assay, where deacetylation of a p53 based peptide was observed by nuclear magnetic resonance spectroscopy. Our work delivers detailed insight into the dynamic loop network of HDAC8 and provides an explanation for a number of experimental observations.

We present the results of large-scale molecular simulations, run over several tens of nanoseconds, of 25-mer sequences of single-stranded ribonucleic acid (RNA) in bulk water and at the surface of three hydrated positively charged MgAl layered double hydroxide (LDH) minerals. The three LDHs differ in surface charge density, through varying the number of isomorphic Al substitutions. Over the course of the simulations, RNA adsorbs tightly to the LDH surface through electrostatic interactions between the charged RNA phosphate groups and the alumina charge sites present in the LDH sheet. The RNA strands arrange parallel to the surface with the base groups aligning normal to the surface and exposed to the bulk aqueous region. This templating effect makes LDH a candidate for amplifying the population of a known RNA sequence from a small number of RNAs. The structure and interactions of RNA at a positively charged, hydroxylated LDH surface were compared with those of RNA at a positively charged calcium montmorillonite surface, allowing us to establish the comparative effect of complexation and water structure at hydroxide and silicate surfaces. The systems were studied by computing radial distribution functions, atom density plots, and radii of gyration, as well as visualization. An observation pertinent to the role of these minerals in prebiotic chemistry is that, for a given charge density on the mineral surface, different genetic sequences of RNA adopt different configurations.

Origins of life studies represent an exciting and highly multidisciplinary research field. In this review we focus on the contributions made by theory, modelling and simulation to addressing fundamental issues in the domain and the advances these approaches have helped to make in the field. Theoretical approaches will continue to make a major impact at the “systems chemistry” level based on the analysis of the remarkable properties of nonlinear catalytic chemical reaction networks, which arise due to the auto-catalytic and cross-catalytic nature of so many of the putative processes associated with self-replication and self-reproduction. In this way, we describe inter alia nonlinear kinetic models of RNA replication within a primordial Darwinian soup, the origins of homochirality and homochiral polymerization. We then discuss state-of-the-art computationally-based molecular modelling techniques that are currently being deployed to investigate various scenarios relevant to the origins of life.

One of the principal targets in human immunodeficiency virus type 1 (HIV-1) therapy is the reverse transcriptase (RT) enzyme. Non-nucleoside RT inhibitors (NNRTIs) are a class of highly specific drugs which bind to a pocket approximately 10 Å from the polymerase active site, inhibiting the enzyme allosterically. It is widely believed that NNRTIs function as “molecular wedges”, disrupting the region between thumb and palm subdomains of the p66 subunit and locking the thumb in a wide-open conformation. Crystal structure data suggest that the binding of NNRTIs forces RT into a wide-open conformation in which the separation between the thumb and fingers subdomains is much higher than in the apo structure. Using ensemble molecular dynamics simulations (aggregate sampling ∼600 ns), we have captured RT bound to the NNRTI efavirenz in a closed conformation similar to that of the apo enzyme, suggesting the constraint of thumb motion is not as complete as previously believed. Rather, our investigation confirms that a conformational distribution across open and closed states must exist in the drug-bound enzyme and that allosteric modulation is effected via the alteration of the kinetic landscape of conformational transitions upon drug-binding. A more detailed understanding of the mechanism of NNRTI inhibition and the effect of binding upon domain motion could aid the design of more effective inhibitors and help identify novel allosteric sites.

Janus kinase 2 (JAK2) is a protein tyrosine kinase implicated in signaling by specific members of the cytokine receptor family. Although it has been established that the JAK2 tyrosine kinase is negatively regulated by the JAK homology 2 (JH2) pseudokinase domain, the underlying mechanism of JH2 mediated regulation remains elusive. To elucidate the regulation of JAK2 kinase, we have built a structural model for the kinase and pseudokinase domains of JAK2. An asymmetric dimer is proposed, in which the kinase domain JH1 occupies a position where it could not be activated. We investigate the dynamic and energetic properties of the dimer by molecular dynamics simulation. JAK2 activation requires the two domains to be dissociated and rearranged in a form such that the JH1 kinase domain can adopt an active conformation. The significance of the above mechanism is emphasized by the finding that the activating V617F mutation destabilizes JH1–JH2 association in the proposed asymmetric dimer. Thus abrogation of the domain–domain interaction seems to be a possible first step for the structural rearrangement of the two domains, resulting in constitutive activation of JAK2 by the V617F mutation.

The development of resistance to different drugs remains a major problem for a wide range of infections. In particular, combinations of specific mutations, which individually demonstrate no effect, exhibit significant cooperativity. Here we show that changes to the energy of ligand binding in different resistant HIV-1 proteases are correlated with the creation of water binding sites in the active site. This correlation is conserved across two drugs (ritonavir and lopinavir). We propose that individual mutations induce changes in flap packing that are insufficient to allow water binding but in combination allow access, leading to the observed cooperative resistance.

The emergence of drug resistance is a major challenge for the effective treatment of HIV. In this article, we explore the application of atomistic molecular dynamics simulations to quantify the level of resistance of a patient-derived HIV-1 protease sequence to the inhibitor lopinavir. A comparative drug ranking methodology was developed to compare drug resistance rankings produced by the Stanford HIVdb, ANRS, and RegaDB clinical decision support systems. The methodology was used to identify a patient sequence for which the three rival online tools produced differing resistance rankings. Mutations at only three positions (L10I, A71IV, and L90M) influenced the resistance level assigned to the sequence. We use ensemble molecular dynamics simulations to elucidate the origin of these discrepancies and the mechanism of resistance. By simulating not only the full patient sequences but also systems containing the constituent mutations, we gain insight into why resistance estimates vary and the interactions between the various mutations. In the same way, we also gain valuable knowledge of the mechanistic causes of resistance. In particular, we identify changes in the relative conformation of the two beta sheets that form the protease dimer interface which suggest an explanation of the relative frequency of different amino acids observed in patients at residue 71.

Since a mineral-mediated origin of life was first hypothesized over 60 years ago, clays have played a significant role in origins of life studies. Such studies have hitherto rarely used computer simulation to understand the possible chemical pathways to the formation of biomolecules. We use molecular dynamics techniques, performed on supercomputing grids, to carry out large-scale simulations of various 25-mer sequences of ribonucleic acid (RNA), in bulk water and with aqueous montmorillonite clay over many tens of nanoseconds. Hitherto, there has only been limited experimental data reported for these systems. Our simulations are found to be in agreement with various experimental observations pertaining to the relative adsorption of RNA on montmorillonite in the presence of charge balancing cations. Over time scales of only a few nanoseconds, specific RNA sequences fold to characteristic secondary structural motifs, which do not form in the corresponding bulk water simulations. Our simulations also show that, in aqueous Ca2+ environments, RNA can tether to the clay surface through a nucleotide base, leaving the 3′-end of the strand exposed, providing a mechanism for the regiospecific adsorption and elongation of RNA oligomers on clay surfaces.

We review the recent advances in large-scale and coarse-grained molecular dynamics applied to clay minerals. Recent advances in local and distributed high performance computational resources together with the development of efficient parallelized algorithms has enabled the simulation of increasingly realistic large-scale models of clay mineral systems. Using this improved technology, it is becoming possible to simulate realistic clay platelet sizes at an atomistic level. This has considerably extended the spatial dimensions of microscopic simulation into a domain normally encountered in mesoscopic simulation. The simulation of large-scale model systems is important to further study complex phenomena, such as the structural and mechanical properties of disordered layered materials such as clays. In order to achieve even larger length and longer time-scales coarse-grained methods are increasingly employed, capturing phenomena such as composite failure modes and intercalation.

The translocation of polynucleotides through transmembrane protein pores is a fundamental biological process with important technological and medical relevance. The translocation process is complex, and it is influenced by a range of factors including the diameter and inner surface of the pore, the secondary structure of the polymer, and the interactions between the polymer and protein. In this paper, we perform nonequilibrium constant velocity-steered molecular dynamics simulations of nucleic acid molecule translocation through the protein nanopore α-hemolysin and use Jarzynski’s identity to determine the associated free energy profiles. With this approach we are able to explain the observed differences in experimental translocation time through the nanopore between polyadenosine and polydeoxycytidine. The translocation of polynucleotides and single nucleotides through α-hemolysin is investigated. These simulations are computationally intensive as they employ models with atomistic level resolution; in addition to their size, these systems are challenging to study due to the time scales of translocation of large asymmetric molecules. Our simulations provide insight into the role of the interactions between the nucleic acid molecules and the protein pore. Mutated protein pores provide confirmation of residue-specific interactions between nucleotides and the protein pore. By harnessing such molecular dynamics simulations, we gain new physicochemical insight into the translocation process.

We investigate the rheological characteristics of ternary amphiphilic gyroid, diamond and primitive cubic phases under applied Couette flow simulated using a kinetic lattice–Boltzmann model and periodic Lees–Edwards boundary conditions. The simulated rheological response of the cubic phases is compared to experimental observations in lyotropic liquid crystals. We relate the variations in the non-Newtonian response and deformation under strain in these cubic phases to their triply bicontinuous cubic morphologies as well as to the differences in the interaction parameters between the three species present in the amphiphilic system. The large system sizes allow simulation of multiple domains which elucidate the correlation between the evolution of the defect texture and the change in the stress field of the cubic phase under applied Couette flow.

Very large-scale molecular dynamics simulations are performed to investigate the effects of montmorillonite clay filler on poly(ethylene) glycol in the formation of clay-polymer nanocomposites. We present the results of MD simulations of intercalated and exfoliated nanocomposites at sizes which approach those of a realistic clay platelet. The simulations allow us to determine the difference between polymer adsorbed on the surface of the clay and that more remote from it. All polymers arrange themselves in layers parallel to the surface, each layer being approximately 4 Å thick. We find the polymer conformation of the inner-most layer is distinct, due to complexation with the counterions found near the charged clay surface. The diffusion of polymer within this layer is much lower than that of others layers, which effectively increases the size of the nanofiller and makes gas permeation more tortuous. We perform non-equilibrium molecular dynamics simulations by imposing a strain on the model and analysing the stress response. By partitioning the stress response into clay and different polymer layers, we find that the inner-most polymer layer has a much higher Young's modulus than the remaining layers in the direction of the polymer chains.

Clay-polymer nanocomposites are a new range of particle-filled composites possessing enhanced mechanical properties. Transmission electron microscopy studies reveal that many single and double clay sheets are bent when immersed in the polymer matrix. To understand the elastic properties of the sheets we investigate their behaviour under compression through molecular dynamics simulations in both aqueous and poly(ethylene oxide) environments. We examine the sheets before and after buckling and find that the elastic properties of the clay sheets are much reduced when buckled in each case. Buckling occurs with only small out-of-plane displacement and we ascertain that any observably bent sheet will have reduced elastic properties. The sheets possess a negative Poisson ratio during buckling, which allows the clay sheet to regain its uncompressed area. We present a simple analytical model which explains the buckling through the competition between in-plane strain energy and out-of-plane displacements. Knowledge of how the elastic properties of the clay sheets change with shape is important for the design of new nanocomposites.

Grid computing is distributed computing performed transparently across multiple administrative domains. Grid middleware, which is meant to enable access to grid resources, is currently widely seen as being too heavyweight and, in consequence, unwieldy for general scientific use. Its heavyweight nature, especially on the client-side, has severely restricted the uptake of grid technology by computational scientists. In this paper, we describe the Application Hosting Environment (AHE) which we have developed to address some of these problems. The AHE is a lightweight, easily deployable environment designed to allow the scientist to quickly and easily run legacy applications on distributed grid resources. It provides a higher level abstraction of a grid than is offered by existing grid middleware schemes such as the Globus Toolkit. As a result, the computational scientist does not need to know the details of any particular underlying grid middleware and is isolated from any changes to it on the distributed resources. The functionality provided by the AHE is ‘application-centric’: applications are exposed as web services with a well-defined standards-compliant interface. This allows the computational scientist to start and manage application instances on a grid in a transparent manner, thus greatly simplifying the user experience. We describe how a range of computational science codes have been hosted within the AHE and how the design of the AHE allows us to implement complex workflows for deployment on grid infrastructure.

The use of computational methods for the study of clay minerals has become an essential adjunct to experimental techniques for the analysis of these poorly ordered materials. Although information may be obtained through conventional methods of analysis regarding macroscopic properties of clay minerals, information about the spatial arrangement of molecules within the interlayers is hard to obtain without the aid of computer simulation. The interpretation of experimental data from techniques such as solid-state nuclear magnetic resonance or neutron diffraction studies is considerably assisted by the application of computer simulations. Using a series of case studies, we review the techniques, applications and insight gained from the use of molecular simulation applied to the study of clay systems (particularly for materials applications). The amount of information that can be gleaned from such simulations continues to grow, and is leading to ever larger-scale and hence more realistic classical and quantum mechanical studies which promise to reveal new and unexpected phenomena.

Coupled models are set to become increasingly important in all aspects of science and engineering as tools with which to study complex systems in an integrated manner. Such coupled, hybrid simulations typically communicate data between the component models of which they are comprised relatively infrequently, and so a Grid is expected to present an ideal architecture on which to run them. In the present paper, we describe a simple, flexible and extensible architecture for a two-component hybrid molecular-continuum coupled model (hybrid MD). We discuss its deployment on distributed resources and the extensions to the RealityGrid computational-steering system to handle coupled models.

Certain systems, such as amphiphile solutions or diblock copolymer melts, may assemble into structures called “mesophases”, with properties intermediate between those of a solid and a liquid. These mesophases can be of very regular structure, but may contain defects and grain boundaries. Different visualization techniques such as volume rendering or isosurfacing of fluid density distributions allow the human eye to detect and track defects in liquid crystals because humans are easily capable of finding imperfections in repetitive spatial structures. However, manual data analysis becomes too time consuming and algorithmic approaches are needed when there are large amounts of data. We present and compare two different approaches we have developed to study defects in gyroid mesophases of amphiphilic ternary fluids. While the first method is based on a pattern recognition algorithm, the second uses the particular structural properties of gyroid mesophases to detect defects.

We report a theoretical investigation of hydrated clay–polymer nanocomposites exchanged with Li+, Na+ and K+. This work is the result of the implementation of Teppen's force field within a highly scalable molecular dynamics (MD) program called Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) with which, we performed large-scale MD simulations. The results show that, in contradiction to the situation pertaining in the absence of polymer, the behaviour of Li+ and Na+ based nanocomposites is quite distinct. Unlike K+ and Na+, the Li+ cations are observed to diffuse within the tetrahedral pockets of the clay sheets as well as the centre of the galleries.

This paper reports a study of the intercalation of single bifunctional poly(ethylene oxide)- and poly(propylene oxide)-based compounds, with and without modified end groups, in Na+-montmorillonite and Na+-hectorite clays as the starting point for the design of new non-exfoliated clay–polymer materials. The terminal functional groups employed on intercalated poly(ethylene oxide) and poly(propylene oxide) include acrylate, methacrylate and amine functions. For some of the monomers, in situ polymerisation was observed. Attempts to rationalise these behaviours are made on the basis of atomistic grand canonical Monte Carlo, canonical Monte Carlo and molecular dynamics simulations. Poly(ethylene oxide) diacrylate monomers, for which in situ polymerisation was observed, interact more strongly with the Na+-montmorillonite clay layers than the other monomers.

We compare two recently developed mesoscale models of binary immiscible and ternary amphiphilic fluids. We describe and compare the algorithms in detail and discuss their stability properties. The simulation results for the cases of self-assembly of ternary droplet phases and binary water-amphiphile sponge phases are compared and discussed. Both models require parallel implementation and deployment on large scale parallel computing resources in order to achieve reasonable simulation times for three-dimensional models. The parallelization strategies and performance on two distinct parallel architectures are compared and discussed. Large scale three-dimensional simulation of multiphase fluids requires the extensive use of high performance visualization techniques in order to enable the large quantities of complex data to be interpreted. We report on our experiences with two commercial visualization products: AVS and VTK. We also discuss the application and use of novel computational steering techniques for the more efficient utilization of high performance computing resources. We close the paper with some suggestions for the future development of both models.

Highlights•MUSCLE 2 simulates multiscale phenomena by coupling heterogeneous submodels.•MUSCLE 2 supports Java, C/C++, Fortran and Python code with MPI, OpenMP and threads.•The total coupling overhead of MUSCLE 2 is in the order of seconds.•The message speed of MUSCLE 2 beats that of file copy and GridFTP but not that of MPI.•A canal system model shows an acceptable overhead of MUSCLE 2 for large problem sizes.We present the Multiscale Coupling Library and Environment: MUSCLE 2. This multiscale component-based execution environment has a simple to use Java, C++, C, Python and Fortran API, compatible with MPI, OpenMP and threading codes. We demonstrate its local and distributed computing capabilities and compare its performance to MUSCLE 1, file copy, MPI, MPWide, and GridFTP. The local throughput of MPI is about two times higher, so very tightly coupled code should use MPI as a single submodel of MUSCLE 2; the distributed performance of GridFTP is lower, especially for small messages. We test the performance of a canal system model with MUSCLE 2, where it introduces an overhead as small as 5% compared to MPI.

Highlights•MUSCLE 2 simulates multiscale phenomena by coupling heterogeneous submodels.•MUSCLE 2 supports Java, C/C++, Fortran and Python code with MPI, OpenMP and threads.•The total coupling overhead of MUSCLE 2 is in the order of seconds.•The message speed of MUSCLE 2 beats that of file copy and GridFTP but not that of MPI.•A canal system model shows an acceptable overhead of MUSCLE 2 for large problem sizes.We present the Multiscale Coupling Library and Environment: MUSCLE 2. This multiscale component-based execution environment has a simple to use Java, C++, C, Python and Fortran API, compatible with MPI, OpenMP and threading codes. We demonstrate its local and distributed computing capabilities and compare its performance to MUSCLE 1, file copy, MPI, MPWide, and GridFTP. The local throughput of MPI is about two times higher, so very tightly coupled code should use MPI as a single submodel of MUSCLE 2; the distributed performance of GridFTP is lower, especially for small messages. We test the performance of a canal system model with MUSCLE 2, where it introduces an overhead as small as 5% compared to MPI.

Secure access to patient data and analysis tools to run on that data will revolutionize the treatment of a wide range of diseases, by using advanced simulation techniques to underpin the clinical decision making process. To achieve these goals, suitable e-Science infrastructures are required to allow clinicians and researchers to trivially access data and launch simulations. In this paper we describe the open source Individualized MEdiciNe Simulation Environment (IMENSE), which provides a platform to securely manage clinical data, and to perform wide ranging analysis on that data, ultimately with the intention of enhancing clinical decision making with direct impact on patient health care. We motivate the design decisions taken in the development of the IMENSE system by considering the needs of researchers in the ContraCancrum project, which provides a paradigmatic case in which clinicians and researchers require coordinated access to data and simulation tools. We show how the modular nature of the IMENSE system makes it applicable to a wide range of biomedical computing scenarios, from within a single hospital to major international research projects.Highlights► Clinicians face making decisions based on analysis of many different types of data. ► Integrating disparate data sources can assist the clinical decision making process. ► The IMENSE system integrates different types of clinical data from multiple sources. ► Clinicians and researchers can use the system to perform complex analysis techniques.

We report on a sophisticated numerical study of a parallel space–time algorithm for the computation of periodic solutions of the driven, incompressible Navier–Stokes equations in the turbulent regime. Efforts to apply the machinery of dynamical systems theory to fluid turbulence depend on the ability to accurately and reliably compute such unstable periodic orbits (UPOs). For example, the UPOs may be used to construct the dynamical zeta function of the system, from which very accurate turbulent averages of observables may be extracted.Though a number of algorithms for computing such orbits have been proposed and tested, in this paper we focus on a space–time variational principle introduced by Lan and Cvitanović in 2004 [15]. This method has not, to our knowledge, been tested on dynamical systems of high dimension because of the formidable storage and computation required. In this paper, we use petascale computation to apply this algorithm to weak hydrodynamic turbulence.We begin with a brief description and reformulation of the space–time algorithm of Lan and Cvitanović. We then describe how to apply this algorithm to the lattice-Boltzmann method for the solution of the Navier–Stokes equations. In particular, we describe the fully parallel implementation of this algorithm using the Message Passing Interface. This implementation, called HYPO4D, has been successfully deployed on a large variety of platforms both in the UK and the USA and has shown very good scalability to tens of thousands of computing cores.Finally, we describe the application of this implementation to the problem of weak homogeneous turbulence driven by an Arnold–Beltrami–Childress force field in three spatial dimensions, at a Reynolds number of 371. We commence by systematically searching for nearly periodic orbits as candidate solutions from which to begin the relaxation; we then apply the variational algorithm until convergence is obtained. Because the algorithm requires storage of the space–time lattice, even the smallest orbits require resources on the order of tens of thousands of computing cores. Using this approach, two UPOs have been identified and some of their properties have been analysed.

T-cell activation requires interaction of T-cell receptors (TCR) with peptide epitopes bound by major histocompatibility complex (MHC) proteins. This interaction occurs at a special cell–cell junction known as the immune or immunological synapse. Fluorescence microscopy has shown that the interplay among one agonist peptide-MHC (pMHC), one TCR and one CD4 provides the minimum complexity needed to trigger transient calcium signalling. We describe a computational approach to the study of the immune synapse. Using molecular dynamics simulation, we report here on a study of the smallest viable model, a TCR–pMHC–CD4 complex in a membrane environment. The computed structural and thermodynamic properties are in fair agreement with experiment. A number of biomolecules participate in the formation of the immunological synapse. Multi-scale molecular dynamics simulations may be the best opportunity we have to reach a full understanding of this remarkable supra-macromolecular event at a cell–cell junction.

Highlights•We have developed the Application Interaction Model, to improve user experience when accessing e-infrastructure.•We present out implementation of this model in version 3.0 of the Application Hosting Environment.•We compare the performance of AHE3 and AHE2 to illustrate the benefits of using the new version of the software.•Finally we describe the projects and users of AHE3.Computer simulation is finding a role in an increasing number of scientific disciplines, concomitant with the rise in available computing power. Marshalling this power facilitates new, more effective and different research than has been hitherto possible. Realizing this inevitably requires access to computational power beyond the desktop, making use of clusters, supercomputers, data repositories, networks and distributed aggregations of these resources. The use of diverse e-infrastructure brings with it the ability to perform distributed multiscale simulations. Accessing one such resource entails a number of usability and security problems; when multiple geographically distributed resources are involved, the difficulty is compounded. In this paper we present a solution, the Application Hosting Environment,3 which provides a Software as a Service layer on top of distributed e-infrastructure resources. We describe the performance and usability enhancements present in AHE version 3, and show how these have led to a high performance, easy to use gateway for computational scientists working in diverse application domains, from computational physics and chemistry, materials science to biology and biomedicine.

We review the recent advances in large-scale and coarse-grained molecular dynamics applied to clay minerals. Recent advances in local and distributed high performance computational resources together with the development of efficient parallelized algorithms has enabled the simulation of increasingly realistic large-scale models of clay mineral systems. Using this improved technology, it is becoming possible to simulate realistic clay platelet sizes at an atomistic level. This has considerably extended the spatial dimensions of microscopic simulation into a domain normally encountered in mesoscopic simulation. The simulation of large-scale model systems is important to further study complex phenomena, such as the structural and mechanical properties of disordered layered materials such as clays. In order to achieve even larger length and longer time-scales coarse-grained methods are increasingly employed, capturing phenomena such as composite failure modes and intercalation.

Origins of life studies represent an exciting and highly multidisciplinary research field. In this review we focus on the contributions made by theory, modelling and simulation to addressing fundamental issues in the domain and the advances these approaches have helped to make in the field. Theoretical approaches will continue to make a major impact at the “systems chemistry” level based on the analysis of the remarkable properties of nonlinear catalytic chemical reaction networks, which arise due to the auto-catalytic and cross-catalytic nature of so many of the putative processes associated with self-replication and self-reproduction. In this way, we describe inter alia nonlinear kinetic models of RNA replication within a primordial Darwinian soup, the origins of homochirality and homochiral polymerization. We then discuss state-of-the-art computationally-based molecular modelling techniques that are currently being deployed to investigate various scenarios relevant to the origins of life.

The purpose of statistical mechanics is to provide a route to the calculation of macroscopic properties of matter from their constituent microscopic components. It is well known that the macrostates emerge as ensemble averages of microstates. However, this is more often stated than implemented in computer simulation studies. Here we consider foundational aspects of statistical mechanics which are overlooked in most textbooks and research articles that purport to compute macroscopic behaviour from microscopic descriptions based on classical mechanics and show how due attention to these issues leads in directions which have not been widely appreciated in the field of molecular dynamics simulation.