Folding simulations for proteins with diverse topologies are accessible in days with a physics-based force field and implicit solvent

Rational design of glycosaminoglycan binding cyclic peptides using cPEPmatch

Cyclic peptides present a robust platform for drug design, offering high specificity and stability due to their conformationally constrained structures. In this study, we introduce an updated version of the Cyclic Peptide Matching program (cPEPmatch) tailored for the identification of cyclic peptides capable of mimicking protein-glycosaminoglycan (GAG) binding sites. We focused on engineering cyclic peptides to replicate the GAG-binding affinity of antithrombin III (ATIII), a protein that plays a crucial role in modulating anticoagulation through interaction with the GAG heparin. By integrating computational and experimental methods, we successfully identified a cyclic peptide binder with promising potential for future optimization. MD simulations and MM-GBSA calculations were used to assess binding efficacy, supplemented by umbrella sampling to approximate free energy landscapes. The binding specificity was further validated through NMR and ITC experiments. Our findings demonstrate that the computationally designed cyclic peptides effectively target GAGs, suggesting their potential as novel therapeutic agents. This study advances our understanding of peptide-GAG interactions and lays the groundwork for future development of cyclic peptide-based therapeutics.

Implicit Solvent with Explicit Ions Generalized Born Model in Molecular Dynamics: Application to DNA

The ion atmosphere surrounding highly charged biomolecules, such as nucleic acids, is crucial for their dynamics, structure, and interactions. Here, we develop an approach for the explicit treatment of ions within an implicit solvent framework suitable for atomistic simulations of biomolecules. The proposed implicit solvent/explicit ions model, GBION, is based on a modified generalized Born (GB) model; it includes separate, modified GB terms for solute-ion and ion-ion interactions. The model is implemented in the AMBER package (version 24), and its performance is thoroughly investigated in atomistic molecular dynamics (MD) simulations of double-stranded DNA on a microsecond time scale. The aggregate characteristics of monovalent (Na+ and K+) and trivalent (Cobalt Hexammine, CoHex3+) counterion distributions around double-stranded DNA predicted by the model are in reasonable agreement with the experiment (where available), all-atom explicit water MD simulations, and the expectation from the Manning condensation theory. The radial distributions of monovalent cations around DNA are reasonably close to the ones obtained using the explicit water model: expressed in units of energy, the maximum deviations of local ion concentrations from the explicit solvent reference are within 1 kBT, comparable to the corresponding deviations expected between different established explicit water models. The proposed GBION model is able to simulate DNA fragments in a large volume of solvent with explicit ions with little additional computational overhead compared with the fully implicit GB treatment of ions. Ions simulated using the developed model explore conformational space at least 2 orders of magnitude faster than in the explicit solvent. These advantages allowed us to observe and explore an unexpected "stacking" mode of DNA condensation in the presence of trivalent counterions (CoHex3+) that was revealed by recent experiments.

Force Field Limitations of All-Atom Continuous Constant pH Molecular Dynamics

All-atom constant pH molecular dynamics simulations offer a powerful tool for understanding pH-mediated and proton-coupled biological processes. As the protonation equilibria of protein sidechains are shifted by electrostatic interactions and desolvation energies, pK a values calculated from the constant pH simulations may be sensitive to the underlying protein force field and water model. Here we investigated the force field dependence of the all-atom particle mesh Ewald (PME) continuous constant pH (PME-CpHMD) simulations of a mini-protein BBL. The replica-exchange titration simulations based on the Amber ff19SB and ff14SB force fields with the respective water models showed significantly overestimated pK a downshifts for a buried histidine (His166) and for two glutamic acids (Glu141 and Glu161) that are involved in salt-bridge interactions. These errors (due to undersolvation of neutral histidines and overstabilization of salt bridges) are consistent with the previously reported pK a's based on the CHARMM c22/CMAP force field, albeit in larger magnitudes. The pK a calculations also demonstrated that ff19SB with OPC water is significantly more accurate than ff14SB with TIP3P water, and the salt-bridge related pK a downshifts can be partially alleviated by the atom-pair specific Lennard-Jones corrections (NBFIX). Together, these data suggest that the accuracies of the protonation equilibria of proteins from constant pH simulations can significantly benefit from improvements of force fields.

Machine-learned molecular mechanics force fields from large-scale quantum chemical data

The development of reliable and extensible molecular mechanics (MM) force fields-fast, empirical models characterizing the potential energy surface of molecular systems-is indispensable for biomolecular simulation and computer-aided drug design. Here, we introduce a generalized and extensible machine-learned MM force field, espaloma-0.3, and an end-to-end differentiable framework using graph neural networks to overcome the limitations of traditional rule-based methods. Trained in a single GPU-day to fit a large and diverse quantum chemical dataset of over 1.1 M energy and force calculations, espaloma-0.3 reproduces quantum chemical energetic properties of chemical domains highly relevant to drug discovery, including small molecules, peptides, and nucleic acids. Moreover, this force field maintains the quantum chemical energy-minimized geometries of small molecules and preserves the condensed phase properties of peptides and folded proteins, self-consistently parametrizing proteins and ligands to produce stable simulations leading to highly accurate predictions of binding free energies. This methodology demonstrates significant promise as a path forward for systematically building more accurate force fields that are easily extensible to new chemical domains of interest.

Dissipative Particle Dynamics Simulation of Protein Folding in Explicit and Implicit Solvents: Coarse-Grained Model for Atomic Resolution

Advancements have been made to dissipative particle dynamics (DPD), a robust coarse-grained (CG) simulation method, to study the folded structures of four miniproteins (1L2Y, 1WN8, 1YRF, and 2I9M) in explicit and implicit solvents. In this endeavor, we aim to establish model parametrization and enhance computational efficiency. Unlike traditional CG models that use empirical force parameters, ex-force parameters (r0(ex), a ~ , δd, δp) of DPD particles constructed for specific research purposes can be obtained from atomistic molecular dynamics simulations. On the other hand, im-force parameters (r0(im), c, σ) can be derived from ex-DPD simulations, according to the underlying thermodynamic theory. Based on a mapping scheme proposed for the modeling of amino acids, all-atom proteins can be converted into a CG model. Both ex-/im-DPDs are then carried out to investigate the folding pathways of the four mini-proteins. Structural analysis of the RMSDs shows that the im-simulated proteins have greater structural similarity to native proteins than the ex-simulated ones. The constructed CG models achieve a resolution of Angstrom (Å), a level normally associated with atomic models. Additionally, speed tests reveal that im-DPD accelerates the simulation process and significantly improves simulation efficiency.

Protein Retrieval via Integrative Molecular Ensembles (PRIME) through Extended Similarity Indices

Molecular dynamics (MD) simulations are ideally suited to describe conformational ensembles of biomolecules such as proteins and nucleic acids. Microsecond-long simulations are now routine, facilitated by the emergence of graphical processing units. Clustering, which groups objects based on structural similarity, is typically used to process ensembles, leading to different states, their populations, and the identification of representative structures. A popular pipeline combines hierarchical clustering for clustering and selecting the cluster centroid as representative of the cluster. Here, we propose to improve on this approach, by developing a module-Protein Retrieval via Integrative Molecular Ensembles (PRIME), that consists of tools to improve the prediction of the representative in the most populated cluster using extended continuous similarity. PRIME is integrated with our Molecular Dynamics Analysis with N-ary Clustering Ensembles (MDANCE) package and can be used as a postprocessing tool for arbitrary clustering algorithms, compatible with several MD suites. PRIME predictions produced structures that when aligned to the experimental structure were better superposed (lower RMSD). A further benefit of PRIME is its linear scaling─rather than the traditional O(N2) traditionally associated with comparisons of elements in a set.

Inclusion of Water Multipoles into the Implicit Solvation Framework Leads to Accuracy Gains

The current practical "workhorses" of the atomistic implicit solvation─the Poisson-Boltzmann (PB) and generalized Born (GB) models─face fundamental accuracy limitations. Here, we propose a computationally efficient implicit solvation framework, the Implicit Water Multipole GB (IWM-GB) model, that systematically incorporates the effects of multipole moments of water molecules in the first hydration shell of a solute, beyond the dipole water polarization already present at the PB/GB level. The framework explicitly accounts for coupling between polar and nonpolar contributions to the total solvation energy, which is missing from many implicit solvation models. An implementation of the framework, utilizing the GAFF force field and AM1-BCC atomic partial charges model, is parametrized and tested against the experimental hydration free energies of small molecules from the FreeSolv database. The resulting accuracy on the test set (RMSE ∼ 0.9 kcal/mol) is 12% better than that of the explicit solvation (TIP3P) treatment, which is orders of magnitude slower. We also find that the coupling between polar and nonpolar parts of the solvation free energy is essential to ensuring that several features of the IWM-GB model are physically meaningful, including the sign of the nonpolar contributions.

Exact Analytical Algorithm for the Solvent-Accessible Surface Area and Derivatives in Implicit Solvent Molecular Simulations on GPUs

In this paper, we present differentiable solvent-accessible surface area (dSASA), an exact geometric method to calculate SASA analytically along with atomic derivatives on GPUs. The atoms in a molecule are first assigned to tetrahedra in groups of four atoms by Delaunay tetrahedralization adapted for efficient GPU implementation, and the SASA values for atoms and molecules are calculated based on the tetrahedralization information and inclusion-exclusion method. The SASA values from the numerical icosahedral-based method can be reproduced with >98% accuracy for both proteins and RNAs. Having been implemented on GPUs and incorporated into AMBER, we can apply dSASA to implicit solvent molecular dynamics simulations with the inclusion of this nonpolar term. The current GPU version of GB/SA simulations has been accelerated up to nearly 20-fold compared to the CPU version, outperforming LCPO, a commonly used, fast algorithm for calculating SASA, as the system size increases. While we focus on the accuracy of the SASA calculations for proteins and nucleic acids, we also demonstrate stable GB/SA MD mini-protein simulations.

Assessing the Performance of Peptide Force Fields for Modeling the Solution Structural Ensembles of Cyclic Peptides

Molecular dynamics simulation is a powerful tool for characterizing the solution structural ensembles of cyclic peptides. However, the ability of simulation to recapitulate experimental results and make accurate predictions largely depends on the force fields used. In our work here, we evaluate the performance of seven state-of-the-art force fields in recapitulating the experimental NMR results in water of 12 benchmark cyclic peptides, consisting of 6 cyclic pentapeptides, 4 cyclic hexapeptides, and 2 cyclic heptapeptides. The results show that RSFF2+TIP3P, RSFF2C+TIP3P, and Amber14SB+TIP3P exhibit similar and the best performance, all recapitulating the NMR-derived structure information on 10 cyclic peptides. Amber19SB+OPC successfully recapitulates the NMR-derived structure information on 8 cyclic peptides. In contrast, OPLS-AA/M+TIP4P, Amber03+TIP3P, and Amber14SBonlysc+GB-neck2 could only recapitulate the NMR-derived structure information on 5 cyclic peptides, the majority of which are not well-structured.

Chimeric nanobody-decorated liposomes by self-assembly

Liposomes as drug vehicles have advantages, such as payload protection, tunable carrying capacity and improved biodistribution. However, due to the dysfunction of targeting moieties and payload loss during preparation, immunoliposomes have yet to be favoured in commercial manufacturing. Here we report a chemical modification-free biophysical approach for producing immunoliposomes in one step through the self-assembly of a chimeric nanobody (cNB) into liposome bilayers. cNB consists of a nanobody against human epidermal growth factor receptor 2 (HER2), a flexible peptide linker and a hydrophobic single transmembrane domain. We determined that 64% of therapeutic compounds can be encapsulated into 100-nm liposomes, and up to 2,500 cNBs can be anchored on liposomal membranes without steric hindrance under facile conditions. Subsequently, we demonstrate that drug-loaded immunoliposomes increase cytotoxicity on HER2-overexpressing cancer cell lines by 10- to 20-fold, inhibit the growth of xenograft tumours by 3.4-fold and improve survival by more than twofold.

Exploring Structural Insights of Aβ42 and α-Synuclein Monomers and Heterodimer: A Comparative Study Using Implicit and Explicit Solvent Simulations

Protein misfolding, aggregation, and fibril formation play a central role in the development of severe neurological disorders, including Alzheimer's and Parkinson's diseases. The structural stability of mature fibrils in these diseases is of great importance, as organisms struggle to effectively eliminate amyloid plaques. To address this issue, it is crucial to investigate the early stages of fibril formation when monomers aggregate into small, toxic, and soluble oligomers. However, these structures are inherently disordered, making them challenging to study through experimental approaches. Recently, it has been shown experimentally that amyloid-β 42 (Aβ42) and α-synuclein (α-Syn) can coassemble. This has motivated us to investigate the interaction between their monomers as a first step toward exploring the possibility of forming heterodimeric complexes. In particular, our study involves the utilization of various Amber and CHARMM force-fields, employing both implicit and explicit solvent models in replica exchange and conventional simulation modes. This comprehensive approach allowed us to assess the strengths and weaknesses of these solvent models and force fields in comparison to experimental and theoretical findings, ensuring the highest level of robustness. Our investigations revealed that Aβ42 and α-Syn monomers can indeed form stable heterodimers, and the resulting heterodimeric model exhibits stronger interactions compared to the Aβ42 dimer. The binding of α-Syn to Aβ42 reduces the propensity of Aβ42 to adopt fibril-prone conformations and induces significant changes in its conformational properties. Notably, in AMBER-FB15 and CHARMM36m force fields with the use of explicit solvent, the presence of Aβ42 significantly increases the β-content of α-Syn, consistent with the experiments showing that Aβ42 triggers α-Syn aggregation. Our analysis clearly shows that although the use of implicit solvent resulted in too large compactness of monomeric α-Syn, structural properties of monomeric Aβ42 and the heterodimer were preserved in explicit-solvent simulations. We anticipate that our study sheds light on the interaction between α-Syn and Aβ42 proteins, thus providing the atom-level model required to assess the initial stage of aggregation mechanisms related to Alzheimer's and Parkinson's diseases.

Exploring the inhibitory potential of KPHs-AL-derived GLLF peptide on pancreatic lipase and cholesterol esterase activities

The effective modulation of pancreatic lipase and cholesterol esterase activities proves critical in maintaining circulatory triglycerides and cholesterol levels within physiological boundaries. In this study, peptides derived from KPHs-AL, produced through the enzymatic hydrolysis of skipjack tuna dark muscle using alkaline protease, have a specific inhibitory effect on pancreatic lipase and cholesterol esterase. It is hypothesized that these peptides target and modulate the activities of enzymes by inducing conformational changes within their binding pockets, potentially impacting the catalytic functions of both pancreatic lipase and cholesterol esterase. Results revealed these peptides including AINDPFIDL, FLGM, GLLF and WGPL, were found to nestle into the binding site groove of pancreatic lipase and cholesterol esterase. Among these, GLLF stood out, demonstrating potent inhibition with IC50 values of 0.1891 mg/mL and 0.2534 mg/mL for pancreatic lipase and cholesterol esterase, respectively. The kinetics studies suggested that GLLF competed effectively with substrates for the enzyme active sites. Spectroscopic analyses, including ultraviolet-visible, fluorescence quenching, and circular dichroism, indicated that GLLF binding induced conformational changes within the enzymes, likely through hydrogen bond formation and hydrophobic interactions, thereby increasing structural flexibility. Molecular docking and molecular dynamics simulations supported these findings, showing GLLF's stable interaction with vital active site residues. These findings position GLLF as a potent inhibitor of key digestive enzymes, offering insights into its role in regulating lipid metabolism and highlighting its potential as functional ingredient.

Exact analytical algorithm for solvent accessible surface area and derivatives in implicit solvent molecular simulations on GPUs

In this paper, we present dSASA (differentiable SASA), an exact geometric method to calculate solvent accessible surface area (SASA) analytically along with atomic derivatives on GPUs. The atoms in a molecule are first assigned to tetrahedra in groups of four atoms by Delaunay tetrahedrization adapted for efficient GPU implementation and the SASA values for atoms and molecules are calculated based on the tetrahedrization information and inclusion-exclusion method. The SASA values from the numerical icosahedral-based method can be reproduced with more than 98% accuracy for both proteins and RNAs. Having been implemented on GPUs and incorporated into the software Amber, we can apply dSASA to implicit solvent molecular dynamics simulations with inclusion of this nonpolar term. The current GPU version of GB/SA simulations has been accelerated up to nearly 20-fold compared to the CPU version, outperforming LCPO, a commonly used, fast algorithm for calculating SASA, as the system size increases. While we focus on the accuracy of the SASA calculations for proteins and nucleic acids, we also demonstrate stable GB/SA MD mini-protein simulations.

Myomerger Induces Membrane Hemifusion and Regulates Fusion Pore Expansion

Membrane fusion is a crucial mechanism in a wide variety of important events in cell biology from viral infection to exocytosis. However, despite many efforts and much progress, cell-cell fusion has remained elusive to our understanding. Along the life of the fusion pore, large conformational changes take place from the initial lipid bilayer bending, passing through the hemifusion intermediates, and ending with the formation of the first nascent fusion pore. In this sense, computer simulations are an ideal technique for describing such complex lipid remodeling at the molecular level. In this work, we studied the role played by the muscle-specific membrane protein Myomerger during the formation of the fusion pore. We have conducted μs length atomistic and coarse-grained molecular dynamics, together with free-energy calculations using ad hoc collective variables. Our results show that Myomerger favors the hemifusion diaphragm-stalk transition, reduces the nucleation-expansion energy difference, and promotes the formation of nonenlarging fusion pores.

Oligo-PROTAC strategy for cell-selective and targeted degradation of activated STAT3

Decoy oligodeoxynucleotides (ODNs) allow targeting undruggable transcription factors, such as STAT3, but their limited potency and lack of delivery methods hampered translation. To overcome these challenges, we conjugated a STAT3-specific decoy to thalidomide, a ligand to cereblon in E3 ubiquitin ligase complex, to generate a proteolysis-targeting chimera (STAT3DPROTAC). STAT3DPROTAC downregulated STAT3 in target cells, but not STAT1 or STAT5. Computational modeling of the STAT3DPROTAC ternary complex predicted two surface lysines, K601 and K626, in STAT3 as potential ubiquitination sites. Accordingly, K601/K626 point mutations in STAT3, as well as proteasome inhibition or cereblon deletion, alleviated STAT3DPROTAC effect. Next, we conjugated STAT3DPROTAC to a CpG oligonucleotide targeting Toll-like receptor 9 (TLR9) to generate myeloid/B cell-selective C-STAT3DPROTAC. Naked C-STAT3DPROTAC was spontaneously internalized by TLR9+ myeloid cells, B cells, and human and mouse lymphoma cells but not by T cells. C-STAT3DPROTAC effectively decreased STAT3 protein levels and also STAT3-regulated target genes critical for lymphoma cell proliferation and/or survival (BCL2L1, CCND2, and MYC). Finally, local C-STAT3DPROTAC administration to human Ly3 lymphoma-bearing mice triggered tumor regression, while control C-STAT3D and C-SCR treatments had limited effects. Our results underscore the feasibility of using a PROTAC strategy for cell-selective, decoy oligonucleotide-based STAT3 targeting of and potentially other tumorigenic transcription factors for cancer therapy.

Improving the Accuracy of Physics-Based Hydration-Free Energy Predictions by Machine Learning the Remaining Error Relative to the Experiment

The accuracy of computational models of water is key to atomistic simulations of biomolecules. We propose a computationally efficient way to improve the accuracy of the prediction of hydration-free energies (HFEs) of small molecules: the remaining errors of the physics-based models relative to the experiment are predicted and mitigated by machine learning (ML) as a postprocessing step. Specifically, the trained graph convolutional neural network attempts to identify the "blind spots" in the physics-based model predictions, where the complex physics of aqueous solvation is poorly accounted for, and partially corrects for them. The strategy is explored for five classical solvent models representing various accuracy/speed trade-offs, from the fast analytical generalized Born (GB) to the popular TIP3P explicit solvent model; experimental HFEs of small neutral molecules from the FreeSolv set are used for the training and testing. For all of the models, the ML correction reduces the resulting root-mean-square error relative to the experiment for HFEs of small molecules, without significant overfitting and with negligible computational overhead. For example, on the test set, the relative accuracy improvement is 47% for the fast analytical GB, making it, after the ML correction, almost as accurate as uncorrected TIP3P. For the TIP3P model, the accuracy improvement is about 39%, bringing the ML-corrected model's accuracy below the 1 kcal/mol threshold. In general, the relative benefit of the ML corrections is smaller for more accurate physics-based models, reaching the lower limit of about 20% relative accuracy gain compared with that of the physics-based treatment alone. The proposed strategy of using ML to learn the remaining error of physics-based models offers a distinct advantage over training ML alone directly on reference HFEs: it preserves the correct overall trend, even well outside of the training set.

WLJP-025p, a homogeneous Lonicera japonica polysaccharide, attenuates atopic dermatitis by regulating the MAPK/NFκB/AP-1 axis via Act1

Atopic dermatitis (AD) is usually treated with steroids, but long-term use is not an effective cure because side effects and disease aggravation. Therefore, more effective and safer treatments are needed. Using dexamethasone as the positive control, the mechanism of action of water-extracted homogeneous honeysuckle Lonicera japonica polysaccharide (WLJP-025p) to alleviate AD was studied. Mice were administered 2,4-dinitrochlorobenzene in their bare back and right ear to mimic an AD model. The efficacy of WLJP-025p in AD was assessed by measuring right ear thickness and skin lesion scores, pathological observation (haematoxylin-eosin and toluidine blue staining), and serum IgE and IL-1β concentrations. The expression of relevant genes and proteins in the serum and back skin was detected using RT-qPCR, ELISA, western blotting, and immunofluorescence. Molecular docking and dynamic simulation of WLJP-025p and Act1 were performed. WLJP-025p considerably alleviated skin hyperplasia and pathological abnormalities in AD mice and inhibited the expression of Act1, Nucleus-P65, Nucleus-AP-1, and MAPK-related proteins in skin tissues. WLJP-025p formed a stable conformation with Act1, inhibited splenic Th17 differentiation, IL-17 release, and upregulated the expression of related skin barrier proteins. In conclusion, WLJP-025p affects the inflammation regulation via the MAPK/NFκB/AP-1 axis by binding to Act1, promotes the recovery of epithelial barrier function, and alleviates AD in mice.

When MELD Meets GaMD: Accelerating Biomolecular Landscape Exploration

We introduce Gaussian accelerated MELD (GaMELD) as a new method for exploring the energy landscape of biomolecules. GaMELD combines the strengths of Gaussian accelerated molecular dynamics (GaMD) and modeling employing limited data (MELD) to navigate complex energy landscapes. MELD uses replica-exchange molecular simulations to integrate limited and uncertain data into simulations via Bayesian inference. MELD has been successfully applied to problems of structure prediction like protein folding and complex structure prediction. However, the computational cost for MELD simulations has limited its broader applicability. The synergy of GaMD and MELD surmounts this limitation efficiently sampling the energy landscape at a lower computational cost (reducing the computational cost by a factor of 2 to six). Effectively, GaMD is used to shift energy distributions along replicas to increase the overlap in energy distributions across replicas, facilitating a random walk in replica space. We tested GaMELD on a benchmark set of 12 small proteins that have been previously studied through MELD and conventional MD. GaMELD consistently achieves accurate predictions with fewer replicas. By increasing the efficacy of replica exchange, GaMELD effectively accelerates convergence in the conformational space, enabling improved sampling across a diverse set of systems.

Mapping the configurational landscape and aggregation phase behavior of the tau protein fragment PHF6

The PHF6 (Val-Gln-Ile-Val-Tyr-Lys) motif, found in all isoforms of the microtubule-associated protein tau, forms an integral part of ordered cores of amyloid fibrils formed in tauopathies and is thought to play a fundamental role in tau aggregation. Because PHF6 as an isolated hexapeptide assembles into ordered fibrils on its own, it is investigated as a minimal model for insight into the initial stages of aggregation of larger tau fragments. Even for this small peptide, however, the large length and time scales associated with fibrillization pose challenges for simulation studies of its dynamic assembly, equilibrium configurational landscape, and phase behavior. Here, we develop an accurate, bottom-up coarse-grained model of PHF6 for large-scale simulations of its aggregation, which we use to uncover molecular interactions and thermodynamic driving forces governing its assembly. The model, not trained on any explicit information about fibrillar structure, predicts coexistence of formed fibrils with monomers in solution, and we calculate a putative equilibrium phase diagram in concentration-temperature space. We also characterize the configurational and free energetic landscape of PHF6 oligomers. Importantly, we demonstrate with a model of heparin that this widely studied cofactor enhances the aggregation propensity of PHF6 by ordering monomers during nucleation and remaining associated with growing fibrils, consistent with experimentally characterized heparin-tau interactions. Overall, this effort provides detailed molecular insight into PHF6 aggregation thermodynamics and pathways and, furthermore, demonstrates the potential of modern multiscale modeling techniques to produce predictive models of amyloidogenic peptides simultaneously capturing sequence-specific effects and emergent aggregate structures.

Rapid Rescoring and Refinement of Ligand-Receptor Complexes Using Replica Exchange Molecular Dynamics with a Monte Carlo Pose Reservoir

Virtual screening (VS) involves generation of poses for a library of ligands and ranking using simplified energy functions and limited flexibility. Top-scored poses are used to rank and prioritize ligands. Here, we adapt the reservoir replica exchange molecular dynamics (res-REMD) method to rerank poses generated through VS. REMD simulations are carried out but with occasional Monte Carlo jumps to alternate VS-generated poses using a Metropolis criterion. The simulations converge within 10 ns for all systems, generating populations of alternate poses in the context of fully flexible ligand and protein side chains. The protocol is applied to four model protein-ligand complexes, where DOCK resulted in two successes and two scoring failures. In all four systems, the most populated cluster from the final ensemble exhibits high similarity to the crystallographic pose with ligand RMSD values under 2.0 Å. Both DOCK failures were rescued. For one DOCK success, the protocol identified the correct pose but also sampled an alternate pose at equal probability. Opportunities for future improvements and extensions are discussed.

Improving Estimation of the Koopman Operator with Kolmogorov-Smirnov Indicator Functions

It has become common to perform kinetic analysis using approximate Koopman operators that transform high-dimensional timeseries of observables into ranked dynamical modes. The key to the practical success of the approach is the identification of a set of observables that form a good basis on which to expand the slow relaxation modes. Good observables are, however, difficult to identify a priori and suboptimal choices can lead to significant underestimations of characteristic time scales. Leveraging the representation of slow dynamics in terms of Hidden Markov Models (HMM), we propose a simple and computationally efficient clustering procedure to infer surrogate observables that form a good basis for slow modes. We apply the approach to an analytically solvable model system as well as on three protein systems of different complexities. We consistently demonstrate that the inferred indicator functions can significantly improve the estimation of the leading eigenvalues of Koopman operators and correctly identify key states and transition time scales of stochastic systems, even when good observables are not known a priori.

Intrinsic disorder and conformational coexistence in auxin coreceptors

AUXIN/INDOLE 3-ACETIC ACID (Aux/IAA) transcriptional repressor proteins and the TRANSPORT INHIBITOR RESISTANT 1/AUXIN SIGNALING F-BOX (TIR1/AFB) proteins to which they bind act as auxin coreceptors. While the structure of TIR1 has been solved, structural characterization of the regions of the Aux/IAA protein responsible for auxin perception has been complicated by their predicted disorder. Here, we use NMR, CD and molecular dynamics simulation to investigate the N-terminal domains of the Aux/IAA protein IAA17/AXR3. We show that despite the conformational flexibility of the region, a critical W-P bond in the core of the Aux/IAA degron motif occurs at a strikingly high (1:1) ratio of cis to trans isomers, consistent with the requirement of the cis conformer for the formation of the fully-docked receptor complex. We show that the N-terminal half of AXR3 is a mixture of multiple transiently structured conformations with a propensity for two predominant and distinct conformational subpopulations within the overall ensemble. These two states were modeled together with the C-terminal PB1 domain to provide the first complete simulation of an Aux/IAA. Using MD to recreate the assembly of each complex in the presence of auxin, both structural arrangements were shown to engage with the TIR1 receptor, and contact maps from the simulations match closely observations of NMR signal-decreases. Together, our results and approach provide a platform for exploring the functional significance of variation in the Aux/IAA coreceptor family and for understanding the role of intrinsic disorder in auxin signal transduction and other signaling systems.

Einstein Model of a Graph to Characterize Protein Folded/Unfolded States

The folded structures of proteins can be accurately predicted by deep learning algorithms from their amino-acid sequences. By contrast, in spite of decades of research studies, the prediction of folding pathways and the unfolded and misfolded states of proteins, which are intimately related to diseases, remains challenging. A two-state (folded/unfolded) description of protein folding dynamics hides the complexity of the unfolded and misfolded microstates. Here, we focus on the development of simplified order parameters to decipher the complexity of disordered protein structures. First, we show that any connected, undirected, and simple graph can be associated with a linear chain of atoms in thermal equilibrium. This analogy provides an interpretation of the usual topological descriptors of a graph, namely the Kirchhoff index and Randić resistance, in terms of effective force constants of a linear chain. We derive an exact relation between the Kirchhoff index and the average shortest path length for a linear graph and define the free energies of a graph using an Einstein model. Second, we represent the three-dimensional protein structures by connected, undirected, and simple graphs. As a proof of concept, we compute the topological descriptors and the graph free energies for an all-atom molecular dynamics trajectory of folding/unfolding events of the proteins Trp-cage and HP-36 and for the ensemble of experimental NMR models of Trp-cage. The present work shows that the local, nonlocal, and global force constants and free energies of a graph are promising tools to quantify unfolded/disordered protein states and folding/unfolding dynamics. In particular, they allow the detection of transient misfolded rigid states.

Oligo-PROTAC strategy for cell-selective and targeted degradation of activated STAT3

Decoy-oligodeoxynucleotides (D-ODNs) can target undruggable transcription factors, such as STAT3. However, challenges in D-ODN delivery and potency hampered their translation. To overcome these limitations, we conjugated STAT3-specific D-ODN to thalidomide (Tha), a known ligand to cereblon (CRBN, a component of E3 ubiquitin ligase) to generate a proteolysis-targeting chimera (STAT3D PROTAC ). STAT3D PROTAC downregulated STAT3, but not STAT1 or STAT5, in target cells. Computational modeling of the STAT3D PROTAC ternary complex predicted two surface lysines on STAT3, K601 and K626 as potential ubiquitination sites for the PROTAC bound E3 ligase. Accordingly, K601/K626 point mutations in STAT3, as well as proteasome inhibitors, and CRBN deletion alleviated STAT3D PROTAC effect. Next, we conjugated STAT3D PROTAC to a CpG ligand targeting Toll-like receptor 9 (TLR9) to generate myeloid/B-cell-selective C-STAT3D PROTAC conjugate. Naked C-STAT3D PROTAC was spontaneously internalized by TLR9 + myeloid cells, B cells as well as human Ly18 and mouse A20 lymphoma cells, but not by T cells. C-STAT3D PROTAC decreased STAT3 levels to 50% at 250 nM and over 85% at 2 µM dosing in myeloid cells. We also observed significantly improved downregulation of STAT3 target genes involved in lymphoma cell proliferation and/or survival ( BCL2L1, CCND2, MYC ). Finally, we assessed the antitumor efficacy of C-STAT3D PROTAC compared to C-STAT3D or scrambled control (C-SCR) against human lymphoma xenotransplants. Local C-STAT3D PROTAC administration triggered lymphoma regression while control treatments had limited effects. Our results underscore feasibility of using PROTAC strategy for cell-selective, decoy oligonucleotide-based targeting of STAT3 and potentially other tumorigenic transcription factors for cancer therapy.

Novel Insights on 6:6 Perfluoroalkyl Phosphonic Acid-Induced Melanin Synthesis Disorders Leading to Pigmentation in Tadpoles

As alternatives to traditional per- and polyfluoroalkyl substances, perfluoroalkyl phosphonic acids (PFPiAs) are widely present in aquatic environments and can potentially harm aquatic organisms. Pigmentation affects the probability of aquatic organisms being preyed on and serves as an important toxic endpoint of development, but little is known about the impacts of PFPiAs on the development of aquatic organisms. In this study, Xenopus laevis embryos were exposed to 6:6 PFPiA (1, 10, and 100 nM) for 14 days. The developed tadpoles exhibited evident pigmentation with increased melanin particle size and density on the skin. Pathological and behavioral experiments revealed that the retinal layers became thinner, reducing the photosensitivity and disturbing the circadian rhythm of the tadpoles. Compared to the control group, the exposed tadpoles showed higher levels but less changes of melanin throughout the light/dark cycle, as well as distinct oxidative damage. Consequently, the expression level of microphthalmia-associated transcription factor (MITF), a key factor inducing melanin synthesis, increased significantly. Molecular docking analysis suggested that 6:6 PFPiA forms strong interactions in the binding pocket of MITF, implying that it could activate MITF directly. The activation of MITF ultimately promotes melanin synthesis, resulting in pigmentation on tadpoles.

Enhanced sampling in explicit solvent by deep learning module in FSATOOL

FSATOOL is an integrated molecular simulation and data analysis program. Its old molecular dynamics engine only supports simulations in vacuum or implicit solvent. In this work, we implement the well-known smooth particle mesh Ewald method for simulations in explicit solvent. The new developed engine is runnable on both CPU and GPU. All the existed analysis modules in the program are compatible with the new engine. Moreover, we also build a complete deep learning module in FSATOOL. Based on the module, we further implement two useful trajectory analysis methods: state-free reversible VAMPnets and time-lagged autoencoder. They are good at searching the collective variables related to the conformational transitions of biomolecules. In FSATOOL, these collective variables can be further used to construct a bias potential for the enhanced sampling purpose. We introduce the implementation details of the methods and present their actual performances in FSATOOL by a few enhanced sampling simulations.

Rethinking the protein folding problem from a new perspective

One of the main concerns of Anfinsen was to reveal the connection between the amino-acid sequence and their biologically active conformation. This search gave rise to two crucial questions in structural biology, namely, why the proteins fold and how a sequence encodes its folding. As to the why, he proposes a plausible answer, namely, the thermodynamic hypothesis. As to the how, this remains an unsolved challenge. Consequently, the protein folding problem is examined here from a new perspective, namely, as an 'analytic whole'. Conceiving the protein folding in this way enabled us to (i) examine in detail why the force-field-based approaches have failed, among other purposes, in their ability to predict the three-dimensional structure of a protein accurately; (ii) propose how to redefine them to prevent these shortcomings, and (iii) conjecture on the origin of the state-of-the-art numerical-methods success to predict the tridimensional structure of proteins accurately.

Basis for Accurate Protein pKa Prediction with Machine Learning

pH regulates protein structures and the associated functions in many biological processes via protonation and deprotonation of ionizable side chains where the titration equilibria are determined by pK_a's. To accelerate pH-dependent molecular mechanism research in the life sciences or industrial protein and drug designs, fast and accurate pK_a prediction is crucial. Here we present a theoretical pK_a data set PHMD549, which was successfully applied to four distinct machine learning methods, including DeepKa, which was proposed in our previous work. To reach a valid comparison, EXP67S was selected as the test set. Encouragingly, DeepKa was improved significantly and outperforms other state-of-the-art methods, except for the constant-pH molecular dynamics, which was utilized to create PHMD549. More importantly, DeepKa reproduced experimental pK_a orders of acidic dyads in five enzyme catalytic sites. Apart from structural proteins, DeepKa was found applicable to intrinsically disordered peptides. Further, in combination with solvent exposures, it is revealed that DeepKa offers the most accurate prediction under the challenging circumstance that hydrogen bonding or salt bridge interaction is partly compensated by desolvation for a buried side chain. Finally, our benchmark data qualify PHMD549 and EXP67S as the basis for future developments of protein pK_a prediction tools driven by artificial intelligence. In addition, DeepKa built on PHMD549 has been proven an efficient protein pK_a predictor and thus can be applied immediately to, for example, pK_a database construction, protein design, drug discovery, and so on.

CryoFold 2.0: Cryo-EM Structure Determination with MELD

Cryo-electron microscopy data are becoming more prevalent and accessible at higher resolution levels, leading to the development of new computational tools to determine the atomic structure of macromolecules. However, while existing tools adapted from X-ray crystallography are suitable for the highest-resolution maps, new tools are needed for lower-resolution levels and to account for map heterogeneity. In this article, we introduce CryoFold 2.0, an integrative physics-based approach that combines Bayesian inference and the ability to handle multiple data sources with the molecular dynamics flexible fitting (MDFF) approach to determine the structures of macromolecules by using cryo-EM data. CryoFold 2.0 is incorporated into the MELD (modeling employing limited data) plugin, resulting in a pipeline that is more computationally efficient and accurate than running MELD or MDFF alone. The approach requires fewer computational resources and shorter simulation times than the original CryoFold, and it minimizes manual intervention. We demonstrate the effectiveness of the approach on eight different systems, highlighting its various benefits.

Exploring the Transferability of Replica Exchange Structure Reservoirs to Accelerate Generation of Ensembles for Alternate Hamiltonians or Protein Mutations

Generating precise ensembles is commonly a prerequisite to understand the energetics of biological processes using Molecular Dynamics (MD) simulations. Previously, we have shown how unweighted reservoirs built from high temperature MD simulations can accelerate convergence of Boltzmann-weighted ensembles by at least 10× with the Reservoir Replica Exchange MD (RREMD) method. Therefore, in this work, we explore whether an unweighted structure reservoir generated with one Hamiltonian (solute force field plus solvent model) can be reused to quickly generate accurately weighted ensembles for Hamiltonians other than the one that was used to generate the reservoir. We also extended this methodology to rapidly estimate the effects of mutations on peptide stability by using a reservoir of diverse structures obtained from wild-type simulations. These results suggest that structures generated via fast methods such as coarse-grained models or structures predicted by Rosetta or deep learning approaches could be integrated into a reservoir to accelerate generation of ensembles using more accurate representations.

Integrative network pharmacology and in silico analyses identify the anti-omicron SARS-CoV-2 potential of eugenol

Eugenol as a natural product is the source of isoniazid, and purified eugenol is extensively used in the cosmetics industry and the productive processes of edible spices. Accumulating evidence suggested that eugenol exerted potent anti-microorganism and anti-inflammation effects. Application of eugenol effectively reduced the risk of atherosclerosis, arterial embolism, and Type 2 diabetes. A previous study confirmed that treatment with eugenol attenuated lung inflammation and improved heart functions in SARS-CoV-2 spike S1-intoxicated mice. In addition to the study, based on a series of public datasets, computational analyses were conducted to characterize the acting targets of eugenol and the functional roles of these targets in COVID-19. The binding capacities of eugenol to conservative sites of SARS-CoV-2 like RNA-dependent RNA polymerase (RdRp) and mutable site as spike (S) protein, were calculated by using molecular docking following the molecular dynamics simulation with RMSD, RMSF, and MM-GBSA methods. The results of network pharmacology indicated that six targets, including PLAT, HMOX1, NUP88, CTSL, ITGB1 andTMPRSS2 were eugenol-SARS-CoV-2 interacting proteins. The omics results of in-silico study further implicated that eugenol increased the expression of SCARB1, HMOX1 and GDF15, especially HMOX1, which were confirmed the potential interacting targets between eugenol and SARS-CoV-2 antigens. Enrichment analyses indicated that eugenol exerted extensive biological effects such as regulating immune infiltration of macrophage, lipid localization, monooxyenase activity, iron ion binding and PPAR signaling. The results of the integrated analysis of eugenol targets and immunotranscription profile of COVID-19 cases shows that eugenol also plays an important role in strengthen of immunologic functions and regulating cytokine signaling. As a complement to the integrated analysis, the results of molecular docking indicated the potential binding interactions between eugenol and four proteins relating to cytokine production/release and the function of T type lymphocytes, including human TLR-4, TCR, NF-κB, JNK and AP-1. Furthermore, results of molecular docking and molecular dynamics (100ns) simulations implicated that stimulated modification of eugenol to the SARS-CoV-2 Omicron Spike-ACE2 complex, especially for human ACE2, and the molecular interaction of eugenol to SARS-CoV-2 RdRp, were no less favorable than two positive controls, molnupiravir and nilotinib. Dynamics (200ns) simulations indicated that the binding capacities and stabilities of eugenol to finger subdomain of RdRp is no less than molnupiravir. However, the simulated binding capacity of eugenol to SARS-CoV-2 wild type RBD and Omicron mutant RBD were less than nilotinib. Eugenol was predicted to have more favor LD50 value and lower cytotoxicity than two positive controls, and eugenol can pass through the blood-brain barrier (BBB). In a brief, eugenol is helpful for attenuating systemic inflammation induced by SARS-CoV-2 infection, due to the direct interaction of eugenol to SARS-CoV-2 proteins and extensive bio-manipulation of pro-inflammatory factors. This study carefully suggests eugenol is a candidate compound of developing drugs and supplement agents against SARS-CoV-2 and its Omicron variants.

Computational modeling of the molecular basis for the calcium-dependence of the mannuronan C-5 epimerase AvAlgE6 from Azotobacter vinelandii

The mannuronan C-5 epimerases catalyze epimerization of β-d-mannuronic acid to α-l-guluronic acid in alginate polymers. The seven extracellular Azotobacter vinelandii epimerases (AvAlgE1-7) are calcium-dependent, and calcium is essential for the structural integrity of their carbohydrate binding R-modules. Ca²⁺ is also found in the crystal structures of the A-modules, where it is suggested to play a structural role. In this study, the structure of the catalytic A-module of the A. vinelandii mannuronan C-5 epimerase AvAlgE6 is used to investigate the role of this Ca²⁺. Molecular dynamics (MD) simulations with and without calcium reveal the possible importance of the bound Ca²⁺ in the hydrophobic packing of β-sheets. In addition, a putative calcium binding site is found in the active site, indicating a potential direct role of this calcium in the catalysis. According to the literature, two of the residues coordinating calcium in this site are essential for the activity. MD simulations of the interaction with bound substrate indicate that the presence of a calcium ion in this binding site increases the binding strength. Further, explicit calculations of the substrate dissociation pathways with umbrella sampling simulations show and energetically higher dissociation barrier when calcium is present. The present study eludes to a putative catalytic role of calcium in the charge neutralizing first step of the enzymatic reaction. In addition to the importance for understanding these enzymes' molecular mechanisms, this could have implications for engineering strategies of the epimerases in industrial alginate processing.

On the plasticity of amyloid formation: The impact of destabilizing small to large substitutions on islet amyloid polypeptide amyloid formation

Amyloids are partially ordered, proteinaceous, β-sheet rich deposits that have been implicated in a wide range of diseases. An even larger set of proteins that do not normally form amyloid in vivo can be induced to do so in vitro. A growing number of structures of amyloid fibrils have been reported and a common feature is the presence of a tightly packed core region in which adjacent monomers pack together in extremely tight interfaces, often referred to as steric zippers. A second common feature of many amyloid fibrils is their polymorphous nature. We examine the consequences of disrupting the tight packing in amyloid fibrils on the kinetics of their formation using the 37 residue polypeptide hormone islet amyloid polypeptide (IAPP, amylin) as a model system. IAPP forms islet amyloid in vivo and is aggressively amyloidogenic in vitro. Six Cryo-EM structures of IAPP amyloid fibrils are available and in all Gly24 is in the core of the structured region and makes tight contacts with other residues. Calculations using the ff14SBonlysc forcefield in Amber20 show that substitutions with larger amino acids significantly disrupt close packing and are predicted to destabilize the various fibril structures. However, Gly to 2-amino butyric acid (2-carbon side chain) and Gly to Leu substitutions actually enhance the rate of amyloid formation. A Pro substitution slows, but does not prevent amyloid formation.

A Suite of Tutorials for the WESTPA 2.0 Rare-Events Sampling Software [Article v2.0]

The weighted ensemble (WE) strategy has been demonstrated to be highly efficient in generating pathways and rate constants for rare events such as protein folding and protein binding using atomistic molecular dynamics simulations. Here we present two sets of tutorials instructing users in the best practices for preparing, carrying out, and analyzing WE simulations for various applications using the WESTPA software. The first set of more basic tutorials describes a range of simulation types, from a molecular association process in explicit solvent to more complex processes such as host-guest association, peptide conformational sampling, and protein folding. The second set ecompasses six advanced tutorials instructing users in the best practices of using key new features and plugins/extensions of the WESTPA 2.0 software package, which consists of major upgrades for larger systems and/or slower processes. The advanced tutorials demonstrate the use of the following key features: (i) a generalized resampler module for the creation of "binless" schemes, (ii) a minimal adaptive binning scheme for more efficient surmounting of free energy barriers, (iii) streamlined handling of large simulation datasets using an HDF5 framework, (iv) two different schemes for more efficient rate-constant estimation, (v) a Python API for simplified analysis of WE simulations, and (vi) plugins/extensions for Markovian Weighted Ensemble Milestoning and WE rule-based modeling for systems biology models. Applications of the advanced tutorials include atomistic and non-spatial models, and consist of complex processes such as protein folding and the membrane permeability of a drug-like molecule. Users are expected to already have significant experience with running conventional molecular dynamics or systems biology simulations.

Structural properties of Aβ (1-40) peptide in protonation stage of one, two, and three: New insights from the histidine protonation behaviors

Structural properties and aggregation tendency can be significantly influenced by histidine behaviors (histidine on N^δ-H state is defined as δ, likewise, N^ε-H: ε and both N^δ-H and N^ε-H: p). In current study, we investigated structural properties of Aβ(1-40) peptide during protonation evolution stage of one, two, and three by total 19 independent replica exchange molecular dynamics simulations using implicit solvent. Our results show that any kind of protonated state will promote β-sheet structure formation in comparison with deprotonated (εεε). With increase in number of protonation, the lowest β-sheet content increased. The highest averaged β-sheet structure content was detected in (δpδ) (46.0 %), (εpp) (36.8 %), and (ppp) (16.0 %) in each protonation stage. With three β-strand structures, (δpδ) shows more stable features and high hydrophobic properties. Further analysis confirmed that H13 and H14 are more important than H6. Specifically, H13 and H14 have a synergistic effect for structural formations by controlling H-bond networks in H13(p) with V39/V40 and H14(p/δ) with G37/G38. Finally, the Pearson correlation coefficient results confirmed that experimental result (ref. 44) is corresponding to our (εpp) system. Our current study will be conducive to understanding the effects of the histidine behaviors, it provides new insights for exploration protein folding and misfolding processes.

Rapamycin-induced autophagy in osteosarcoma cells is mediated via the biglycan/Wnt/β-catenin signaling axis

Biglycan is a class I secreted small leucine-rich proteoglycan (SLRP), which regulates signaling pathways connected to bone pathologies. Autophagy is a vital catabolic process with a dual role in cancer progression. Here, we show that biglycan inhibits autophagy in two osteosarcoma cell lines (P ≤ 0.001), while rapamycin-induced autophagy decreases biglycan expression in MG63 osteosarcoma cells and abrogates the biglycan-induced cell growth increase (P ≤ 0.001). Rapamycin also inhibits β-catenin translocation to the nucleus, inhibiting the Wnt pathway (P ≤ 0.001) and reducing biglycan's colocalization with the Wnt coreceptor LRP6 (P ≤ 0.05). Furthermore, biglycan exhibits protective effects against the chemotherapeutic drug doxorubicin in MG63 OS cells through an autophagy-dependent manner (P ≤ 0.05). Cotreatment of these cells with rapamycin and doxorubicin enhances cells response to doxorubicin by decreasing biglycan (P ≤ 0.001) and β-catenin (P ≤ 0.05) expression. Biglycan deficiency leads to increased caspase-3 activation (P ≤ 0.05), suggesting increased apoptosis of biglycan-deficient cells treated with doxorubicin. Computational models of LRP6 and biglycan complexes suggest that biglycan changes the receptor's ability to interact with other signaling molecules by affecting the interdomain bending angles in the receptor structure. Biglycan binding to LRP6 activates the Wnt pathway and β-catenin nuclear translocation by disrupting β-catenin degradation complex (P ≤ 0.01 and P ≤ 0.05). Interestingly, this mechanism is not followed in moderately differentiated, biglycan-nonexpressing U-2OS OS cells. To sum up, biglycan exhibits protective effects against the doxorubicin in MG63 OS cells by activating the Wnt signaling pathway and inhibiting autophagy.

Metamorphic proteins under a computational microscope: Lessons from a fold-switching RfaH protein

Metamorphic proteins constitute unexpected paradigms of the protein folding problem, as their sequences encode two alternative folds, which reversibly interconvert within biologically relevant timescales to trigger different cellular responses. Once considered a rare aberration, metamorphism may be common among proteins that must respond to rapidly changing environments, exemplified by NusG-like proteins, the only transcription factors present in every domain of life. RfaH, a specialized paralog of bacterial NusG, undergoes an all-α to all-β domain switch to activate expression of virulence and conjugation genes in many animal and plant pathogens and is the quintessential example of a metamorphic protein. The dramatic nature of RfaH structural transformation and the richness of its evolutionary history makes for an excellent model for studying how metamorphic proteins switch folds. Here, we summarize the structural and functional evidence that sparked the discovery of RfaH as a metamorphic protein, the experimental and computational approaches that enabled the description of the molecular mechanism and refolding pathways of its structural interconversion, and the ongoing efforts to find signatures and general properties to ultimately describe the protein metamorphome.

CHARMM-GUI Implicit Solvent Modeler for Various Generalized Born Models in Different Simulation Programs

Implicit solvent models are widely used because they are advantageous to speed up simulations by drastically decreasing the number of solvent degrees of freedom, which allows one to achieve long simulation time scales for large system sizes. CHARMM-GUI, a web-based platform, has been developed to support the setup of complex multicomponent molecular systems and prepare input files. This study describes an Implicit Solvent Modeler (ISM) in CHARMM-GUI for various generalized Born (GB) implicit solvent simulations in different molecular dynamics programs such as AMBER, CHARMM, GENESIS, NAMD, OpenMM, and Tinker. The GB models available in ISM include GB-HCT, GB-OBC, GB-neck, GBMV, and GBSW with the CHARMM and Amber force fields for protein, DNA, RNA, glycan, and ligand systems. Using the system and input files generated by ISM, implicit solvent simulations of protein, DNA, and RNA systems produce similar results for different simulation packages with the same input information. Protein-ligand systems are also considered to further validate the systems and input files generated by ISM. Simple ligand root-mean-square deviation (RMSD) and molecular mechanics generalized Born surface area (MM/GBSA) calculations show that the performance of implicit simulations is better than docking and can be used for early stage ligand screening. These reasonable results indicate that ISM is a useful and reliable tool to provide various implicit solvent simulation applications.

Fundamental Redesign of the TIGER2hs Kernel to Address Severe Parameter Sensitivity

Replica exchange molecular dynamics simulations are one of the most popular approaches to enhance conformational sampling of molecular systems. Applications range from protein folding to protein-protein or other host-guest interactions, as well as binding free energy calculations. While these methods are computationally expensive, highly accurate results can be obtained. We recently developed TIGER2hs, an improved version of the temperature intervals with global exchange of replicas (TIGER2) algorithm. This method combines the replica-based enhanced sampling in an explicit solvent with a hybrid solvent energy evaluation. During the exchange attempts, bulk water is replaced by an implicit solvent model, allowing sampling with significantly less replicas than parallel tempering (REMD). This enables accurate enhanced sampling calculations with only a fraction of computational resources compared to REMD. Our latest results highlight several issues with sampling imbalance and parameter sensitivity within the original TIGER2 exchange algorithms that affect the overall state populations. A high sensitivity on replica number and maximum temperature is eliminated by changing to a pairwise exchange kernel (PE) without additional sorting. Simulations are controlled by adjusting the average temperature change per exchange ⟨ΔT/χ⟩ to below 30 K to mimic a controlled temperature mixing of replicas similar to REMD. Thus, this parameter provides an applicable property for selecting combinations of replica number and maximum temperature to adjust simulations for best accuracy, with flexible resource investment. This increases the robustness of the method and ensures results in excellent agreement with REMD, as demonstrated for three different peptides.

Mechanical Stability of Ribonuclease A Heavily Depends on the Redox Environment

Disulfide bonds are covalent bonds that connect nonlocal fragments of proteins, and they are unique post-translational modifications of proteins. They require the oxidizing environment to be stable, which occurs for example during oxidative stress; however, in a cell the reductive environment is maintained, lowering their stability. Despite many years of research on disulfide bonds, their role in the protein life cycle is not fully understood and seems to strictly depend on a system or process in which they are involved. In this article, coarse-grained UNited RESidue (UNRES), and all-atom Assisted Model Building with Energy Refinement (AMBER) force fields were applied to run a series of steered molecular dynamics (SMD) simulations of one of the most studied, but still not fully understood, proteins─ribonuclease A (RNase A). SMD simulations were performed to study the mechanical stability of RNase A in different oxidative-reductive environments. As disulfide bonds (and any other covalent bonds) cannot break/form in any classical all-atom force field, we applied additional restraints between sulfur atoms of reduced cysteines which were able to mimic the breaking of the disulfide bonds. On the other hand, the coarse-grained UNRES force field enables us to study the breaking/formation of the disulfide bonds and control the reducing/oxidizing environment owing to the presence of the designed distance/orientation-dependent potential. This study reveals that disulfide bonds have a strong influence on the mechanical stability of RNase A only in a highly oxidative environment. However, the local stability of the secondary structure seems to play a major factor in the overall stability of the protein. Both our thermal unfolding and mechanical stretching studies show that the most stable disulfide bond is Cys65-Cys72. The breaking of disulfide bonds Cys26-Cys84 and Cys58-Cys110 is associated with large force peaks. They are structural bridges, which are mostly responsible for stabilizing the RNase A conformation, while the presence of the remaining two bonds (Cys65-Cys72 and Cys40-Cys95) is most likely connected with the enzymatic activity rather than the structural stability of RNase A in the cytoplasm. Our results prove that disulfide bonds are indeed stabilizing fragments of the proteins, but their role is strongly redox environment-dependent.

Deciphering the Folding Mechanism of Proteins G and L and Their Mutants

Much of our understanding of folding mechanisms comes from interpretations of experimental ϕ and ψ value analysis, relating the differences in stability of the transition state ensemble (TSE) and folded state. We introduce a unified approach combining simulations and Bayesian inference to provide atomistic detail for the folding mechanism of proteins G and L and their mutants. Proteins G and L fold to similar topologies despite low sequence similarity, but differ in their folding pathways. A fast folding redesign of protein G, NuG2, switches folding pathways and folds through a similar pathway with protein L. A redesign of protein L also leads to faster folding, respecting the original folding pathway. Our Bayesian inference approach starts from the same prior on all systems and correctly identifies the folding mechanism for each of the four proteins, a success of the force field and sampling strategy. The approach is computationally efficient and correctly identifies the TSE and intermediate structures along the folding pathway in good agreement with experiments. We complement our findings by using two orthogonal approaches that differ in computational cost and interpretability. Adaptive sampling MD combined with the Markov state model provides a kinetic model that confirms the more complex folding mechanism of protein G and its mutant. Finally, a novel fragment decomposition approach using AlphaFold identifies preferences for secondary structure element combinations that follow the order of events observed in the folding pathways.

Explicit solvent repulsive scaling replica exchange molecular dynamics (RS-REMD) in molecular modeling of protein-glycosaminoglycan complexes

Glycosaminoglcyans (GAGs), linear anionic periodic polysaccharides, are crucial for many biologically relevant functions in the extracellular matrix. By interacting with proteins GAGs mediate processes such as cancer development, cell proliferation and the onset of neurodegenerative diseases. Despite this eminent importance of GAGs, they still represent a limited focus for the computational community in comparison to other classes of biomolecules. Therefore, there is a lack of modeling tools designed specifically for docking GAGs. One has to rely on existing docking software developed mostly for small drug molecules substantially differing from GAGs in their basic physico-chemical properties. In this study, we present an updated protocol for docking GAGs based on the Repulsive Scaling Replica Exchange Molecular Dynamics (RS-REMD) that includes explicit solvent description. The use of this water model improved docking performance both in terms of its accuracy and speed. This method represents a significant computational progress in GAG-related research.

Accelerating the Ensemble Convergence of RNA Hairpin Simulations with a Replica Exchange Structure Reservoir

RNA is a key participant in many biological processes, but studies of RNA using computer simulations lag behind those of proteins, largely due to less-developed force fields and the slow dynamics of RNA. Generating converged RNA ensembles for force field development and other studies remains a challenge. In this study, we explore the ability of replica exchange molecular dynamics to obtain well-converged conformational ensembles for two RNA hairpin systems in an implicit solvent. Even for these small model systems, standard REMD remains computationally costly, but coupling to a pre-generated structure library using the reservoir REMD approach provides a dramatic acceleration of ensemble convergence for both model systems. Such precise ensembles could facilitate RNA force field development and validation and applications of simulation to more complex RNA systems. The advantages and remaining challenges of applying R-REMD to RNA are investigated in detail.

Generalized Born Implicit Solvent Models Do Not Reproduce Secondary Structures of De Novo Designed Glu/Lys Peptides

We test a range of standard generalized Born (GB) models and protein force fields for a set of five experimentally characterized, designed peptides comprising alternating blocks of glutamate and lysine, which have been shown to differ significantly in α-helical content. Sixty-five combinations of force fields and GB models are evaluated in >800 μs of molecular dynamics simulations. GB models generally do not reproduce the experimentally observed α-helical content, and none perform well for all five peptides. These results illustrate that these models are not usefully predictive in this context. These peptides provide a useful test set for simulation methods.

In Silico Modeling and Scoring of PROTAC-Mediated Ternary Complex Poses

Proteolysis targeting chimeras (PROTACs) are molecules that induce protein degradation via formation of ternary complexes between an E3 ubiquitin ligase and a target protein. The rational design of PROTACs requires accurate knowledge of the native configuration of the PROTAC-induced ternary complex. This study demonstrates that native and non-native ternary complex poses can be distinguished based on the pose occupancy time in MD, where native poses exhibit longer occupancy times at both room and higher temperatures. Candidate poses are generated by MD sampling and pre-ranked by classic MM/GBSA. A specific heating scheme is then applied to accelerate ternary pose departure, with the pose occupancy time and fraction being measured. This scoring identifies the native pose in all systems tested. Its success is partially attributed to the dynamic nature of pose departure analyses, which accounts for entropic effects typically neglected in the faster static scoring methods, while entropy plays a greater role in protein-protein than in protein-ligand systems.

Accelerating Protein Folding Molecular Dynamics Using Inter-Residue Distances from Machine Learning Servers

Recently, predicting the native structures of proteins has become possible using computational molecular physics (CMP)─physics-based force fields sampled with proper statistics─but only for small proteins. Algorithms with better scaling are needed. We describe ML x MELD x MD, a molecular dynamics (MD) method that inputs residue contacts derived from machine learning (ML) servers into MELD, a Bayesian accelerator that preserves detailed-balance statistics. Contacts are derived from trRosetta-predicted distance histograms (distograms) and are integrated into MELD's atomistic MD as spatial restraints through parametrized potential functions. In the CASP14 blind prediction event, ML x MELD x MD predicted 13 native structures to better than 4.5 Å error, including for 10 proteins in the range of 115-250 amino acids long. Also, the scaling of simulation time vs protein length is much better than unguided MD: tsim ∼ e0.023N for ML x MELD x MD vs tsim ∼ e0.168N for MD alone. This shows how machine learning information can be leveraged to advance physics-based modeling of proteins.

Biglycan Interacts with Type I Insulin-like Receptor (IGF-IR) Signaling Pathway to Regulate Osteosarcoma Cell Growth and Response to Chemotherapy

Simple Summary

Osteosarcoma (OS) is an aggressive, primary bone cancer. OS cells produce altered osteoid whose components participate in signaling correlated to the development of this cancer. Biglycan (BGN), a proteoglycan, is correlated to aggressive OS type and resistance to chemotherapy. A constitutive signaling of insulin-like growth factor receptor I (IGF-IR) signaling in sarcoma progression was established. We showed that biglycan binds IGF-IR resulting in prolonged IGF-IR activation, nuclear translocation, and growth response of the poorly-differentiated MG63 cells correlated to increased aggressiveness markers expression and enhanced chemoresistance. This mechanism is not valid in moderately and well-differentiated, biglycan non-expressing U-2OS and Saos-2 OS cells.

Abstract

Osteosarcoma (OS) is a mesenchymally derived, aggressive bone cancer. OS cells produce an aberrant nonmineralized or partly mineralized extracellular matrix (ECM) whose components participate in signaling pathways connected to specific pathogenic phenotypes of this bone cancer. The expression of biglycan (BGN), a secreted small leucine-rich proteoglycan (SLRP), is correlated to aggressive OS phenotype and resistance to chemotherapy. A constitutive signaling of IGF-IR signaling input in sarcoma progression has been established. Here, we show that biglycan activates the IGF-IR signaling pathway to promote MG63 biglycan-secreting OS cell growth by forming a complex with the receptor. Computational models of IGF-IR and biglycan docking suggest that biglycan binds IGF-IR dimer via its concave surface. Our binding free energy calculations indicate the formation of a stable complex. Biglycan binding results in prolonged IGF-IR activation leading to protracted IGF-IR-dependent cell growth response of the poorly-differentiated MG63 cells. Moreover, biglycan facilitates the internalization (p ≤ 0.01, p ≤ 0.001) and sumoylation-enhanced nuclear translocation of IGF-IR (p ≤ 0.05) and its DNA binding in MG63 cells (p ≤ 0.001). The tyrosine kinase activity of the receptor mediates this mechanism. Furthermore, biglycan downregulates the expression of the tumor-suppressor gene, PTEN (p ≤ 0.01), and increases the expression of endothelial–mesenchymal transition (EMT) and aggressiveness markers vimentin (p ≤ 0.01) and fibronectin (p ≤ 0.01) in MG63 cells. Interestingly, this mechanism is not valid in moderately and well-differentiated, biglycan non-expressing U-2OS and Saos-2 OS cells. Furthermore, biglycan exhibits protective effects against the chemotherapeutic drug, doxorubicin, in MG63 OS cells (p ≤ 0.01). In conclusion, these data indicate a potential direct and adjunct therapeutical role of biglycan in osteosarcoma.

Searching for Low Probability Opening Events in a DNA Sliding Clamp

The β subunit of E. coli DNA polymererase III is a DNA sliding clamp associated with increasing the processivity of DNA synthesis. In its free form, it is a circular homodimer structure that can accomodate double-stranded DNA in a nonspecific manner. An open state of the clamp must be accessible before loading the DNA. The opening mechanism is still a matter of debate, as is the effect of bound DNA on opening/closing kinetics. We use a combination of atomistic, coarse-grained, and enhanced sampling strategies in both explicit and implicit solvents to identify opening events in the sliding clamp. Such simulations of large nucleic acid and their complexes are becoming available and are being driven by improvements in force fields and the creation of faster computers. Different models support alternative opening mechanisms, either through an in-plane or out-of-plane opening event. We further note some of the current limitations, despite advances, in modeling these highly charged systems with implicit solvent.