Conventional Molecular Dynamics Research Articles

Dissecting the Effect of Temperature on Hyperthermophilic Pf2001 Esterase Dimerization by Molecular Dynamics.

Pf2001 esterase (Pf2001) from Pyrococcus furiosus has hyperthermophilic properties and exerts a biocatalytic function in a dimeric state. Crystal structures revealed that the structural rearrangement of the cap domain is responsible for the Pf2001 dimer formation. However, the details of the cap domain remodeling and the effects of temperature on the dimerization process remain elusive at the molecular level, taking into account that experimental methods are difficult to capture the dynamic process of dimerization to some extent. Herein, four dimer models based on the monomeric crystal structure (PDB ID: 5G59) were constructed to investigate the conformational transition details and temperature effects in the dimerization by conventional molecular dynamics and accelerated molecular dynamics simulations. Our simulation results indicate that the monomer undergoes a conformational change into a "preparatory state" at high temperatures, which is more favorable for its transformation into a stable dimer. The subsequent free energy landscape analysis further identifies four intermediate states (from separated state to dimeric state) and discloses that a more accessible α-helix driven by stronger hydrophobic interactions induces a rearrangement of the cap domain, displaying a "tic-tac-toe" activation feature that is important for stabilizing the dimer interface and facilitating the formation of hydrophobic pockets. In addition, the electrostatic potential surface analysis illustrates that the weaker electrostatic repulsion (Lys and Arg) in the dimer interface at high temperatures is also a key factor for dimer stabilization. Altogether, our results can provide molecular-level insight into the dimer formation process of hyperthermophilic esterase and would be useful to understand the enzymatic specificity of α/β-hydrolase.

Dihydroartemisinin binds human PI3K-β affinity pocket and forces flat conformation in P-loop MET783: A molecular dynamics study

Artemisinin and its semi-synthetic derivatives are not only indicated for malaria but also cancer, inflammatory and autoimmune diseases. Its inflammatory and immunosuppressive target is PI3K/AKT pathways. The structural and kinetic aspect of the PI3K inhibition was investigated in the current study using computational approaches. Binding energies of dihydroartemisinin (DHA) to p110-PI3K-β was computed using the MMPBSA method in comparison with the standard inhibitor (GD9). Kinetic parameter (Kon/Koff) was also evaluated for the complexes using adaptive sampling protocols and Markov state model analysis. p110-PI3K- β dynamics and community network analysis were also performed following conventional Molecular dynamics simulation. The results showed −63.99 ± 1.53 and −74.14 ± 3.47 (Kj/mol) binding energies for DHA and GD9 respectively. Kon/Koff estimates for DHA and GD9 are 12.4, and 2.13 (M−1) respectively. Analysis of the trajectories showed that DHA selectively partitions into p110-PI3K- β affinity pocket, forces open conformation, and kept catalytic pocket-M783 in a flat conformation whilst forcing large displacement around the C2-domain. In conclusion, DHA is a high affinity (slow-binding, slow-dissociating), flat-conformation p110-PI3K- β inhibitor.

Comprehensive Approach to Simulating Large Scale Conformational Changes in Biological Systems Utilizing a Path Collective Variable and New Barrier Restraint.

Conformational sampling of complex biomolecules is an emerging frontier in drug discovery. Advances in lab-based structural biology and related computational approaches like AlphaFold have made great strides in obtaining static protein structures for biologically relevant targets. However, biology is in constant motion, and many important biological processes rely on conformationally driven events. Conventional molecular dynamics (MD) simulations run on standard hardware are impractical for many drug design projects, where conformationally driven biological events can take microseconds to milliseconds or longer. An alternative approach is to focus the search on a limited region of conformational space defined by a putative reaction coordinate (i.e., path collective variable). The search space is typically limited by applying restraints, which can be guided by insights about the underlying biological process of interest. The challenge is striking a balance between the degree to which the system is constrained and still allowing for natural motions along the path. A plethora of restraints exist to limit the size of conformational search space, although each has drawbacks when simulating complex biological motions. In this work, we present a three-stage procedure to construct realistic path collective variables (PCVs) and introduce a new kind of barrier restraint that is particularly well suited for complex conformationally driven biological events, such as allosteric modulations and conformational signaling. The PCV presented here is all-atom (as opposed to C-alpha or backbone only) and is derived from all-atom MD trajectory frames. The new restraint relies on a barrier function (specifically, the scaled reciprocal function), which we show is particularly beneficial in the context of molecular dynamics, where near-hard-wall restraints are needed with zero tolerance to restraint violation. We have implemented our PCV and barrier restraint within a hybrid sampling framework that combines well-tempered metadynamics and extended-Lagrangian adaptive biasing force (meta-eABF). We use three particular examples of high pharmaceutical interest to demonstrate the value of this approach: (1) sampling the distance from ubiquitin to a protein of interest within the supramolecular cullin-RING ligase complex, (2) stabilizing the wild-type conformation of the oncogenic mutant JAK2-V617F pseudokinase domain, and (3) inducing an activated state of the stimulator of interferon genes (STING) protein observed upon ligand binding. For examples 2 and 3, we present statistical analysis of meta-eABF free energy estimates and, for each case, code for reproducing this work.

Deep Boosted Molecular Dynamics: Accelerating Molecular Simulations with Gaussian Boost Potentials Generated Using Probabilistic Bayesian Deep Neural Network.

We have developed a new deep boosted molecular dynamics (DBMD) method. Probabilistic Bayesian neural network models were implemented to construct boost potentials that exhibit Gaussian distribution with minimized anharmonicity, thereby allowing for accurate energetic reweighting and enhanced sampling of molecular simulations. DBMD was demonstrated on model systems of alanine dipeptide and the fast-folding protein and RNA structures. For alanine dipeptide, 30 ns DBMD simulations captured up to 83-125 times more backbone dihedral transitions than 1 μs conventional molecular dynamics (cMD) simulations and were able to accurately reproduce the original free energy profiles. Moreover, DBMD sampled multiple folding and unfolding events within 300 ns simulations of the chignolin model protein and identified low-energy conformational states comparable to previous simulation findings. Finally, DBMD captured a general folding pathway of three hairpin RNAs with the GCAA, GAAA, and UUCG tetraloops. Based on a deep learning neural network, DBMD provides a powerful and generally applicable approach to boosting biomolecular simulations. DBMD is available with open source in OpenMM at https://github.com/MiaoLab20/DBMD/.

Using Density Functional Theory for Testing the Robustness of Mobile-Proton Molecular Dynamics Simulations on Electrosprayed Ions: Structural Implications for Gaseous Proteins.

Current experiments only provide low-resolution information on gaseous protein ions generated by electrospray ionization (ESI). Molecular dynamics (MD) simulations can yield complementary insights. Unfortunately, conventional MD does not capture the mobile nature of protons in gaseous proteins. Mobile-proton MD (MPMD) overcomes this limitation. Earlier MPMD data at 300 K indicated that protein ions generated by "native" ESI retain solution-like structures with a hydrophobic core and zwitterionic exterior [Bakhtiari, M.; Konermann, L. J. Phys. Chem. B 2019, 123, 1784-1796]. MPMD redistributes protons using electrostatic and proton affinity calculations. The robustness of this approach has never been scrutinized. Here, we close this gap by benchmarking MPMD against density functional theory (DFT) at the B3LYP/6-31G* level, which is well suited for predicting proton affinities. The computational cost of DFT necessitated the use of small peptides. The MPMD energetic ranking of proton configurations was found to be consistent with DFT single-point energies, implying that MPMD can reliably identify favorable protonation sites. Peptide MPMD runs converged to DFT-optimized structures only when applying 300-500 K temperature cycling, which was necessary to prevent trapping in local minima. Temperature cycling MPMD was then applied to gaseous protein ions. Native ubiquitin converted to slightly expanded structures with a zwitterionic core and a nonpolar exterior. Our data suggest that such inside-out protein structures are intrinsically preferred in the gas phase, and that they form in ESI experiments after moderate collisional excitation. This is in contrast to native ESI (with minimal collisional excitation, simulated by MPMD at 300 K), where kinetic trapping promotes the survival of solution-like structures. In summary, this work validates the MPMD approach for simulations on gaseous peptides and proteins.

Studying the impact of phase behavior in the morphology of molecular dynamics models of bitumen

Conventional Molecular Dynamics (MD) models of bitumen are built by homogeneously mixing molecules in a volume without considering that the molecules in bitumen are known to exhibit phase behavior and form distinctive molecular arrangements. These are known to have a significant impact in the behavior of bitumen, and considering their existence is paramount in producing improved representations of bitumen using computational models. This study explores whether MD models of bitumen that are conventionally assumed to be in equilibrium can still undergo significant phase separation over considerably long simulation times. It also aims to establish a more formal pathway to build and study models with highly heterogeneous arrangements of their molecules. Moreover, it aims to evaluate whether the presence of distinct morphologies have a significant impact in numerous physical properties of bitumen. The study shows that conventional and widely used models of bitumen exhibit significant molecular rearrangements over long times (>360 ns). It also shows that building heterogeneous morphologies is possible and result in energetically favorable conformations. Moreover, it proves that studying properties regularly used to validate MD models of bitumen (e.g., density) are insufficient in assessing the impact of different morphologies; more thorough methods are required to evaluate them.

The α-Synuclein Monomer May Have Different Misfolding Mechanisms in the Induction of α-Synuclein Fibrils with Different Polymorphs

The aggregation of alpha-synuclein (α-Syn) is closely related to the occurrence of some neurodegenerative diseases such as Parkinson's disease. The misfolding of α-Syn monomer plays a key role in the formation of aggregates and extension of fibril. However, the misfolding mechanism of α-Syn remains elusive. Here, three different α-Syn fibrils (isolated from a diseased human brain, generated by in vitro cofactor-tau induction, and obtained by in vitro cofactor-free induction) were selected for the study. The misfolding mechanisms of α-Syn were uncovered by studying the dissociation of the boundary chains based on the conventional molecular dynamics (MD) and Steered MD simulations. The results showed that the dissociation paths of the boundary chains in the three systems were different. According to the reverse process of dissociation, we concluded that in the human brain system, the binding of the monomer and template starts from the C-terminal and gradually misfolds toward the N-terminal. In the cofactor-tau system, the monomer binding starts from residues 58-66 (contain β3), followed by the C-terminal coil (residues 67-79). Then, the N-terminal coil (residues 36-41) and residues 50-57 (contain β2) bind to the template, followed by residues 42-49 (contain β1). In the cofactor-free system, two misfolding paths were found. One is that the monomer binds to the N/C-terminal (β1/β6) and then binds to the remaining residues. The other one is that the monomer binds sequentially from the C- to N-terminal, similar to the human brain system. Furthermore, in the human brain and cofactor-tau systems, electrostatic interactions (especially from residues 58-66) are the main driving force during the misfolding process, whereas in the cofactor-free system, the contributions of electrostatic and van der Waals interactions are comparable. These results may provide a deeper understanding for the misfolding and aggregation mechanism of α-Syn.

Nonequilibrium molecular dynamics for accelerated computation of ion–ion correlated conductivity beyond Nernst–Einstein limitation

Condensed matters with high ionic conductivities are crucial in various solid devices such as solid-state batteries. The conduction is characterized by the cooperative ionic motion associated with the high carrier density. However, the high cost of computing correlated ionic conductivities has forced almost all ab initio molecular dynamics (MD) to rely on the Nernst–Einstein dilute-solution approximation, which ignores the cross-correlation effect. Here we develop a chemical color-diffusion nonequilibrium MD (CCD-NEMD) method, which enables to calculate the correlated conductivities with fewer sampling steps than the conventional MD. This CCD-NEMD is demonstrated to well evaluate the conductivities in the representative solid electrolyte bulk Li10GeP2S12 and Li7La3Zr2O12. We also applied CCD-NEMD to the grain boundary of Li7La3Zr2O12 and demonstrated its applicability for calculating interfacial local conductivities, which is essential for investigating grain boundaries and composite electrolytes. CCD-NEMD can provide further accurate understanding of ionics with ionic correlations and promote developing solid devices.

Mechanistic Insights into the Binding of Different Positron Emission Tomography Tracers to Chronic Traumatic Encephalopathy Tau Protofibrils.

Chronic traumatic encephalopathy (CTE) is a neurodegenerative disease associated with exposure to repetitive head impacts, and it is neuropathologically defined as the accumulation of abnormally hyperphosphorylated tau (p-tau). Early detection of p-tau in the brain is of great value in the prevention and treatment of CTE. Previous experimental studies reported that positron emission tomography (PET) technique using several tau tracers are available for imaging certain neurodegenerative diseases. However, few studies have focused on the development of CTE tau tracers. In this work, we performed conventional molecular docking and molecular dynamics simulations to address the binding properties and mechanisms of PET tracers (18F-PM-PBB3, 18F-CBD-2115, 18F-PI-2620, 18F-RO-948, 18F-MK-6240, and 18F-flortaucipir) to CTE tau protofibrils. The results show that the hydrophobic cavity and the top of the concave structure of CTE tau protofibrils are the preferred binding sites for the six tracers, and 18F-PM-PBB3 has the most competitive binding affinity to CTE tau protofibrils. Further investigation into the binding patterns of the six tracers to the CTE tau protofibrils showed that 18F-CBD-2115 and 18F-PM-PBB3 have a high number of H-bonds and hydrophobic contacts with tau protofibrils, resulting in strong hydrogen bonding and hydrophobic interactions; 18F-flortaucipir/18F-PI-2620 and 18F-PI-2620/18F-RO-948 form more intense π-π and cation-π interactions with tau protofibrils, respectively. Subsequently, we conducted a detailed analysis of the binding mechanism of 18F-PM-PBB3 to CTE tau protofibrils. The benzothiazole ring of 18F-PM-PBB3 exhibits stronger π-π stacking and cation-π interactions with tau protofibrils than the pyridine ring and forms a more concentrated T-shaped π-π stacking pattern. This study contributes to understanding the binding mechanism of PET tracers to CTE tau protofibrils and provides new insights into the design of potential novel tracers.

Protex-A Python utility for proton exchange in molecular dynamics simulations.

Protex is an open-source program that enables proton exchanges of solvent molecules during molecular dynamics simulations. While conventional molecular dynamics simulations do not allow for bond breaking or formation, protex offers an easy-to-use interface to augment these simulations and define multiple proton sites for (de-)protonation using a single topology approach with two different λ-states. Protex was successfully applied to a protic ionic liquid system, where each molecule is prone to (de-)protonation. Transport properties were calculated and compared to experimental values and simulations without proton exchange.

Molecular Modeling of Self-Assembling Peptides.

Peptide epitopes mediate as many as 40% of protein-protein interactions and fulfill signaling, inhibition, and activation roles within the cell. Beyond protein recognition, some peptides can self- or coassemble into stable hydrogels, making them a readily available source of biomaterials. While these 3D assemblies are routinely characterized at the fiber level, there are missing atomistic details about the assembly scaffold. Such atomistic detail can be useful in the rational design of more stable scaffold structures and with improved accessibility to functional motifs. Computational approaches can in principle reduce the experimental cost of such an endeavor by predicting the assembly scaffold and identifying novel sequences that adopt said structure. Yet, inaccuracies in physical models and inefficient sampling have limited atomistic studies to short (two or three amino acid) peptides. Given recent developments in machine learning and advances in sampling strategies, we revisit the suitability of physical models for this task. We use the MELD (Modeling Employing Limited Data) approach to drive self-assembly in combination with generic data in cases where conventional MD is unsuccessful. Finally, despite recent developments in machine learning algorithms for protein structure and sequence predictions, we find the algorithms are not yet suited for studying the assembly of short peptides.

Experimental study of a standalone membrane water desalination unit fully powered by solar energy

In this study, a fully solar-driven water desalination system is designed, fabricated, and tested. The proposed system is designed for remote, isolated areas with limited electricity and freshwater access. For commercial PVDF membranes, the performance of the MD unit is evaluated at various feed temperatures, flow rates, and salt concentrations. In addition, the effect of MD cell design was examined to determine the effect of the obstacles on the feed channel. Increasing the temperature by 20 degrees from 55 to 75 °C enhances productivity by 165 % from 6.2 to 16.4 kg/m2.h. At different operating conditions, the new design of the MD cell with cylindrical obstacles attained higher flux productivity compared to the conventional MD cell design. The maximum attained freshwater output flux was around 15 kg/m2.h at a feed inlet temperature of 65 °C and feed mass flowrate of 66 kg/h for the MD cell with cylindrical obstacles. In the presence of obstacles, a reduction in the specific energy consumption (SEC) by 35–15 % is associated with an increase in the gain output ratio (GOR) by 50–20 % compared to the case of no obstacles. Observations indicate that a small PV module with an area of 0.377 m2 is sufficient to operate the membrane unit for contentious 15 h. Collectively, using solar MD units collectively has the potential to be used for decentralized desalination of water.

A steered molecular dynamics investigation on the role of glycation on the energy dissipation capacity of osteopontin at mineral interfaces in bone.

Type-2 Diabetes (T2D) is associated with a 3-fold increase in bone fracture risk. T2D impairs bone quality, whereby the intrinsic material properties of the bone matrix are altered, through glycation of collagenous and non-collagenous proteins (NCPs). There is a distinct lack of understanding of how altered protein configurations impair whole-bone biomechanics, particularly in the case of NCPs since their structural roles in bone are only now emerging. In this study, the effect of glycation on the energy dissipation behaviour of OPN, the most abundant NCP, on HydroxyApatite (HAp) mineral surfaces was investigated. The adsorption modes of OPN on the surface was obtained by conventional MD by comparing the residues residence time for each protein orientation with their average values. Then, energy dissipation of OPN was obtained through Steered Molecular Dynamics (SMD) by pulling simulations perpendicular to the surface. Two different glycation possibilities were considered whereby (a) Glycation occurred prior to Adsorption (GA) or (b) Adsorption took place before Glycation (AG). The simulation results showed that AG has the same adsorption modes as the healthy OPN, while the GA experienced a different mode. It was found that energy dissipation in absence of glycation was EOPN = 2646 kcal/mol. For glycated groups energy dissipation values were EOPN-AG = 2324 kcal/mol and EOPN-GA = 773 kcal/mol. The simulation results showed that OPN had reduced capacity for energy dissipation for either glycation case. However, it was notable that if OPN glycation occurred before its adsorption on HAp, its energy dissipation capacity was greatly reduced, which was caused by changes to the adsorption mode itself. The current study for the first time showed that the glycation sequence could negatively affects energy dissipation behaviour of OPN, which could have consequence for bone toughness.

Optimization of hydrogen-deuterium exchange ensemble reweighting (HDXer) to aid the biophysical characterization of proteins.

Hydrogen-deuterium exchange mass spectrometry (HDX-MS) experiments provide unique insights into native protein structure and dynamics in solution but have limited spatial resolution. Molecular dynamics (MD) simulations, in contrast, allow for the exploration of protein conformational landscapes with an atomistic level of resolution. However, forcefield inaccuracies, crystallographic biases, and timescale limitations may compromise the validity of the final structural ensemble. Hydrogen Deuterium eXchange Ensemble Reweighting (HDXer) is a maximum entropy approach that was recently developed to (1) quantitatively predict hydrogen-deuterium exchange from a modeled conformational ensemble and (2) refine the candidate ensemble to best fit an experimental HDX-MS dataset. Here, we compare how various computational and experimental parameters impact the performance of HDXer, using extracellular regulated kinase 2 (ERK2) and hen egg white lysozyme (HEWL) as model proteins. Conventional molecular dynamics simulations of ERK2 and HEWL were performed in triplicate using the CHARMM36m, AMBER FF19SB, and Drude forcefields. After predicting the HDX rates of each candidate ensemble, HDXer was used to fit the simulated structures to the corresponding HDX-MS datasets. To assess the influence of simulation length, truncated and strided copies of the full-length simulations were generated. Other factors, including peptide charge state selection, ion mobility mass spectrometry, and independent HDX-MS data acquisitions, were also tested. We find that while some variables (e.g., high vs. low peptide charge states) do not have a significant effect, other parameters (e.g., early vs. late timepoint selection) directly correlate with the disagreement between calculated and experimental deuterium uptake. The results from this benchmarking inform the development of a recommended HDXer workflow and best practices for optimal data acquisition.

Determining the confromational ensemble of the HIV-1 LTR G-quadruplexes through gaussian-accelerated molecular dynamics using the drude polarizable force field.

Human immunodeficiency virus (HIV) impacts nearly 37.6 million people globally, but there is no effective cure for the infections it causes. HIV-1 is a retrovirus containing a single-stranded RNA genome. After reverse transcription, the 5’ long terminal repeat (LTR) acts as a promoter for viral replication. Two DNA G-quadruplexes (GQs), known as LTR-III and LTR-IV, form within this promoter and inversely regulate viral gene expression. The LTR-III GQ downregulates viral expression and is a quadruplex: duplex hybrid whereas LTR-IV upregulates viral expression and is a parallel GQ containing a single-nucleotide bulge. Their distinct topologies and functional roles suggest their potential to act as selective drug targets for HIV-1.To gain a greater understanding of structure and dynamics of the LTR-III and LTR-IV GQs, we performed conventional molecular dynamics simulations using the Drude polarizable force field. The loop nucleotides in both structures were highly flexible, sampling a variety of conformational states. The LTR-IV loop also formed base pair interactions that limited ion access to its loop cavity and further highlighted two distinct ion binding sites near the backbone tetrad guanines. Within the LTR-III duplex loop, an additional Watson-Crick base pair formed between Ade4:Thy14 and persisted for 40% of the snapshots collected in the total 4-μs simulation time. Analysis of the electric field exerted on Thy14 indicates a preference for base pairing with Ade4 compared to being unpaired as in the NMR ensemble. Gaussian-accelerated molecular dynamics simulations were also performed to identify high and low energy conformational states for subsequent drug targeting. Our results emphasize the importance of electronic polarization in GQ dynamics and can guide future development of antiviral therapeutics.

Development of CHARMM additive potential energy parameters for α-methyl amino acids.

Potential energy parameters for α-methyl amino acids were generated with ab initio calculations on α-methyl-N-acetylalanyl-N′-methylamide (the α-methyl “alanine dipeptide”) which served as an input to a grid-based correction to the backbone torsional potential (known as CMAP) consistent with the CHARMM36m additive protein force field. The new parameters were validated by comparison with experimentally determined helicities of the 22 residue C-terminal peptide (H10) from apolipoprotein A1 and five α-methylated variants in water and 0.3:0.7 trifluoroethanol (TFE)/water. Conventional molecular dynamics simulation totaling 30 μs for each peptide is in overall good agreement with the experiment, including the increased helicity in 30% TFE. An additional 500 ns of simulation using two-dimensional dihedral biasing (bpCMAP) replica exchange reduced left-handed conformations, increased right-handed helices, and thereby mostly decreased agreement with the experiment. Analysis of side chain–side chain salt bridges suggests that the overestimation of the helical content may be, in part, due to such interactions. The increased helicity of the peptides in 30% TFE arises from decreased hydrogen bonding of the backbone atoms to water and a concomitant increase in intramolecular backbone hydrogen bonds. The present inclusion of parameters for amino acids methylated at the α-carbon can be extended to all amino acids that contain two β-carbons such as isovaline and stapled peptides.

Exploring the conformational dynamics of the disordered Lid region of HDM2 N-terminal domain with molecular dynamics simulation.

HDM2 is the negative regulator of the tumor suppressor p53, which plays an important role in controlling intracellular homeostasis. The p53 binding pocket of HDM2 is located at the N-terminal domain consisting of the two regions, Lid and Core. The former is a disordered region located at the N-terminus of the domain, while the latter includes the p53 binding site. Lid dynamics modulate the relative orientation of Lid and Core and thereby expose the p53 binding site. In addition, the dynamics of Lid would be affected by macromolecular crowding environments in the cell. We therefore aim to investigate the conformational dynamics of Lid in HDM2 using molecular dynamics (MD) simulations in dilute solution as well as protein-crowded solution. As the first step of this work, we performed conventional MD simulations of HDM2 for 1 μs in a dilute solution using the three force fields, namely, CHARMM36m, AMBER ff19SB, and AMBER ff99SBws-STQ. Since Lid is a highly disordered region, selecting the appropriate force field is essential for comparing the simulation results with experiments and understanding the nature of Lid conformational dynamics. All the MD simulations were performed using GENESIS software on Supercomputer HOKUSAI BigWaterfall at RIKEN in Japan. Two structures derived from the NMR structure (PDB ID: 1Z1M) were used as initial structures in the simulations. Based on the difference in the distance between Lid and Core, the two structures are named “bound” and “unbound” states, respectively. The MD simulation results showed that the secondary structure in Lid and the distance between Lid and Core are strongly dependent on the force fields used. We discuss which force field is the most reliable for studying conformational dynamics of Lid in a dilute and protein-crowded solution.

Theoretical Studies on Selectivity of HPK1/JAK1 Inhibitors by Molecular Dynamics Simulations and Free Energy Calculations.

Hematopoietic progenitor kinase 1 (HPK1) is a negative regulator of T cell receptor, which has been regarded as a potential target for immunotherapy. Yu et al. observed the off-target effect of the high-throughput screening HPK1 kinase inhibitor hits on JAK1 kinase. The off-target effect is usually due to the lack of specificity of the drug, resulting in toxic side effects. Therefore, exploring the mechanisms to selectively inhibit HPK1 is critical for developing effective and safe inhibitors. In this study, two indazole compounds as HPK1 inhibitors with different selectivity towards JAK1 were used to investigate the selectivity mechanism using multiple computational methods, including conventional molecular dynamics simulations, binding free energy calculations and umbrella sampling simulations. The results indicate that the salt bridge between the inhibitor and residue Asp101 of HPK1 favors their selectivity towards HPK1 over JAK1. Information obtained from this study can be used to discover and design more potent and selective HPK1 inhibitors for immunotherapy.

Predicting surfactant phase behavior with a molecularly informed field theory

HypothesisThe computational study of surfactants and self-assembly is challenging because 1) models need to reflect chemistry-specific interactions, and 2) self-assembled structures are difficult to equilibrate with conventional molecular dynamics. We propose to overcome these challenges with a multiscale simulation approach where relative entropy minimization transfers chemically-detailed information from all-atom (AA) simulations to coarse-grained (CG) models that can be simulated using field-theoretic methods. Field-theoretic simulations are not limited by intrinsic physical time scales like diffusion and allow for rigorous equilibration via free energy minimization. This approach should enable the study of properties that are difficult to obtain by particle-based simulations. Simulation WorkWe apply this workflow to sodium dodecylsulfate. To ensure chemical fidelity we present an AA force field calibrated against interfacial tension experiments. We generate CG models from AA simulation trajectories and show that particle-based and field-theoretic simulations of the CG model reproduce AA simulations and experimental measurements. FindingsThe workflow captures the complex balance of interactions in a multicomponent system ultimately described by an atomistic model. The resulting CG models can study complex 3D phases like double or alternating gyroids, and reproduce salt effects on properties like aggregation number and shape transitions.

Deciphering the Shared and Specific Drug Resistance Mechanisms of Anaplastic Lymphoma Kinase via Binding Free Energy Computation.

Anaplastic lymphoma kinase (ALK), a tyrosine receptor kinase, has been proven to be associated with the occurrence of numerous malignancies. Although there have been already at least 3 generations of ALK inhibitors approved by FDA or in clinical trials, the occurrence of various mutations seriously attenuates the effectiveness of the drugs. Unfortunately, most of the drug resistance mechanisms still remain obscure. Therefore, it is necessary to reveal the bottom reasons of the drug resistance mechanisms caused by the mutations. In this work, on the basis of verifying the accuracy of 2 main kinds of binding free energy calculation methodologies [end-point method of Molecular Mechanics with Poisson-Boltzmann/Generalized Born and Surface Area (MM/PB(GB)SA) and alchemical method of Thermodynamic Integration (TI)], we performed a systematic analysis on the ALK systems to explore the underlying shared and specific drug resistance mechanisms, covering the one-drug-multiple-mutation and multiple-drug-one-mutation cases. Through conventional molecular dynamics (cMD) simulation in conjunction with MM/PB(GB)SA and umbrella sampling (US) in conjunction with contact network analysis (CNA), the resistance mechanisms of the in-pocket, out-pocket, and multiple-site mutations were revealed. Especially for the out-pocket mutation, a possible transfer chain of the mutation effect was revealed, and the reason why different drugs exhibited various sensitivities to the same mutation was also uncovered. The proposed mechanisms may be prevalent in various drug resistance cases.