Coarse-Grained Simulation of mRNA-Loaded Lipid Nanoparticle Self-Assembly.

Ionizable lipid-containing lipid nanoparticles (LNPs) have enabled the delivery of RNA for a range of therapeutic applications. In order to optimize safe, targeted, and effective LNP-based RNA delivery platforms, an understanding of the role of composition and pH in their structural properties and self-assembly is crucial, yet there have been few computational studies of such phenomena. Here we present a coarse-grained model of ionizable lipid and mRNA-containing LNPs. Our model allows access to the large length- and time-scales necessary for LNP self-assembly and is mapped and parametrized with reference to all-atom structures and simulations of the corresponding components at compositions typical of LNPs used for mRNA delivery. Our simulations reveal insights into the dynamics of self-assembly of such mRNA-encapsulating LNPs, as well as the subsequent pH change-driven LNP morphology and release of mRNA.

Epik: pKa and Protonation State Prediction through Machine Learning.

Epik version 7 is a software program that uses machine learning for predicting the pKa values and protonation state distribution of complex, druglike molecules. Using an ensemble of atomic graph convolutional neural networks (GCNNs) trained on over 42,000 pKa values across broad chemical space from both experimental and computed origins, the model predicts pKa values with 0.42 and 0.72 pKa unit median absolute and root mean square errors, respectively, across seven test sets. Epik version 7 also generates protonation states and recovers 95% of the most populated protonation states compared to previous versions. Requiring on average only 47 ms per ligand, Epik version 7 is rapid and accurate enough to evaluate protonation states for crucial molecules and prepare ultra-large libraries of compounds to explore vast regions of chemical space. The simplicity and time required for the training allow for the generation of highly accurate models customized to a program's specific chemistry.

Shearing friction behaviour of synthetic polymers compared to a functionalized polysaccharide on biomimetic surfaces: models for the prediction of performance of eco-designed formulations.

The substitution of natural, bio-based and/or biodegradable polymers for those of petrochemical origin in consumer formulations has become an active area of research and development as the sourcing and destiny of material components becomes a more critical factor in product design. These polymers often differ from their petroleum-based counterparts in topology, raw material composition and solution behaviour. Effective and efficient reformulation that maintains comparable cosmetic performance to existing products requires a deep understanding of the differences in frictional behaviour between polymers as a function of their molecular structure. In this work, we simulate the tribological behaviour of three topologically distinct polymers in solution with surfactants and in contact with hair-biomimetic patterned surfaces. We compare a generic functionalized polysaccharide to two performant polymers used in shampoo formulations: a strongly positively charged polyelectrolyte and a zwitterionic copolymer. Topological differences are expected to affect rheological properties, as well as their direct interaction with structured biological substrates. Using a refined Martini-style coarse-grained model we describe the polymer-dependent differences in aggregation behaviour as well as selective interactions with a biomimetic model hair surface. Additionally, we introduce a formalism to characterize the response of the solution to shear as an initial study on lubrication properties, which define the sensorial performance of these systems in cosmetics (i.e., manageability, touch, etc.). The tools and techniques presented in this work illustrate the strength of molecular simulation in eco-design of formulation as a complement to experiment. These efforts help advance our understanding of how we can relate complex atomic-scale solution behaviour to relevant macroscopic properties. We expect these techniques to play an increasingly important role in advancing strategies for green polymer formulation design by providing an understanding for how new polymers could reach and even exceed the level of performance of existing polymers.

Molecular-Level Examination of Amorphous Solid Dispersion Dissolution.

Amorphous solid dispersions (ASDs) are commonly used to orally deliver small-molecule drugs that are poorly water-soluble. ASDs consist of drug molecules in the amorphous form which are dispersed in a hydrophilic polymer matrix. Producing a high-performance ASD is critical for effective drug delivery and depends on many factors such as solubility of the drug in the matrix and the rate of drug release in aqueous medium (dissolution), which is linked to bioperformance. Often, researchers perform a large number of design iterations to achieve this objective. A detailed molecular-level understanding of the mechanisms behind ASD dissolution behavior would aid in the screening, designing, and optimization of ASD formulations and would minimize the need for testing a wide variety of prototype formulations. Molecular dynamics and related types of simulations, which model the collective behavior of molecules in condensed phase systems, can provide unique insights into these mechanisms. To study the effectiveness of these simulation techniques in ASD formulation dissolution, we carried out dissipative particle dynamics simulations, which are particularly an efficient form of molecular dynamics calculations. We studied two stages of the dissolution process: the early-stage of the dissolution process, which focuses on the dissolution at the ASD/water interface, and the late-stage of the dissolution process, where significant drug release would have occurred and there would be a mixture of drug and polymer molecules in a predominantly aqueous environment. Experimentally, we used Fourier transform infrared spectroscopy to study the interactions between drugs, polymers, and water in the dry and wet states and the chromatographic technique to study the rate of drug and polymer release. Both experiments and simulations provided evidence of polymer microstructures and drug-polymer interactions as important factors for the dissolution behavior of the investigated ASDs, consistent with previous work by Pudlas et al. (Eur. J. Pharm. Sci.2015, 67, 21-31). As experimental and simulation results are consistent and complementary, it is clear that there is significant potential for combined experimental and computational research for a detailed understanding of ASD formulations and, hence, formulation optimization.

A New Mixed All-Atom/Coarse-Grained Model: Application to Melittin Aggregation in Aqueous Solution.

We introduce a new mixed resolution, all-atom/coarse-grained approach (AACG), for modeling peptides in aqueous solution and apply it to characterizing the aggregation of melittin. All of the atoms in peptidic components are represented, while a single site is used for each water molecule. With the full flexibility of the peptide retained, our AACG method achieves speedups by a factor of 3-4 for CPU time reduction and another factor of roughly 7 for diffusion. An Ewald treatment permits the inclusion of long-range electrostatic interactions. These characteristics fit well with the requirements for studying peptide association and aggregation, where the system sizes and time scales require considerable computational resources with all-atom models. In particular, AACG is well suited for biologics since changes in peptide shape and long-range electrostatics may play an important role. The application of AACG to melittin, a 26-residue peptide with a well-known propensity to aggregate in solution, serves as an initial demonstration of this technology for studying peptide aggregation. We observed the formation of melittin aggregates during our simulations and characterized the time-evolution of aggregate size distribution, buried surface areas, and residue contacts. Key interactions including π-cation and π-stacking involving TRP19 were also examined. Our AACG simulations demonstrated a clear salt effect and a moderate temperature effect on aggregation and support the molten globule model of melittin aggregates. As a showcase, this work illustrates the useful role for AACG in investigations of peptide aggregation and its potential to guide formulation and design of biologics.

Capturing the multiscale dynamics of membrane protein complexes with all-atom, mixed-resolution, and coarse-grained models.

The structures and dynamics of protein complexes are often challenging to model in heterogeneous environments such as biological membranes. Herein, we meet this fundamental challenge at attainable cost with all-atom, mixed-resolution, and coarse-grained models of vital membrane proteins. We systematically simulated five complex models formed by two distinct G protein-coupled receptors (GPCRs) in the lipid-bilayer membrane on the ns-to-µs timescales. These models, which suggest the swinging motion of an intracellular loop, for the first time, provide the molecular details for the regulatory role of such a loop. For the models at different resolutions, we observed consistent structural stability but various levels of speed-ups in protein dynamics. The mixed-resolution and coarse-grained models show two and four times faster protein diffusion than the all-atom models, in addition to a 4- and 400-fold speed-up in the simulation performance. Furthermore, by elucidating the strengths and challenges of combining all-atom models with reduced resolution models, this study can serve as a guide to simulating other complex systems in heterogeneous environments efficiently.

Melittin Aggregation in Aqueous Solutions: Insight from Molecular Dynamics Simulations.

Melittin is a natural peptide that aggregates in aqueous solutions with paradigmatic monomer-to-tetramer and coil-to-helix transitions. Since little is known about the molecular mechanisms of melittin aggregation in solution, we simulated its self-aggregation process under various conditions. After confirming the stability of a melittin tetramer in solution, we observed­for the first time in atomistic detail­that four separated melittin monomers aggregate into a tetramer. Our simulated dependence of melittin aggregation on peptide concentration, temperature, and ionic strength is in good agreement with prior experiments. We propose that melittin mainly self-aggregates via a mechanism involving the sequential addition of monomers, which is supported by both qualitative and quantitative evidence obtained from unbiased and metadynamics simulations. Moreover, by combining computer simulations and a theory of the electrical double layer, we provide evidence to suggest why melittin aggregation in solution likely stops at the tetramer, rather than forming higher-order oligomers. Overall, our study not only explains prior experimental results at the molecular level but also provides quantitative mechanistic information that may guide the engineering of melittin for higher efficacy and safety.

Macrocycle conformational sampling with MacroModel.

Sampling low energy conformations of macrocycles is challenging due to the large size of many of these molecules and the constraints imposed by the macrocycle. We present a new conformational search method (implemented in MacroModel) that uses brief MD simulations followed by minimization and normal-mode search steps. The method was parametrized using a set of 100 macrocycles from the PDB and CSD. It was then tested on a publicly available data set for which there are published results using alternative methods; we found that when the same force field is used (in this case MMFFs in vacuum), our method tended to identify conformations with lower energies than what the other methods identified. The performance on a new set of 50 macrocycles from the PDB and CSD was also quite good; the mean and median RMSD values for just the ring atoms were 0.60 and 0.33 Å, respectively. However, the RMSD values for macrocycles with more than 30 ring-atoms were quite a bit larger compared to the smaller macrocycles. Possible origins for this and ideas for improving the performance on very large macrocycles are discussed.

Ligand-dependent activation and deactivation of the human adenosine A(2A) receptor.

G-protein-coupled receptors (GPCRs) are membrane proteins with critical functions in cellular signal transduction, representing a primary class of drug targets. Acting by direct binding, many drugs modulate GPCR activity and influence the signaling pathways associated with numerous diseases. However, complete details of ligand-dependent GPCR activation/deactivation are difficult to obtain from experiments. Therefore, it remains unclear how ligands modulate a GPCR's activity. To elucidate the ligand-dependent activation/deactivation mechanism of the human adenosine A2A receptor (AA2AR), a member of the class A GPCRs, we performed large-scale unbiased molecular dynamics and metadynamics simulations of the receptor embedded in a membrane. At the atomic level, we have observed distinct structural states that resemble the active and inactive states. In particular, we noted key structural elements changing in a highly concerted fashion during the conformational transitions, including six conformational states of a tryptophan (Trp246(6.48)). Our findings agree with a previously proposed view that, during activation, this tryptophan residue undergoes a rotameric transition that may be coupled to a series of coherent conformational changes, resulting in the opening of the G-protein binding site. Further, metadynamics simulations provide quantitative evidence for this mechanism, suggesting how ligand binding shifts the equilibrium between the active and inactive states. Our analysis also proposes that a few specific residues are associated with agonism/antagonism, affinity, and selectivity, and suggests that the ligand-binding pocket can be thought of as having three distinct regions, providing dynamic features for structure-based design. Additional simulations with AA2AR bound to a novel ligand are consistent with our proposed mechanism. Generally, our study provides insights into the ligand-dependent AA2AR activation/deactivation in addition to what has been found in crystal structures. These results should aid in the discovery of more effective and selective GPCR ligands.

Modeling Local Structural Rearrangements Using FEP/REST: Application to Relative Binding Affinity Predictions of CDK2 Inhibitors.

Accurate and reliable calculation of protein-ligand binding affinities remains a hotbed of computer-aided drug design research. Despite the potentially large impact FEP (free energy perturbation) may have in drug design projects, practical applications of FEP in industrial contexts have been limited. In this work, we use a recently developed method, FEP/REST (free energy perturbation/replica exchange with solute tempering), to calculate the relative binding affinities for a set of congeneric ligands binding to the CDK2 receptor. We compare the FEP/REST results with traditional FEP/MD (molecular dynamics) results and MM/GBSA (molecular mechanics/Generalized Born Surface Area model) results and examine why FEP/REST performed notably better than these other methods, as well as why certain ligand mutations lead to large increases of the binding affinity while others do not. We also introduce a mathematical framework for assessing the consistency and reliability of the calculations using cycle closures in FEP mutation paths.

Comments on the article "Evaluation of pK(a) estimation methods on 211 druglike compounds".

The recent article "Evaluation of pK(a) Estimation Methods on 211 Druglike Compounds" ( Manchester, J.; et al. J. Chem Inf. Model. 2010, 50, 565-571 ) reports poor results for the program Epik. Here, we highlight likely sources for the poor performance and describe work done to improve the performance. Running Epik in the mode intended to calculate pK(a) values for sequentially adding/removing protons, as needed to reproduce the experimental conditions, improves the root mean squared error (RMSE) from 3.0 to 2.18 for the 85 public compounds available from the paper. Despite this improvement, there are still other programs in the Manchester paper that outperform Epik. The primary reason is that the public portion of the data set is not diverse and Epik is missing a few key functional groups in this data set that are heavily represented. We show that incorporation of these missing functional groups into the Epik training set improves the RMSE for the public compounds to 1.04. Furthermore, these enhancements help improve the overall performance of Epik on a large druglike test set.

Towards the comprehensive, rapid, and accurate prediction of the favorable tautomeric states of drug-like molecules in aqueous solution.

Generating the appropriate protonation states of drug-like molecules in solution is important for success in both ligand- and structure-based virtual screening. Screening collections of millions of compounds requires a method for determining tautomers and their energies that is sufficiently rapid, accurate, and comprehensive. To maximise enrichment, the lowest energy tautomers must be determined from heterogeneous input, without over-enumerating unfavourable states. While computationally expensive, the density functional theory (DFT) method M06-2X/aug-cc-pVTZ(-f) [PB-SCRF] provides accurate energies for enumerated model tautomeric systems. The empirical Hammett-Taft methodology can very rapidly extrapolate substituent effects from model systems to drug-like molecules via the relationship between pK(T) and pK(a). Combining the two complementary approaches transforms the tautomer problem from a scientific challenge to one of engineering scale-up, and avoids issues that arise due to the very limited number of measured pK(T) values, especially for the complicated heterocycles often favoured by medicinal chemists for their novelty and versatility. Several hundreds of pre-calculated tautomer energies and substituent pK(a) effects are tabulated in databases for use in structural adjustment by the program Epik, which treats tautomers as a subset of the larger problem of the protonation states in aqueous ensembles and their energy penalties. Accuracy and coverage is continually improved and expanded by parameterizing new systems of interest using DFT and experimental data. Recommendations are made for how to best incorporate tautomers in molecular design and virtual screening workflows.

ConfGen: a conformational search method for efficient generation of bioactive conformers.

We describe the methodology, parametrization, and application of a conformational search method, called ConfGen, designed to efficiently generate bioactive conformers. We define efficiency as the ability to generate a bioactive conformation within a small total number of conformations using a reasonable amount of computer time. The method combines physics-based force field calculations with empirically derived heuristics designed to achieve efficient searching and prioritization of the ligand's conformational space. While many parameter settings are supported, four modes spanning a range of speed and quality trades-offs are defined and characterized. The validation set used to test the method is composed of ligands from 667 crystal structures covering a broad array of target and ligand classes. With the fastest mode, ConfGen uses an average of 0.5 s per ligand and generates only 14.3 conformers per ligand, at least one of which lies within 2.0 A root-mean-squared deviation of the crystal structure for 96% of the ligands. The most computationally intensive mode raises this recovery rate to 99%, while taking 8 s per ligand. Combining multiple search modes to "fill-in" holes in the conformation space or energy minimizing using an all-atom force field each lead to improvements in the recovery rates at higher resolutions. Overall, ConfGen is at least as good as competing programs at high resolution and demonstrates higher efficiency at resolutions sufficient for many downstream applications, such as pharmacophore modeling.

A role for a specific cholesterol interaction in stabilizing the Apo configuration of the human A(2A) adenosine receptor.

The function of G-protein-coupled receptors is tightly modulated by the lipid environment. Long-timescale molecular dynamics simulations (totaling approximately 3 mus) of the A(2A) receptor in cholesterol-free bilayers, with and without the antagonist ZM241385 bound, demonstrate the instability of helix II in the apo receptor in cholesterol-poor membrane regions. We directly observe that the effect of cholesterol binding is to stabilize helix II against a buckling-type deformation, perhaps rationalizing the observation that the A(2A) receptor couples to G protein only in the presence of cholesterol (Zezula and Freissmuth, 2008). The results suggest a mechanism by which the A(2A) receptor may function as a coincidence detector, activating only in the presence of both cholesterol and agonist. We also observed a previously hypothesized conformation of the tryptophan "rotameric switch" on helix VI in which a phenylalanine on helix V positions the tryptophan out of the ligand binding pocket.

Epik: a software program for pK( a ) prediction and protonation state generation for drug-like molecules.

Epik is a computer program for predicting pK(a) values for drug-like molecules. Epik can use this capability in combination with technology for tautomerization to adjust the protonation state of small drug-like molecules to automatically generate one or more of the most probable forms for use in further molecular modeling studies. Many medicinal chemicals can exchange protons with their environment, resulting in various ionization and tautomeric states, collectively known as protonation states. The protonation state of a drug can affect its solubility and membrane permeability. In modeling, the protonation state of a ligand will also affect which conformations are predicted for the molecule, as well as predictions for binding modes and ligand affinities based upon protein-ligand interactions. Despite the importance of the protonation state, many databases of candidate molecules used in drug development do not store reliable information on the most probable protonation states. Epik is sufficiently rapid and accurate to process large databases of drug-like molecules to provide this information. Several new technologies are employed. Extensions to the well-established Hammett and Taft approaches are used for pK(a) prediction, namely, mesomer standardization, charge cancellation, and charge spreading to make the predicted results reflect the nature of the molecule itself rather just for the particular Lewis structure used on input. In addition, a new iterative technology for generating, ranking and culling the generated protonation states is employed.

Transmembrane peptide-induced lipid sorting and mechanism of Lalpha-to-inverted phase transition using coarse-grain molecular dynamics.

Molecular dynamics results are presented for a coarse-grain model of 1,2-di-n-alkanoyl-sn-glycero-3-phosphocholine, water, and a capped cylindrical model of a transmembrane peptide. We first demonstrate that different alkanoyl-length lipids are miscible in the liquid-disordered lamellar (Lalpha) phase. The transmembrane peptide is constructed of hydrophobic sites with hydrophilic caps. The hydrophobic length of the peptide is smaller than the hydrophobic thickness of a bilayer consisting of an equal mixture of long and short alkanoyl tail lipids. When incorporated into the membrane, a meniscus forms in the vicinity of the peptide and the surrounding area is enriched in the short lipid. The meniscus region draws water into it. In the regions that are depleted of water, the bilayers can fuse. The lipid headgroups then rearrange to solvate the newly formed water pores, resulting in an inverted phase. This mechanism appears to be a viable pathway for the experimentally observed Lalpha-to-inverse hexagonal (HII) peptide-induced phase transition.

Preparation and physical characterization of pure beta-carotene.

Pure all-trans beta-carotene has been prepared on the 10's of grams scale by isothermal Fractional Dissolution (FD) of commercial laboratory samples in tetrahydrofuran (THF). beta-Carotene purified in this way is black, with a faint brownish tinge. The electronic spectra of black samples extend into the near infrared, with end-absorption past 750 nm. Black samples react directly with dioxygen under mild conditions to yield the familiar orange or red powders. Pure beta-carotene rigorously obeys Beer's Law in octane over the entire UV-Vis spectral range, while commercial laboratory samples and recrystallized samples do not. NMR self-diffusion coefficient data demonstrate that beta-carotene exists as simple molecular solutions in octane and toluene. The anomalously high crystallinity of beta-carotene can be attributed (from analysis using molecular mechanics) to the facts that: (1) the number of theoretically possible conformers of beta-carotene is extremely small, and (2) only a small fraction of these (ca. 12%, or 127) may actually exist in fluid phases.