Quality Threshold Clustering of Molecular Dynamics: A Word of Caution

NF2 regulates IP3R-mediated Ca2+ signal and apoptosis in meningiomas

Meningiomas are the most common primary intracranial tumors and account for nearly 30% of all nervous system tumors. Approximately half of meningioma patients exhibit neurofibromin 2 (NF2) gene inactivation. Here, NF2 was shown to interact with the endoplasmic reticulum (ER) calcium (Ca2+) channel inositol 1,4,5-trisphosphate receptor 1 (IP3R1) in IOMM-Lee, a high-grade malignant meningioma cell line, and the F1 subdomain of NF2 plays a critical role in this interaction. Functional assays indicated that NF2 promotes the phosphorylation of IP3R (Ser 1756) and IP3R-mediated endoplasmic reticulum (ER) Ca2+ release by binding to IP3R1, which results in Ca2+-dependent apoptosis. Knockout of NF2 decreased Ca2+ release and promoted resistance to apoptosis, which was rescued by wild-type NF2 overexpression but not by F1 subdomain deletion truncation overexpression. The effects of NF2 defects on the development of tumors were further studied in mouse models. The decreased expression level of NF2 caused by NF2 gene knockout or mutation affects the activity of the IP3R channel, which reduces Ca2+-dependent apoptosis, thereby promoting the development of tumors. We elucidated the interaction patterns of NF2 and IP3R1, revealed the molecular mechanism through which NF2 regulates IP3R1-mediated Ca2+ release, and elucidated the new pathogenic mechanism of meningioma-related NF2 variants. Our study broadens the current understanding of the biological function of NF2 and provides ideas for drug screening of NF2-associated meningioma.

k-Means NANI: An Improved Clustering Algorithm for Molecular Dynamics Simulations

One of the key challenges of k-means clustering is the seed selection or the initial centroid estimation since the clustering result depends heavily on this choice. Alternatives such as k-means++ have mitigated this limitation by estimating the centroids using an empirical probability distribution. However, with high-dimensional and complex data sets such as those obtained from molecular simulation, k-means++ fails to partition the data in an optimal manner. Furthermore, stochastic elements in all flavors of k-means++ will lead to a lack of reproducibility. K-means N-Ary Natural Initiation (NANI) is presented as an alternative to tackle this challenge by using efficient n-ary comparisons to both identify high-density regions in the data and select a diverse set of initial conformations. Centroids generated from NANI are not only representative of the data and different from one another, helping k-means to partition the data accurately, but also deterministic, providing consistent cluster populations across replicates. From peptide and protein folding molecular simulations, NANI was able to create compact and well-separated clusters as well as accurately find the metastable states that agree with the literature. NANI can cluster diverse data sets and be used as a standalone tool or as part of our MDANCE clustering package.

Identification of Small Molecule Ligand Binding Sites On and In the ARNT PAS-B Domain

Transcription factors are generally challenging to target with small molecule inhibitors due to their structural plasticity and lack of catalytic sites. Notable exceptions include several naturally ligand-regulated transcription factors, including our prior work with the heterodimeric HIF-2 transcription factor which showed that small molecule binding within an internal pocket of the HIF-2α PAS-B domain can disrupt its interactions with its dimerization partner, ARNT. Here, we explore the feasibility of similarly targeting small molecules to the analogous ARNT PAS-B domain itself, potentially opening a promising route to simultaneously modulate several ARNT-mediated signaling pathways. Using solution NMR screening of an in-house fragment library, we previously identified several compounds that bind ARNT PAS-B and, in certain cases, antagonize ARNT association with the TACC3 transcriptional coactivator. However, these ligands have only modest binding affinities, complicating characterization of their binding sites. We address this challenge by combining NMR, MD simulations, and ensemble docking to identify ligand-binding 'hotspots' on and within the ARNT PAS-B domain. Our data indicate that the two ARNT/TACC3 inhibitors, KG-548 and KG-655, bind to a β-sheet surface implicated in both HIF-2 dimerization and coactivator recruitment. Furthermore, while KG-548 binds exclusively to the β-sheet surface, KG-655 can additionally bind within a water-accessible internal cavity in ARNT PAS-B. Finally, KG-279, while not a coactivator inhibitor, exemplifies ligands that preferentially bind only to the internal cavity. All three ligands promoted ARNT PAS-B homodimerization, albeit to varying degrees. Taken together, our findings provide a comprehensive overview of ARNT PAS-B ligand-binding sites and may guide the development of more potent coactivator inhibitors for cellular and functional studies.

Dynamics of Protein-RNA Interfaces Using All-Atom Molecular Dynamics Simulations

Facing the current challenges posed by human health diseases requires the understanding of cell machinery at a molecular level. The interplay between proteins and RNA is key for any physiological phenomenon, as well protein-RNA interactions. To understand these interactions, many experimental techniques have been developed, spanning a very wide range of spatial and temporal resolutions. In particular, the knowledge of tridimensional structures of protein-RNA complexes provides structural, mechanical, and dynamical pieces of information essential to understand their functions. To get insights into the dynamics of protein-RNA complexes, we carried out all-atom molecular dynamics simulations in explicit solvent on nine different protein-RNA complexes with different functions and interface size by taking into account the bound and unbound forms. First, we characterized structural changes upon binding and, for the RNA part, the change in the puckering. Second, we extensively analyzed the interfaces, their dynamics and structural properties, and the structural waters involved in the binding, as well as the contacts mediated by them. Based on our analysis, the interfaces rearranged during the simulation time showing alternative and stable residue-residue contacts with respect to the experimental structure.

Structural mechanism of voltage-gated sodium channel slow inactivation

Voltage-gated sodium (NaV) channels mediate a plethora of electrical activities. NaV channels govern cellular excitability in response to depolarizing stimuli. Inactivation is an intrinsic property of NaV channels that regulates cellular excitability by controlling the channel availability. The fast inactivation, mediated by the Ile-Phe-Met (IFM) motif and the N-terminal helix (N-helix), has been well-characterized. However, the molecular mechanism underlying NaV channel slow inactivation remains elusive. Here, we demonstrate that the removal of the N-helix of NaVEh (NaVEhΔN) results in a slow-inactivated channel, and present cryo-EM structure of NaVEhΔN in a potential slow-inactivated state. The structure features a closed activation gate and a dilated selectivity filter (SF), indicating that the upper SF and the inner gate could serve as a gate for slow inactivation. In comparison to the NaVEh structure, NaVEhΔN undergoes marked conformational shifts on the intracellular side. Together, our results provide important mechanistic insights into NaV channel slow inactivation.

Small angle X-ray scattering analysis of thermophilic cytochrome P450 CYP119 and the effects of the N-terminal histidine tag

Combining size exclusion chromatography-small angle X-ray scattering (SEC-SAXS) and molecular dynamics (MD) analysis is a promising approach to investigate protein behavior in solution, particularly for understanding conformational changes due to substrate binding in cytochrome P450s (CYPs). This study investigates conformational changes in CYP119, a thermophilic CYP from Sulfolobus acidocaldarius that exhibits structural flexibility similar to mammalian CYPs. Although the crystal structure of ligand-free (open state) and ligand-bound (closed state) forms of CYP119 is known, the overall structure of the enzyme in solution has not been explored until now. It was found that theoretical scattering profiles from the crystal structures of CYP119 did not align with the SAXS data, but conformers from MD simulations, particularly starting from the open state (46 % of all frames), agreed well. Interestingly, a small percentage of closed-state conformers also fit the data (9 %), suggesting ligand-free CYP119 samples ligand-bound conformations. Ab initio SAXS models for N-His tagged CYP119 revealed a tail-like unfolded structure impacting protein flexibility, which was confirmed by in silico modeling. SEC-SAXS analysis of N-His CYP119 indicated pentameric structures in addition to monomers in solution, affecting the stability and activity of the enzyme. This study adds insights into the conformational dynamics of CYP119 in solution.

k-Means NANI: an improved clustering algorithm for Molecular Dynamics simulations

One of the key challenges of k-means clustering is the seed selection or the initial centroid estimation since the clustering result depends heavily on this choice. Alternatives such as k-means++ have mitigated this limitation by estimating the centroids using an empirical probability distribution. However, with high-dimensional and complex datasets such as those obtained from molecular simulation, k-means++ fails to partition the data in an optimal manner. Furthermore, stochastic elements in all flavors of k-means++ will lead to a lack of reproducibility. K-means N-Ary Natural Initiation (NANI) is presented as an alternative to tackle this challenge by using efficient n-ary comparisons to both identify high-density regions in the data and select a diverse set of initial conformations. Centroids generated from NANI are not only representative of the data and different from one another, helping k-means to partition the data accurately, but also deterministic, providing consistent cluster populations across replicates. From peptide and protein folding molecular simulations, NANI was able to create compact and well-separated clusters as well as accurately find the metastable states that agree with the literature. NANI can cluster diverse datasets and be used as a standalone tool or as part of our MDANCE clustering package.

Dynamics and Conformations of a Full-Length CRESS-DNA Replicase

Circular Rep-encoding single-stranded DNA (CRESS-DNA) viruses encode for a Replicase (Rep) that is essential for viral replication. Rep is a helicase with three domains: an endonuclease, an oligomeric, and an ATPase domain (ED, OD, and AD). Our recent cryo-EM structure of the porcine circovirus 2 (PCV2) Rep provided the first structure of a CRESS-DNA Rep. The structure visualized the ED to be highly mobile, Rep to form a homo-hexamer, bound ssDNA and nucleotides, and the AD to adopt a staircase arrangement around the ssDNA. We proposed a hand-over-hand mechanism by the ADs for ssDNA translocation. The hand-over-hand mechanism requires extensive movement of the AD. Here, we scrutinize this mechanism using all-atom Molecular Dynamics (MD) simulation of Rep in three states: (1) Rep bound to ssDNA and ADP, (2) Rep bound to ssDNA, and (3) Rep by itself. Each of the 700 nsec simulations converges within 200 nsec and provides important insight into the dynamics of Rep, the dynamics of Rep in the presence of these biomolecules, and the importance of ssDNA and ADP in driving the AD to adopt the staircase arrangement around the ssDNA. To the best of our knowledge, this is the first example of an all-atom MD simulation of a CRESS-DNA Rep. This study sets the basis of further MD studies aimed at obtaining a chemical understanding of how Rep uses nucleotide binding and hydrolysis to translocate ssDNA.

Structural basis of human Slo2.2 channel gating and modulation

The sodium-activated Slo2.2 channel is abundantly expressed in the brain, playing a critical role in regulating neuronal excitability. The Na+-binding site and the underlying mechanisms of Na+-dependent activation remain unclear. Here, we present cryoelectron microscopy (cryo-EM) structures of human Slo2.2 in closed, open, and inhibitor-bound form at resolutions of 2.6-3.2 Å, revealing gating mechanisms of Slo2.2 regulation by cations and a potent inhibitor. The cytoplasmic gating ring domain of the closed Slo2.2 harbors multiple K+ and Zn2+ sites, which stabilize the channel in the closed conformation. The open Slo2.2 structure reveals at least two Na+-sensitive sites where Na+ binding induces expansion and rotation of the gating ring that opens the inner gate. Furthermore, a potent inhibitor wedges into a pocket formed by pore helix and S6 helix and blocks the pore. Together, our results provide a comprehensive structural framework for the investigation of Slo2.2 channel gating, Na+ sensation, and inhibition.

Embedding Nonrigid Solutes in an Averaged Environment: A Case Study on Rhodopsins

Many simulation methods concerning solvated molecules are based on the assumption that the solvated species and the solvent can be characterized by some representative structures of the solute and some embedding potential corresponding to this structure. While the averaging of the solvent configurations to obtain an embedding potential has been studied in great detail, this hinges on a single solute structure representation. This assumption is re-examined and generalized for conformationally flexible solutes and tested on 4 nonrigid systems. In this generalized approach, the solute is characterized by a set of representative structures and the corresponding embedding potentials. The representative structures are identified by means of subdividing the statistical ensemble, which in this work is generated by a constant-temperature molecular dynamics simulation. The embedding potential defined in the Frozen-Density Embedding Theory is used to characterize the average effect of the solvent in each subensemble. The numerical examples concern the vertical excitation energies of protonated retinal Schiff bases in protein environments. It is comprehensively shown that subensemble averaging leads to huge computational savings compared with explicit averaging of the excitation energies in the whole ensemble while introducing only minor errors in the case of the systems examined.

Integrated NMR/Molecular Dynamics Determination of the Ensemble Conformation of a Thermodynamically Stable CUUG RNA Tetraloop

Both experimental and theoretical structure determinations of RNAs have remained challenging due to the intrinsic dynamics of RNAs. We report here an integrated nuclear magnetic resonance/molecular dynamics (NMR/MD) structure determination approach to describe the dynamic structure of the CUUG tetraloop. We show that the tetraloop undergoes substantial dynamics, leading to averaging of the experimental data. These dynamics are particularly linked to the temperature-dependent presence of a hydrogen bond within the tetraloop. Interpreting the NMR data by a single structure represents the low-temperature structure well but fails to capture all conformational states occurring at a higher temperature. We integrate MD simulations, starting from structures of CUUG tetraloops within the Protein Data Bank, with an extensive set of NMR data, and provide a structural ensemble that describes the dynamic nature of the tetraloop and the experimental NMR data well. We thus show that one of the most stable and frequently found RNA tetraloops displays substantial dynamics, warranting such an integrated structural approach.

Clustering molecular dynamics conformations of the CC'-loop of the PD-1 immuno-checkpoint receptor

Molecular mechanisms within the checkpoint receptor PD-1 are essential for its activation by PD-L1 as well as for blocking such an activation via checkpoint inhibitors. We use molecular dynamics to scrutinize patterns of atomic motion in PD-1 without a ligand. Molecular dynamics is performed for the whole extracellular domain of PD-1, and the analysis focuses on its CC'-loop and some adjacent Cα-atoms. We extend previous work by applying common nearest neighbor clustering (Cnn) and compare the performance of this method with Daura clustering as well as UMAP dimension reduction and subsequent agglomerative linkage clustering. As compared to Daura clustering, we found Cnn less sensitive to cutoff selection and better able to return representative clusters for sets of different 3D atomic conformations. Interestingly, Cnn yields results quite similar to UMAP plus linkage clustering.

Structural Perturbation of Monomers Determines the Amyloid Aggregation Propensity of Calcitonin Variants

Human calcitonin (hCT) is a polypeptide hormone that participates in calcium-phosphorus metabolism. Irreversible aggregation of 32-amino acid hCT into β-sheet-rich amyloid fibrils impairs physiological activity and increases the risk of medullary carcinoma of the thyroid. Amyloid-resistant hCT derivatives substituting critical amyloidogenic residues are of particular interest for clinical applications as therapeutic drugs against bone-related diseases. Uncovering the aggregation mechanism of hCT at the molecular level, therefore, is important for the design of amyloid-resistant hCT analogues. Here, we investigated the aggregation dynamics of hCT, non-amyloidogenic salmon calcitonin (sCT), and two hCT analogues with reduced aggregation tendency─TL-hCT and phCT─using long timescale discrete molecular dynamics simulations. Our results showed that hCT monomers mainly adopted unstructured conformations with dynamically formed helices around the central region. hCT self-assembled into helix-rich oligomers first, followed by a conformational conversion into β-sheet-rich oligomers with β-sheets formed by residues 10-30 and stabilized by aromatic and hydrophobic interactions. Our simulations confirmed that TL-hCT and phCT oligomers featured more helices and fewer β-sheets than hCT. Substitution of central aromatic residues with leucine in TL-hCT and replacing C-terminal hydrophobic residue with hydrophilic amino acid in phCT only locally suppressed β-sheet propensities in the central region and C-terminus, respectively. Having mutations in both central and C-terminal regions, sCT monomers and dynamically formed oligomers predominantly adopted helices, confirming that both central aromatic and C-terminal hydrophobic residues played important roles in the fibrillization of hCT. We also observed the formation of β-barrel intermediates, postulated as the toxic oligomers in amyloidosis, for hCT but not for sCT. Our computational study depicts a complete picture of the aggregation dynamics of hCT and the effects of mutations. The design of next-generation amyloid-resistant hCT analogues should consider the impact on both amyloidogenic regions and also take into account the amplification of transient β-sheet population in monomers upon aggregation.

Examining unsupervised ensemble learning using spectroscopy data of organic compounds

One solution to the challenge of choosing an appropriate clustering algorithm is to combine different clusterings into a single consensus clustering result, known as cluster ensemble (CE). This ensemble learning strategy can provide more robust and stable solutions across different domains and datasets. Unfortunately, not all clusterings in the ensemble contribute to the final data partition. Cluster ensemble selection (CES) aims at selecting a subset from a large library of clustering solutions to form a smaller cluster ensemble that performs as well as or better than the set of all available clustering solutions. In this paper, we investigate four CES methods for the categorization of structurally distinct organic compounds using high-dimensional IR and Raman spectroscopy data. Single quality selection (SQI) forms a subset of the ensemble by selecting the highest quality ensemble members. The Single Quality Selection (SQI) method is used with various quality indices to select subsets by including the highest quality ensemble members. The Bagging method, usually applied in supervised learning, ranks ensemble members by calculating the normalized mutual information (NMI) between ensemble members and consensus solutions generated from a randomly sampled subset of the full ensemble. The hierarchical cluster and select method (HCAS-SQI) uses the diversity matrix of ensemble members to select a diverse set of ensemble members with the highest quality. Furthermore, a combining strategy can be used to combine subsets selected using multiple quality indices (HCAS-MQI) for the refinement of clustering solutions in the ensemble. The IR + Raman hybrid ensemble library is created by merging two complementary "views" of the organic compounds. This inherently more diverse library gives the best full ensemble consensus results. Overall, the Bagging method is recommended because it provides the most robust results that are better than or comparable to the full ensemble consensus solutions.

Substitution of PINK1 Gly411 modulates substrate receptivity and turnover

The ubiquitin (Ub) kinase-ligase pair PINK1-PRKN mediates the degradation of damaged mitochondria by macroautophagy/autophagy (mitophagy). PINK1 surveils mitochondria and upon stress accumulates on the mitochondrial surface where it phosphorylates serine 65 of Ub to activate PRKN and to drive mitochondrial turnover. While loss of either PINK1 or PRKN is genetically linked to Parkinson disease (PD) and activating the pathway seems to have great therapeutic potential, there is no formal proof that stimulation of mitophagy is always beneficial. Here we used biochemical and cell biological methods to study single nucleotide variants in the activation loop of PINK1 to modulate the enzymatic function of this kinase. Structural modeling and in vitro kinase assays were used to investigate the molecular mechanism of the PINK1 variants. In contrast to the PD-linked PINK1^G411S mutation that diminishes Ub kinase activity, we found that the PINK1^G411A variant significantly boosted Ub phosphorylation beyond levels of PINK1 wild type. This resulted in augmented PRKN activation, mitophagy rates and increased viability after mitochondrial stress in midbrain-derived, gene-edited neurons. Mechanistically, the G411A variant stabilizes the kinase fold of PINK1 and transforms Ub to adopt the preferred, C-terminally retracted conformation for improved substrate turnover. In summary, we identify a critical role of residue 411 for substrate receptivity that may now be exploited for drug discovery to increase the enzymatic function of PINK1. The genetic substitution of Gly411 to Ala increases mitophagy and may be useful to confirm neuroprotection in vivo and might serve as a critical positive control during therapeutic development.Abbreviations: ATP: adenosine triphosphate; CCCP: carbonyl cyanide m-chlorophenyl hydrazone; Ub-CR: ubiquitin with C-terminally retracted tail; CTD: C-terminal domain (of PINK1); ELISA: enzyme-linked immunosorbent assay; HCI: high-content imaging; IB: immunoblot; IF: immunofluorescence; NPC: neuronal precursor cells; MDS: molecular dynamics simulation; PD: Parkinson disease; p-S65-Ub: ubiquitin phosphorylated at Ser65; RMSF: root mean scare fluctuation; TOMM: translocase of outer mitochondrial membrane; TVLN: ubiquitin with T66V and L67N mutation, mimics Ub-CR; Ub: ubiquitin; WT: wild-type.

MDSCAN: RMSD-based HDBSCAN clustering of long molecular dynamics

MOTIVATION

The term clustering designates a comprehensive family of unsupervised learning methods allowing to group similar elements into sets called clusters. Geometrical clustering of molecular dynamics (MD) trajectories is a well-established analysis to gain insights into the conformational behavior of simulated systems. However, popular variants collapse when processing relatively long trajectories because of their quadratic memory or time complexity. From the arsenal of clustering algorithms, HDBSCAN stands out as a hierarchical density-based alternative that provides robust differentiation of intimately related elements from noise data. Although a very efficient implementation of this algorithm is available for programming-skilled users (HDBSCAN*), it cannot treat long trajectories under the de facto molecular similarity metric RMSD.

RESULTS

Here, we propose MDSCAN, an HDBSCAN-inspired software specifically conceived for non-programmers users to perform memory-efficient RMSD-based clustering of long MD trajectories. Methodological improvements over the original version include the encoding of trajectories as a particular class of vantage-point tree (decreasing time complexity), and a dual-heap approach to construct a quasi-minimum spanning tree (reducing memory complexity). MDSCAN was able to process a trajectory of 1 million frames using the RMSD metric in about 21 h with <8 GB of RAM, a task that would have taken a similar time but more than 32 TB of RAM with the accelerated HDBSCAN* implementation generally used.

AVAILABILITY AND IMPLEMENTATION

The source code and documentation of MDSCAN are free and publicly available on GitHub (https://github.com/LQCT/MDScan.git) and as a PyPI package (https://pypi.org/project/mdscan/).

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

On Quality Thresholds for the Clustering of Molecular Structures

It has been recently suggested that diametral (so-called quality) similarity thresholds are superior to radial ones for the clustering of molecular three-dimensional structures (González-Alemán et al., 2020). The argument has been made for two clustering algorithms available in various software packages for the analysis of molecular structures from ensembles generated by computer simulations, attributed to Daura et al. (1999) (radial threshold) and Heyer et al. (1999) (diametral threshold). Here, we compare these two algorithms using the root-mean-squared difference (rmsd) between the Cartesian coordinates of selected atoms as pairwise similarity metric. We discuss formally the relation between these two methods and illustrate their behavior with two examples, a set of points in two dimensions and the coordinates of the tau polypeptide along a trajectory extracted from a replica-exchange molecular-dynamics simulation (Shea and Levine, 2016). We show that the two methods produce equally sized clusters as long as adequate choices are made for the respective thresholds. The real issue is not whether the threshold is radial or diametral but how to choose in either case a threshold value that is physically meaningful. We will argue that, when clustering molecular structures with the rmsd as a metric, the simplest best guess for a threshold is actually radial in nature.

Identifying RO9021 as a Potential Inhibitor of PknG from Mycobacterium tuberculosis: Combinative Computational and In Vitro Studies

Tuberculosis (TB) is an infectious disease caused by Mycobacterium tuberculosis (Mtb). Despite being considered curable and preventable, the increase of antibiotic resistance is becoming a serious public health problem. Mtb is a pathogen capable of surviving in macrophages, causing long-term latent infection where the mycobacterial serine/threonine protein kinase G (PknG) plays a protective role. Therefore, PknG is an important inhibitory target to prevent Mtb from entering the latency stage. In this study, we use a pharmacophore-based virtual screening and biochemical assays to identify the compound RO9021 (CHEMBL3237561) as a PknG inhibitor. In detail, 1.5 million molecules were screened using a scalable cloud-based setup, identifying 689 candidates, which were further subjected to additional screening employing molecular docking. Molecular docking spotted 62 compounds with estimated binding affinities of -7.54 kcal/mol (s.d. = 0.77 kcal/mol). Finally, 14 compounds were selected for in vitro experiments considering previously reported biological activities and commercial availability. In vitro assays of PknG activity showed that RO9021 inhibits the kinase activity similarly to AX20017, a known inhibitor. The inhibitory effect was found to be dose dependent with a relative IC₅₀ value of 4.4 ± 1.1 μM. Molecular dynamics simulations predicted that the PknG-RO9021 complex is stable along the tested timescale. Altogether, our study indicates that RO9021 is a noteworthy drug candidate for further developing new anti-TB drugs that hold excellent reported pharmacokinetic parameters.

Designing of disruptor molecules to restrain the protein-protein interaction network of VANG1/SCRIB/NOS1AP using fragment-based drug discovery techniques

Governing protein-protein interaction networks are the cynosure of cell signaling and oncogenic networks. Multifarious processes when aligned with one another can result in a dysregulated output which can result in cancer progression. In the current research, one such network of proteins comprising VANG1/SCRIB/NOS1AP, which is responsible for cell migration, is targeted. The proteins are modeled using in-silico approaches, and the interaction is visualized utilizing protein-protein docking. Designing drugs for the convoluted protein network can serve as a challenging task that can be overcome by fragment-based drug designing, a recent game-changer in the computational drug discovery strategy for protein interaction networks. The model is exposed to the extraction of hotspots, also known as the restrained regions for small molecular hits. The hotspot regions are subjected to a library of generated fragments, which are then recombined and rejoined to develop small molecular disruptors of the macromolecular assemblage. Rapid screening methods using pharmacokinetic tools and 2D interaction studies resulted in four molecules that could serve the purpose of a disruptor. The final validation is executed by long-range simulations of 100 ns and exploring the stability of the complex using several parameters leading to the emergence of two novel molecules VNS003 and VNS005 that could be used as the disruptors of the protein assembly VANG1/SCRIB/NOS1AP. Also, the molecules were explored as single protein targets approbated via molecular docking and 100 ns molecular dynamics simulation. This concluded VNS003 as the most suitable inhibitor module capable of acting as a disruptor of a macromolecular assembly as well as acting on individual protein chains, thus leading to the primary hindrance in the formation of the protein interaction complex.

Successes and challenges in multiscale modelling of artificial metalloenzymes: the case study of POP-Rh2 cyclopropanase

Molecular modelling applications in metalloenzyme design are still scarce due to a series of challenges. On top of that, the simulations of metal-mediated binding and the identification of catalytic competent geometries require both large conformational exploration and simulation of fine electronic properties. Here, we demonstrate how the incorporation of new tools in multiscale strategies, namely substrate diffusion exploration, allows taking a step further. As a showcase, the enantioselective profiles of the most outstanding variants of an artificial Rh₂-based cyclopropanase (GSH, HFF and RFY) developed by Lewis and co-workers (Nat. Commun., 2015, 6, 7789 and Nat. Chem., 2018, 10, 318-324) have been rationalized. DFT calculations on the free-cofactor-mediated process identify the carbene insertion and the cyclopropanoid formation as crucial events, the latter being the enantiodetermining step, which displays up to 8 competitive orientations easily altered by the protein environment. The key intermediates of the reaction were docked into the protein scaffold showing that some mutated residues have direct interaction with the cofactor and/or the co-substrate. These interactions take the form of a direct coordination of Rh in GSH and HFF and a strong hydrophobic patch with the carbene moiety in RFY. Posterior molecular dynamics sustain that the cofactor induces global re-arrangements of the protein. Finally, massive exploration of substrate diffusion, based on the GPathFinder approach, defines this event as the origin of the enantioselectivity in GSH and RFY. For HFF, fine molecular dockings suggest that it is likely related to local interactions upon diffusion. This work shows how modelling of long-range mutations on the catalytic profiles of metalloenzymes may be unavoidable and software simulating substrate diffusion should be applied.

Mechanistic Insights into the Polymorphic Associations and Cross-Seeding of Aβ and hIAPP in the Presence of Histidine Tautomerism: An All-Atom Molecular Dynamic Study

Hundreds of millions of people around the world have been affected by Type 2 diabetes (T2D) which is a metabolic disorder. Clinical research has revealed T2D as a possible risk factor for Alzheimer’s disease (AD) development (and vice versa). Amyloid-β (Aβ) and human islet amyloid polypeptide are the main pathological species in AD and T2D, respectively. However, the mechanisms by which these two amyloidogenic peptides co-aggregate are largely uninvestigated. Herein, for the first time, we present the cross-seeding between Amylin1-37 and Aβ40 considering the particular effect of the histidine tautomerism at atomic resolution applying the all-atom molecular dynamics (MD) simulations for heterodimeric complexes. The results via random seed MD simulations indicated that the Aβ40(δδδ) isomer in cross-talking with Islet(ε) and Islet(δ) isomers could retain or increase the β-sheet content in its structure that may make it more prone to further aggregation and exhibit higher toxicity. The other tautomeric isomers which initially did not have a β-sheet structure in their monomeric forms did not show any generated β-sheet, except for one seed of the Islet(ε) and Aβ40(εεε) heterodimers complex that displayed a small amount of formed β-sheet. This computational research may provide a different point of view to examine all possible parameters that may contribute to the development of AD and T2D and provide a better understanding of the pathological link between these two severe diseases.

A Computational Study of RNA Tetraloop Thermodynamics, Including Misfolded States

An important characteristic of RNA folding is the adoption of alternative configurations of similar stability, often referred to as misfolded configurations. These configurations are considered to compete with correctly folded configurations, although their rigorous thermodynamic and structural characterization remains elusive. Tetraloop motifs found in large ribozymes are ideal systems for an atomistically detailed computational quantification of folding free energy landscapes and the structural characterization of their constituent free energy basins, including nonnative states. In this work, we studied a group of closely related 10-mer tetraloops using a combined parallel tempering and metadynamics technique that allows a reliable sampling of the free energy landscapes, requiring only knowledge that the stem folds into a canonical A-RNA configuration. We isolated and analyzed unfolded, folded, and misfolded populations that correspond to different free energy basins. We identified a distinct misfolded state that has a stability very close to that of the correctly folded state. This misfolded state contains a predominant population that shares the same structural features across all tetraloops studied here and lacks the noncanonical A-G base pair in its loop portion. Further analysis performed with biased trajectories showed that although this competitive misfolded state is not an essential intermediate, it is visited in most of the transitions from unfolded to correctly folded states. Moreover, the tetraloops can transition from this misfolded state to the correctly folded state without requiring extensive unfolding.

BitQT: a graph-based approach to the quality threshold clustering of molecular dynamics

MOTIVATION

Classical Molecular Dynamics (MD) is a standard computational approach to model time-dependent processes at the atomic level. The inherent sparsity of increasingly huge generated trajectories demands clustering algorithms to reduce other post-simulation analysis complexity. The Quality Threshold (QT) variant is an appealing one from the vast number of available clustering methods. It guarantees that all members of a particular cluster will maintain a collective similarity established by a user-defined threshold. Unfortunately, its high computational cost for processing big data limits its application in the molecular simulation field.

RESULTS

In this work, we propose a methodological parallel between QT clustering and another well-known algorithm in the field of Graph Theory, the Maximum Clique Problem. Molecular trajectories are represented as graphs whose nodes designate conformations, while unweighted edges indicate mutual similarity between nodes. The use of a binary-encoded RMSD matrix coupled to the exploitation of bitwise operations to extract clusters significantly contributes to reaching a very affordable algorithm compared to the few implementations of QT for MD available in the literature. Our alternative provides results in good agreement with the exact one while strictly preserving the collective similarity of clusters.

AVAILABILITY AND IMPLEMENTATION

The source code and documentation of BitQT are free and publicly available on GitHub (https://github.com/LQCT/BitQT.git) and ReadTheDocs (https://bitqt.readthedocs.io/en/latest/), respectively.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Reweighting of molecular simulations with explicit-solvent SAXS restraints elucidates ion-dependent RNA ensembles

Small-angle X-ray scattering (SAXS) experiments are increasingly used to probe RNA structure. A number of forward models that relate measured SAXS intensities and structural features, and that are suitable to model either explicit-solvent effects or solute dynamics, have been proposed in the past years. Here, we introduce an approach that integrates atomistic molecular dynamics simulations and SAXS experiments to reconstruct RNA structural ensembles while simultaneously accounting for both RNA conformational dynamics and explicit-solvent effects. Our protocol exploits SAXS pure-solute forward models and enhanced sampling methods to sample an heterogenous ensemble of structures, with no information towards the experiments provided on-the-fly. The generated structural ensemble is then reweighted through the maximum entropy principle so as to match reference SAXS experimental data at multiple ionic conditions. Importantly, accurate explicit-solvent forward models are used at this reweighting stage. We apply this framework to the GTPase-associated center, a relevant RNA molecule involved in protein translation, in order to elucidate its ion-dependent conformational ensembles. We show that (a) both solvent and dynamics are crucial to reproduce experimental SAXS data and (b) the resulting dynamical ensembles contain an ion-dependent fraction of extended structures.

Histidine Tautomerism Driving Human Islet Amyloid Polypeptide Aggregation in the Early Stages of Diabetes Mellitus Progression: Insight at the Atomistic Level

Early oligomerization of human islet amyloid polypeptide (hIAPP), which is accountable for β-cell death, has been implicated in the progression of type 2 diabetes mellitus. Some researches have shown the connection between hIAPP and Alzheimer's disease as well. However, the mechanism of peptide accumulation and associated cytotoxicity remains unclear. Due to the unique properties and significant role of histidine in protein sequences, here for the first time, the tautomeric effect of histidine at the early stages of amylin misfolding was investigated via molecular dynamics simulations. Considering Tau and Pi tautomeric forms of histidine (Tau and Pi tautomers are denoted as ϵ and δ, respectively), simulations were performed on two possible isomers of amylin. Our analysis revealed a higher probability of transient α-helix generation in the δ isomer in monomeric form. In dimeric forms, the δδ and δϵ conformations showed an elevated amount of α-helix and lower coil in comparison to the ϵϵ dimer. Due to the significant role of α-helix in membrane disruption and transition to β-sheet structure, these results may imply a noticeable contribution of the δ isomer and the δδ and δϵ dimers rather than ϵ and ϵϵ conformations in the early stages of diabetes initiation. Our results may aid in elucidating the hIAPP self-association process in the etiology of amyloidosis.

Hybrid Cyclobutane/Proline-Containing Peptidomimetics: The Conformational Constraint Influences Their Cell-Penetration Ability

A new family of hybrid β,γ-peptidomimetics consisting of a repetitive unit formed by a chiral cyclobutane-containing trans-β-amino acid plus a N^α-functionalized trans-γ-amino-l-proline joined in alternation were synthesized and evaluated as cell penetrating peptides (CPP). They lack toxicity on the human tumoral cell line HeLa, with an almost negligible cell uptake. The dodecapeptide showed a substantial microbicidal activity on Leishmania parasites at 50 µM but with a modest intracellular accumulation. Their previously published γ,γ-homologues, with a cyclobutane γ-amino acid, showed a well-defined secondary structure with an average inter-guanidinium distance of 8–10 Å, a higher leishmanicidal activity as well as a significant intracellular accumulation. The presence of a very rigid cyclobutane β-amino acid in the peptide backbone precludes the acquisition of a defined conformation suitable for their cell uptake ability. Our results unveiled the preorganized charge-display as a relevant parameter, additional to the separation among the charged groups as previously described. The data herein reinforce the relevance of these descriptors in the design of CPPs with improved properties.

Rounding Out the Understanding of ACD Toxicity with the Discovery of Cyclic Forms of Actin Oligomers

Actin is an essential element of both innate and adaptive immune systems and can aid in motility and translocation of bacterial pathogens, making it an attractive target for bacterial toxins. Pathogenic Vibrio and Aeromonas genera deliver actin cross-linking domain (ACD) toxin into the cytoplasm of the host cell to poison actin regulation and promptly induce cell rounding. At early stages of toxicity, ACD covalently cross-links actin monomers into oligomers (AOs) that bind through multivalent interactions and potently inhibit several families of actin assembly proteins. At advanced toxicity stages, we found that the terminal protomers of linear AOs can get linked together by ACD to produce cyclic AOs. When tested against formins and Ena/VASP, linear and cyclic AOs exhibit similar inhibitory potential, which for the cyclic AOs is reduced in the presence of profilin. In coarse-grained molecular dynamics simulations, profilin and WH2-motif binding sites on actin subunits remain exposed in modeled AOs of both geometries. We speculate, therefore, that the reduced toxicity of cyclic AOs is due to their reduced configurational entropy. A characteristic feature of cyclic AOs is that, in contrast to the linear forms, they cannot be straightened to form filaments (e.g., through stabilization by cofilin), which makes them less susceptible to neutralization by the host cell.

Rationalizing the Decavanadate(V) and Oxidovanadium(IV) Binding to G-Actin and the Competition with Decaniobate(V) and ATP

The experimental data collected over the past 15 years on the interaction of decavanadate(V) (V₁₀O₂₈^6-; V₁₀), a polyoxometalate (POM) with promising anticancer and antibacterial action, with G-actin, were rationalized by using several computational approaches (docking, density functional theory (DFT), and molecular dynamics (MD)). Moreover, a comparison with the isostructural and more stable decaniobate(V) (Nb₁₀O₂₈^6-; Nb₁₀) was carried out. Four binding sites were identified, named α, β, γ, and δ, the site α being the catalytic nucleotide site located in the cleft of the enzyme at the interface of the subdomains II and IV. It was observed that the site α is preferred by V₁₀, whereas Nb₁₀ is more stable at the site β; this indicates that, differently from other proteins, G-actin could contemporaneously bind the two POMs, whose action would be synergistic. Both decavanadate and decaniobate induce conformational rearrangements in G-actin, larger for V₁₀ than Nb₁₀. Moreover, the binding mode of oxidovanadium(IV) ion, V^IVO²⁺, formed upon the reduction of decavanadate(V) by the -SH groups of accessible cysteine residues, is also found in the catalytic site α with (His161, Asp154) coordination; this adduct overlaps significantly with the region where ATP is bound, accounting for the competition between V₁₀ and its reduction product V^IVO²⁺ with ATP, as previously observed by EPR spectroscopy. Finally, the competition with ATP was rationalized: since decavanadate prefers the nucleotide site α, Ca²⁺-ATP displaces V₁₀ from this site, while the competition is less important for Nb₁₀ because this POM shows a higher affinity for β than for site α. A relevant consequence of this paper is that other metallodrug-protein systems, in the absence or presence of eventual inhibitors and/or competition with molecules of the organism, could be studied with the same approach, suggesting important elements for an explanation of the biological data and a rational drug design.

Investigating Potential Inhibitory Effect of Uncaria tomentosa (Cat's Claw) against the Main Protease 3CLpro of SARS-CoV-2 by Molecular Modeling

COVID-19 is a disease caused by severe acute respiratory syndrome coronavirus 2. Presently, there is no effective treatment for COVID-19. As part of the worldwide efforts to find efficient therapies and preventions, it has been reported the crystalline structure of the SARS-CoV-2 main protease Mpro (also called 3CLpro) bound to a synthetic inhibitor, which represents a major druggable target. The druggability of Mpro could be used for discovering drugs to treat COVID-19. A multilevel computational study was carried out to evaluate the potential antiviral properties of the components of the medicinal herb Uncaria tomentosa (Cat's claw), focusing on the inhibition of Mpro. The in silico approach starts with protein-ligand docking of 26 Cat's claw key components, followed by ligand pathway calculations, molecular dynamics simulations, and MM-GBSA calculation of the free energy of binding for the best docked candidates. The structural bioinformatics approaches led to identification of three bioactive compounds of Uncaria tomentosa (speciophylline, cadambine, and proanthocyanidin B2) with potential therapeutic effects by strong interaction with 3CLpro. Additionally, in silico drug-likeness indices for these components were calculated and showed good predicted therapeutic profiles of these phytochemicals. Our findings suggest the potential effectiveness of Cat's claw as complementary and/or alternative medicine for COVID-19 treatment.

Ligand Pose Predictions for Human G Protein-Coupled Receptors: Insights from the Amber-Based Hybrid Molecular Mechanics/Coarse-Grained Approach

Human G protein-coupled receptors (hGPCRs) are the most frequent targets of Food and Drug Administration (FDA)-approved drugs. Structural bioinformatics, along with molecular simulation, can support structure-based drug design targeting hGPCRs. In this context, several years ago, we developed a hybrid molecular mechanics (MM)/coarse-grained (CG) approach to predict ligand poses in low-resolution hGPCR models. The approach was based on the GROMOS96 43A1 and PRODRG united-atom force fields for the MM part. Here, we present a new MM/CG implementation using, instead, the Amber 14SB and GAFF all-atom potentials for proteins and ligands, respectively. The new implementation outperforms the previous one, as shown by a variety of applications on models of hGPCR/ligand complexes at different resolutions, and it is also more user-friendly. Thus, it emerges as a useful tool to predict poses in low-resolution models and provides insights into ligand binding similarly to all-atom molecular dynamics, albeit at a lower computational cost.

Commensal bacteria stimulate antitumor responses via T cell cross-reactivity

Recent studies show gut microbiota modulate antitumor immune responses; one proposed mechanism is cross-reactivity between antigens expressed in commensal bacteria and neoepitopes. We found that T cells targeting an epitope called SVYRYYGL (SVY), expressed in the commensal bacterium Bifidobacterium breve (B. breve), cross-react with a model neoantigen, SIYRYYGL (SIY). Mice lacking B. breve had decreased SVY-reactive T cells compared with B. breve-colonized mice, and the T cell response was transferable by SVY immunization or by cohousing mice without Bifidobacterium with ones colonized with Bifidobacterium. Tumors expressing the model SIY neoantigen also grew faster in mice lacking B. breve compared with Bifidobacterium-colonized animals. B. breve colonization also shaped the SVY-reactive TCR repertoire. Finally, SVY-specific T cells recognized SIY-expressing melanomas in vivo and led to decreased tumor growth and extended survival. Our work demonstrates that commensal bacteria can stimulate antitumor immune responses via cross-reactivity and how bacterial antigens affect the T cell landscape.

BitClust: Fast Geometrical Clustering of Long Molecular Dynamics Simulations

The growing computational capacity allows the investigation of large biomolecular systems by increasingly extensive molecular dynamics simulations. The resulting huge trajectories demand efficient partition methods to discern relevant structural dissimilarity. Clustering algorithms are available to address this task, but their implementations still need to be improved to gain in computational speed and to reduce the consumption of random access memory. We propose the BitClust code which, based on a combination of Python and C programming languages, performs fast structural clustering of long molecular trajectories. BitClust takes advantage of bitwise operations applied to a bit-encoded pairwise similarity matrix. Our approach allowed us to process a half-million frame trajectory in 6 h using less than 35 GB, a task that is not affordable with any of the similar alternatives.