Computational methods in drug discovery significantly reduce both time and experimental costs. Nonetheless, certain computational tasks in drug discovery can be daunting with classical computing techniques which can be potentially overcome using quantum computing. A crucial task within this domain involves the functional classification of proteins. However, a challenge lies in adequately representing lengthy protein sequences given the limited number of qubits available in existing noisy quantum computers. We show that protein sequences can be thought of as sentences in natural language processing and can be parsed using the existing Quantum Natural Language framework into parameterized quantum circuits of reasonable qubits, which can be trained to solve various protein-related machine-learning problems. We classify proteins based on their sub-cellular locations—a pivotal task in bioinformatics that is key to understanding biological processes and disease mechanisms. Leveraging the quantum-enhanced processing capabilities, we demonstrate that Quantum Tensor Networks (QTN) can effectively handle the complexity and diversity of protein sequences. We present a detailed methodology that adapts QTN architectures to the nuanced requirements of protein data, supported by comprehensive experimental results. We demonstrate two distinct QTNs, inspired by classical recurrent neural networks (RNN) and convolutional neural networks (CNN), to solve the binary classification task mentioned above. Our top-performing quantum model has achieved a 94% accuracy rate, which is comparable to the performance of a classical model that uses the ESM2 protein language model embeddings. It’s noteworthy that the ESM2 model is extremely large, containing 8 million parameters in its smallest configuration, whereas our best quantum model requires only around 800 parameters. We demonstrate that these hybrid models exhibit promising performance, showcasing their potential to compete with classical models of similar complexity.

Molecular dynamics (MD) is the primary computational method by which modern structural biology explores macromolecule structure and function. Boltzmann generators have been proposed as an alternative to MD, by replacing the integration of molecular systems over time with the training of generative neural networks. This neural network approach to MD enables convergence to thermodynamic equilibrium faster than traditional MD; however, critical gaps in the theory and computational feasibility of Boltzmann generators significantly reduce their usability. Here, we develop a mathematical foundation to overcome these barriers; we demonstrate that the Boltzmann generator approach is sufficiently rapid to replace traditional MD for complex macromolecules, such as proteins in specific applications, and we provide a comprehensive toolkit for the exploration of molecular energy landscapes with neural networks.

Traditional drug discovery pipelines can require multiple years and billions of dollars of investment. Deep generative and discriminative models are widely adopted to assist in drug development. Classical machines cannot efficiently reproduce the atypical patterns of quantum computers, which may improve the quality of learned tasks. We propose a suite of quantum machine learning techniques: incorporating generative adversarial networks (GAN), convolutional neural networks (CNN) and variational auto-encoders (VAE) to generate small drug molecules, classify binding pockets in proteins, and generate large drug molecules, respectively.

Combining structural proteomics experimental data with computational methods is a powerful tool for protein structure prediction. Here, we apply a recently-developed approach for de novo protein structure determination based on the incorporation of short-distance crosslinking data as constraints in discrete molecular dynamics simulations (CL-DMD) for the determination of conformational ensemble of the intrinsically disordered protein α-synuclein in the solution. The predicted structures were in agreement with hydrogen-deuterium exchange, circular dichroism, surface modification, and long-distance crosslinking data. We found that α-synuclein is present in solution as an ensemble of rather compact globular conformations with distinct topology and inter-residue contacts, which is well-represented by movements of the large loops and formation of few transient secondary structure elements. Non-amyloid component and C-terminal regions were consistently found to contain β-structure elements and hairpins.

Exposure to the environmental toxin β-methylamino-L-alanine (BMAA) is linked to amyotrophic lateral sclerosis (ALS), but its disease-promoting mechanism remains unknown. We propose that incorporation of BMAA into the ALS-linked protein Cu,Zn superoxide dismutase (SOD1) upon translation promotes protein misfolding and aggregation, which has been linked to ALS onset and progression. Using molecular simulation and predictive energetic computation, we demonstrate that substituting any serine with BMAA in SOD1 results in structural destabilization and aberrant dynamics, promoting neurotoxic SOD1 aggregation. We propose that translational incorporation of BMAA into SOD1 is directly responsible for its toxicity in neurodegeneration, and BMAA modification of SOD1 may serve as a biomarker of ALS.

Mutations in the Cystic Fibrosis Transmembrane Conductance Regulator (CFTR) gene affect CFTR protein biogenesis or its function as a chloride channel, resulting in dysregulation of epithelial fluid transport in the lung, pancreas and other organs in cystic fibrosis (CF). Development of pharmaceutical strategies to treat CF requires understanding of the mechanisms underlying channel function. However, incomplete 3D structural information on the unique ABC ion channel, CFTR, hinders elucidation of its functional mechanism and correction of cystic fibrosis causing mutants. Several CFTR homology models have been developed using bacterial ABC transporters as templates but these have low sequence similarity to CFTR and are not ion channels. Here, we refine an earlier model in an outward (OWF) and develop an inward (IWF) facing model employing an integrated experimental-molecular dynamics simulation (200 ns) approach. Our IWF structure agrees well with a recently solved cryo-EM structure of a CFTR IWF state. We utilize cysteine cross-linking to verify positions and orientations of residues within trans-membrane helices (TMHs) of the OWF conformation and to reconstruct a physiologically relevant pore structure. Comparison of pore profiles of the two conformations reveal a radius sufficient to permit passage of hydrated Cl- ions in the OWF but not the IWF model. To identify structural determinants that distinguish the two conformations and possible rearrangements of TMHs within them responsible for channel gating, we perform cross-linking by bifunctional reagents of multiple predicted pairs of cysteines in TMH 6 and 12 and 6 and 9. To determine whether the effects of cross-linking on gating observed are the result of switching of the channel from open to close state, we also treat the same residue pairs with monofunctional reagents in separate experiments. Both types of reagents prevent ion currents indicating that pore blockage is primarily responsible.

UNLABELLED: A key to understanding RNA function is to uncover its complex 3D structure. Experimental methods used for determining RNA 3D structures are technologically challenging and laborious, which makes the development of computational prediction methods of substantial interest. Previously, we developed the iFoldRNA server that allows accurate prediction of short (<50 nt) tertiary RNA structures starting from primary sequences. Here, we present a new version of the iFoldRNA server that permits the prediction of tertiary structure of RNAs as long as a few hundred nucleotides. This substantial increase in the server capacity is achieved by utilization of experimental information such as base-pairing and hydroxyl-radical probing. We demonstrate a significant benefit provided by integration of experimental data and computational methods. AVAILABILITY AND IMPLEMENTATION: http://ifoldrna.dokhlab.org CONTACT: [email protected].

The increased risk of developing Alzheimer's disease (AD) is associated with the APOE gene, which encodes for three variants of Apolipoprotein E, namely E2, E3, E4, differing only by two amino acids at positions 112 and 158. ApoE4 is known to be the strongest risk factor for AD onset, while ApoE3 and ApoE2 are considered to be the AD-neutral and AD-protective isoforms, respectively. It has been hypothesized that the ApoE isoforms may contribute to the development of AD by modifying the homeostasis of ApoE physiological partners and AD-related proteins in an isoform-specific fashion. Here we find that, despite the high sequence similarity among the three ApoE variants, only ApoE4 exhibits a misfolded intermediate state characterized by isoform-specific domain-domain interactions in molecular dynamics simulations. The existence of an ApoE4-specific intermediate state can contribute to the onset of AD by altering multiple cellular pathways involved in ApoE-dependent lipid transport efficiency or in AD-related protein aggregation and clearance. We present what we believe to be the first structural model of an ApoE4 misfolded intermediate state, which may serve to elucidate the molecular mechanism underlying the role of ApoE4 in AD pathogenesis. The knowledge of the structure for the ApoE4 folding intermediate provides a new platform for the rational design of alternative therapeutic strategies to fight AD.

MOTIVATION: Increasing use of structural modeling for understanding structure-function relationships in proteins has led to the need to ensure that the protein models being used are of acceptable quality. Quality of a given protein structure can be assessed by comparing various intrinsic structural properties of the protein to those observed in high-resolution protein structures. RESULTS: In this study, we present tools to compare a given structure to high-resolution crystal structures. We assess packing by calculating the total void volume, the percentage of unsatisfied hydrogen bonds, the number of steric clashes and the scaling of the accessible surface area. We assess covalent geometry by determining bond lengths, angles, dihedrals and rotamers. The statistical parameters for the above measures, obtained from high-resolution crystal structures enable us to provide a quality-score that points to specific areas where a given protein structural model needs improvement. AVAILABILITY AND IMPLEMENTATION: We provide these tools that appraise protein structures in the form of a web server Gaia (http://chiron.dokhlab.org). Gaia evaluates the packing and covalent geometry of a given protein structure and provides quantitative comparison of the given structure to high-resolution crystal structures. CONTACT: [email protected] SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

MOTIVATION: Identifying the location of binding sites on proteins is of fundamental importance for a wide range of applications, including molecular docking, de novo drug design, structure identification and comparison of functional sites. Here we present Erebus, a web server that searches the entire Protein Data Bank for a given substructure defined by a set of atoms of interest, such as the binding scaffolds for small molecules. The identified substructure contains atoms having the same names, belonging to same amino acids and separated by the same distances (within a given tolerance) as the atoms of the query structure. The accuracy of a match is measured by the root-mean-square deviation or by the normal weight with a given variance. Tests show that our approach can reliably locate rigid binding scaffolds of drugs and metal ions. AVAILABILITY AND IMPLEMENTATION: We provide this service through a web server at http://erebus.dokhlab.org.

Determining the forces that conserve amino acid positions in proteins across species is a fundamental pursuit of molecular evolution. Evolutionary conservation is driven by either a protein's function or its thermodynamic stability. Highly conserved histone proteins offer a platform to evaluate these driving forces. While the conservation of histone H3 and H4 "tail" domains and surface residues are driven by functional importance, the driving force behind the conservation of buried histone residues has not been examined. Using a computational approach, we determined the thermodynamically preferred amino acids at each buried position in H3 and H4. In agreement with what is normally observed in proteins, we find a significant correlation between thermodynamic stability and evolutionary conservation in the buried residues in H4. In striking contrast, we find that thermodynamic stability of buried H3 residues does not correlate with evolutionary conservation. Given that these H3 residues are not post-translationally modified and only regulate H3-H3 and H3-H4 stabilizing interactions, our data imply an unknown function responsible for driving conservation of these buried H3 residues.

Polyglutamine (polyQ) expansion in exon1 (XN1) of the huntingtin protein is linked to Huntington's disease. When the number of glutamines exceeds a threshold of approximately 36-40 repeats, XN1 can readily form amyloid aggregates similar to those associated with disease. Many experiments suggest that misfolding of monomeric XN1 plays an important role in the length-dependent aggregation. Elucidating the misfolding of a XN1 monomer can help determine the molecular mechanism of XN1 aggregation and potentially help develop strategies to inhibit XN1 aggregation. The flanking sequences surrounding the polyQ region can play a critical role in determining the structural rearrangement and aggregation mechanism of XN1. Few experiments have studied XN1 in its entirety, with all flanking regions. To obtain structural insights into the misfolding of XN1 toward amyloid aggregation, we perform molecular dynamics simulations on monomeric XN1 with full flanking regions, a variant missing the polyproline regions, which are hypothesized to prevent aggregation, and an isolated polyQ peptide (Q(n)). For each of these three constructs, we study glutamine repeat lengths of 23, 36, 40 and 47. We find that polyQ peptides have a positive correlation between their probability to form a beta-rich misfolded state and their expansion length. We also find that the flanking regions of XN1 affect its probability to form a beta-rich state compared to the isolated polyQ. Particularly, the polyproline regions form polyproline type II helices and decrease the probability of the polyQ region to form a beta-rich state. Additionally, by lengthening polyQ, the first N-terminal 17 residues are more likely to adopt a beta-sheet conformation rather than an alpha-helix conformation. Therefore, our molecular dynamics study provides a structural insight of XN1 misfolding and elucidates the possible role of the flanking sequences in XN1 aggregation.

Ryanodine receptors (RyRs) are ion channels that regulate muscle contraction by releasing calcium ions from intracellular stores into the cytoplasm. Mutations in skeletal muscle RyR (RyR1) give rise to congenital diseases such as central core disease. The absence of high-resolution structures of RyR1 has limited our understanding of channel function and disease mechanisms at the molecular level. Here, we report a structural model of the pore-forming region of RyR1. Molecular dynamics simulations show high ion binding to putative pore residues D4899, E4900, D4938, and D4945, which are experimentally known to be critical for channel conductance and selectivity. We also observe preferential localization of Ca(2+) over K(+) in the selectivity filter of RyR1. Simulations of RyR1-D4899Q mutant show a loss of preference to Ca(2+) in the selectivity filter as seen experimentally. Electrophysiological experiments on a central core disease mutant, RyR1-G4898R, show constitutively open channels that conduct K(+) but not Ca(2+). Our simulations with G4898R likewise show a decrease in the preference of Ca(2+) over K(+) in the selectivity filter. Together, the computational and experimental results shed light on ion conductance and selectivity of RyR1 at an atomistic level.

UNLABELLED: Three-dimensional RNA structure prediction and folding is of significant interest in the biological research community. Here, we present iFoldRNA, a novel web-based methodology for RNA structure prediction with near atomic resolution accuracy and analysis of RNA folding thermodynamics. iFoldRNA rapidly explores RNA conformations using discrete molecular dynamics simulations of input RNA sequences. Starting from simplified linear-chain conformations, RNA molecules (<50 nt) fold to native-like structures within half an hour of simulation, facilitating rapid RNA structure prediction. All-atom reconstruction of energetically stable conformations generates iFoldRNA predicted RNA structures. The predicted RNA structures are within 2-5 A root mean squre deviations (RMSDs) from corresponding experimentally derived structures. RNA folding parameters including specific heat, contact maps, simulation trajectories, gyration radii, RMSDs from native state, fraction of native-like contacts are accessible from iFoldRNA. We expect iFoldRNA will serve as a useful resource for RNA structure prediction and folding thermodynamic analyses. AVAILABILITY: http://iFoldRNA.dokhlab.org.

Cysteine-rich intestinal protein 1 (CRIP1) has been identified as a novel marker for early detection of cancers. Here we report on the use of phage display in combination with molecular modeling to identify a high-affinity ligand for CRIP1. Panning experiments using a circularized C7C phage library yielded several consensus sequences with modest binding affinities to purified CRIP1. Two sequence motifs, A1 and B5, having the highest affinities for CRIP1, were chosen for further study. With peptide structure information and the NMR structure of CRIP1, the higher-affinity A1 peptide was computationally redesigned, yielding a novel peptide, A1M, whose affinity was predicted to be much improved. Synthesis of the peptide and saturation and competitive binding studies demonstrated approximately a 10-28-fold improvement in the affinity of A1M compared to that of either A1 or B5 peptide. These techniques have broad application to the design of novel ligand peptides.

The absence of a functional ATP Binding Cassette (ABC) protein called the Cystic Fibrosis Transmembrane Conductance Regulator (CFTR) from apical membranes of epithelial cells is responsible for cystic fibrosis (CF). Over 90% of CF patients carry at least one mutant allele with deletion of phenylalanine at position 508 located in the N-terminal nucleotide binding domain (NBD1). Biochemical and cell biological studies show that the DeltaF508 mutant exhibits inefficient biosynthetic maturation and susceptibility to degradation probably due to misfolding of NBD1 and the resultant misassembly of other domains. However, little is known about the direct effect of the Phe508 deletion on the NBD1 folding, which is essential for rational design strategies of cystic fibrosis treatment. Here we show that the deletion of Phe508 alters the folding dynamics and kinetics of NBD1, thus possibly affecting the assembly of the complete CFTR. Using molecular dynamics simulations, we find that meta-stable intermediate states appearing on wild type and mutant folding pathways are populated differently and that their kinetic accessibilities are distinct. The structural basis of the increased misfolding propensity of the DeltaF508 NBD1 mutant is the perturbation of interactions in residue pairs Q493/P574 and F575/F578 found in loop S7-H6. As a proof-of-principle that the S7-H6 loop conformation can modulate the folding kinetics of NBD1, we virtually design rescue mutations in the identified critical interactions to force the S7-H6 loop into the wild type conformation. Two redesigned NBD1-DeltaF508 variants exhibited significantly higher folding probabilities than the original NBD1-DeltaF508, thereby partially rescuing folding ability of the NBD1-DeltaF508 mutant. We propose that these observed defects in folding kinetics of mutant NBD1 may also be modulated by structures separate from the 508 site. The identified structural determinants of increased misfolding propensity of NBD1-DeltaF508 are essential information in correcting this pathogenic mutant.

Nuclear receptor ligand binding domains (LBDs) convert ligand binding events into changes in gene expression by recruiting transcriptional coregulators to a conserved activation function-2 (AF-2) surface. While most nuclear receptor LBDs form homo- or heterodimers, the human nuclear receptor pregnane X receptor (PXR) forms a unique and essential homodimer and is proposed to assemble into a functional heterotetramer with the retinoid X receptor (RXR). How the homodimer interface, which is located 30 A from the AF-2, would affect function at this critical surface has remained unclear. By using 20- to 30-ns molecular dynamics simulations on PXR in various oligomerization states, we observed a remarkably high degree of correlated motion in the PXR-RXR heterotetramer, most notably in the four helices that create the AF-2 domain. The function of such correlation may be to create "active-capable" receptor complexes that are ready to bind to transcriptional coactivators. Indeed, we found in additional simulations that active-capable receptor complexes involving other orphan or steroid nuclear receptors also exhibit highly correlated AF-2 domain motions. We further propose a mechanism for the transmission of long-range motions through the nuclear receptor LBD to the AF-2 surface. Taken together, our findings indicate that long-range motions within the LBD scaffold are critical to nuclear receptor function by promoting a mobile AF-2 state ready to bind coactivators.

UNLABELLED: We built a novel web-based platform for performing discrete molecular dynamics simulations of proteins. In silico protein folding involves searching for minimal frustration in the vast conformational landscape. Conventional approaches for simulating protein folding insufficiently address the problem of simulations in relevant time and length scales necessary for a mechanistic understanding of underlying biomolecular phenomena. Discrete molecular dynamics (DMD) offers an opportunity to bridge the size and timescale gaps and uncover the structural and biological properties of experimentally undetectable protein dynamics. The iFold server supports large-scale simulations of protein folding, thermal denaturation, thermodynamic scan, simulated annealing and p(fold) analysis using DMD and coarse-grained protein model with structure-based Gō-interactions between amino acids. AVAILABILITY: http://ifold.dokhlab.org

Diverse proteins with similar structures are grouped into families of homologs and analogs, if their sequence similarity is higher or lower, respectively, than 20%-30%. It was suggested that protein homologs and analogs originate from a common ancestor and diverge in their distinct evolutionary time scales, emerging as a consequence of the physical properties of the protein sequence space. Although a number of studies have determined key signatures of protein family organization, the sequence-structure factors that differentiate the two evolution-related protein families remain unknown. Here, we stipulate that subtle structural changes, which appear due to accumulating mutations in the homologous families, lead to distinct packing of the protein core and, thus, novel compositions of core residues. The latter process leads to the formation of distinct families of homologs. We propose that such differentiation results in the formation of analogous families. To test our postulate, we developed a molecular modeling and design toolkit, Medusa, to computationally design protein sequences that correspond to the same fold family. We find that analogous proteins emerge when a backbone structure deviates only 1-2 angstroms root-mean-square deviation from the original structure. For close homologs, core residues are highly conserved. However, when the overall sequence similarity drops to approximately 25%-30%, the composition of core residues starts to diverge, thereby forming novel families of protein homologs. This direct observation of the formation of protein homologs within a specific fold family supports our hypothesis. The conservation of amino acids in designed sequences recapitulates that of the naturally occurring sequences, thereby validating our computational design methodology.

In PLoS Computational Biology, volume 2, issue 7: DOI: 10.1371/journal.pcbi.0020085 The references to the figure parts in the legend of figure 3 were incorrect. The correct caption is as follows: Figure 3. The Sequence Identity for the Constructed Homologous Structures Three different protein folds are studied: HPR domain (A,D), ROSSMAN fold (B,E), and SH3 domain (C,F). (A,B,C) The sequence identities of the redesigned proteins using the flexible-backbone design simulation are presented as the function of the backbone-RMSD from the reference protein. (D,E,F) The sequence identity of the core is also plotted against the overall sequence identity. The “twilight zone” of sequence identity (20%–30%) corresponds to regions between horizontal (A,B,C) or vertical (D,E,F) lines.

Expansion of polyglutamine (polyQ) tracts in proteins results in protein aggregation and is associated with cell death in at least nine neurodegenerative diseases. Disease age of onset is correlated with the polyQ insert length above a critical value of 35-40 glutamines. The aggregation kinetics of isolated polyQ peptides in vitro also shows a similar critical-length dependence. While recent experimental work has provided considerable insights into polyQ aggregation, the molecular mechanism of aggregation is not well understood. Here, using computer simulations of isolated polyQ peptides, we show that a mechanism of aggregation is the conformational transition in a single polyQ peptide chain from random coil to a parallel beta-helix. This transition occurs selectively in peptides longer than 37 glutamines. In the beta-helices observed in simulations, all residues adopt beta-strand backbone dihedral angles, and the polypeptide chain coils around a central helical axis with 18.5 +/- 2 residues per turn. We also find that mutant polyQ peptides with proline-glycine inserts show formation of antiparallel beta-hairpins in their ground state, in agreement with experiments. The lower stability of mutant beta-helices explains their lower aggregation rates compared to wild type. Our results provide a molecular mechanism for polyQ-mediated aggregation.