Computational optimization of the SARS-CoV-2 receptor-binding-motif affinity for human ACE2.

The coronavirus severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which is responsible for the coronavirus disease 2019 pandemic, and the closely related SARS-CoV coronavirus enter cells by binding at the human angiotensin converting enzyme 2 (hACE2). The stronger hACE2 affinity of SARS-CoV-2 has been connected with its higher infectivity. In this work, we study hACE2 complexes with the receptor-binding domains (RBDs) of the human SARS-CoV-2 and human SARS-CoV viruses, using all-atom molecular dynamics simulations and computational protein design with a physics-based energy function. The molecular dynamics simulations identify charge-modifying substitutions between the CoV-2 and CoV RBDs, which either increase or decrease the hACE2 affinity of the SARS-CoV-2 RBD. The combined effect of these mutations is small, and the relative affinity is mainly determined by substitutions at residues in contact with hACE2. Many of these findings are in line and interpret recent experiments. Our computational protein design calculations redesign positions 455, 493, 494, and 501 of the SARS-CoV-2 receptor binding motif, which contact hACE2 in the complex and are important for ACE2 recognition. Sampling is enhanced by an adaptive importance sampling Monte Carlo method. Sequences with increased affinity replace CoV-2 glutamine by a negative residue at position 493; serine by a nonpolar or aromatic residue or an asparagine at position 494; and asparagine by valine or threonine at position 501. Substitutions at positions 455 and 501 have a smaller effect on affinity. Substitutions suggested by our design are seen in viral sequences encountered in other species, including bat and pangolin. Our results might be used to identify potential virus strains with higher human infectivity and assist in the design of peptide-based or peptidomimetic compounds with the potential to inhibit SARS-CoV-2 binding at hACE2.

Recognition of LD motifs by the focal adhesion targeting domains of focal adhesion kinase and proline-rich tyrosine kinase 2-beta: Insights from molecular dynamics simulations.

The focal adhesion kinase (FAK) and the proline-rich tyrosine kinase 2-beta (PYK2) are implicated in cancer progression and metastasis and represent promising biomarkers and targets for cancer therapy. FAK and PYK2 are recruited to focal adhesions (FAs) via interactions between their FA targeting (FAT) domains and conserved segments (LD motifs) on the proteins Paxillin, Leupaxin, and Hic-5. A promising new approach for the inhibition of FAK and PYK2 targets interactions of the FAK domains with proteins that promote localization at FAs. Advances toward this goal include the development of surface plasmon resonance, heteronuclear single quantum coherence nuclear magnetic resonance (HSQC-NMR) and fluorescence polarization assays for the identification of fragments or compounds interfering with the FAK-Paxillin interaction. We have recently validated this strategy, showing that Paxillin mimicking polypeptides with 2 to 3 LD motifs displace FAK from FAs and block kinase-dependent and independent functions of FAK, including downstream integrin signaling and FA localization of the protein p130Cas. In the present work we study by all-atom molecular dynamics simulations the recognition of peptides with the Paxillin and Leupaxin LD motifs by the FAK-FAT and PYK2-FAT domains. Our simulations and free-energy analysis interpret experimental data on binding of Paxillin and Leupaxin LD motifs at FAK-FAT and PYK2-FAT binding sites, and assess the roles of consensus LD regions and flanking residues. Our results can assist in the design of effective inhibitory peptides of the FAK-FAT: Paxillin and PYK2-FAT:Leupaxin complexes and the construction of pharmacophore models for the discovery of potential small-molecule inhibitors of the FAK-FAT and PYK2-FAT focal adhesion based functions.

Physics-Based Computational Protein Design: An Update.

We describe methods for physics-based protein design and some recent applications from our work. We present the physical interpretation of a MC simulation in sequence space and show that sequences and conformations form a well-defined statistical ensemble, explored with Monte Carlo and Boltzmann sampling. The folded state energy combines molecular mechanics for solutes with continuum electrostatics for solvent. We usually assume one or a few fixed protein backbone structures and discrete side chain rotamers. Methods based on molecular dynamics, which introduce additional backbone and side chain flexibility, are under development. The redesign of a PDZ domain and an aminoacyl-tRNA synthetase enzyme were successful. We describe a versatile, adaptive, Wang-Landau MC method that can be used to design for substrate affinity, catalytic rate, catalytic efficiency, or the specificity of these properties. The methods are transferable to all biomolecules, can be systematically improved, and give physical insights.

Hybrid MC/MD for protein design.

Computational protein design relies on simulations of a protein structure, where selected amino acids can mutate randomly, and mutations are selected to enhance a target property, such as stability. Often, the protein backbone is held fixed and its degrees of freedom are modeled implicitly to reduce the complexity of the conformational space. We present a hybrid method where short molecular dynamics (MD) segments are used to explore conformations and alternate with Monte Carlo (MC) moves that apply mutations to side chains. The backbone is fully flexible during MD. As a test, we computed side chain acid/base constants or pKa's in five proteins. This problem can be considered a special case of protein design, with protonation/deprotonation playing the role of mutations. The solvent was modeled as a dielectric continuum. Due to cost, in each protein we allowed just one side chain position to change its protonation state and the other position to change its type or mutate. The pKa's were computed with a standard method that scans a range of pH values and with a new method that uses adaptive landscape flattening (ALF) to sample all protonation states in a single simulation. The hybrid method gave notably better accuracy than standard, fixed-backbone MC. ALF decreased the computational cost a factor of 13.

Comparing pairwise-additive and many-body generalized Born models for acid/base calculations and protein design.

Generalized Born (GB) solvent models are common in acid/base calculations and protein design. With GB, the interaction between a pair of solute atoms depends on the shape of the protein/solvent boundary and, therefore, the positions of all solute atoms, so that GB is a many-body potential. For compute-intensive applications, the model is often simplified further, by introducing a mean, native-like protein/solvent boundary, which removes the many-body property. We investigate a method for both acid/base calculations and protein design that uses Monte Carlo simulations in which side chains can explore rotamers, bind/release protons, or mutate. The fluctuating protein/solvent dielectric boundary is treated in a way that is numerically exact (within the GB framework), in contrast to a mean boundary. Its originality is that it captures the many-body character while retaining the residue-pairwise complexity given by a fixed boundary. The method is implemented in the Proteus protein design software. It yields a slight but systematic improvement for acid/base constants in nine proteins and a significant improvement for the computational design of three PDZ domains. It eliminates a source of model uncertainty, which will facilitate the analysis of other model limitations. © 2017 Wiley Periodicals, Inc.

Simple models for nonpolar solvation: Parameterization and testing.

Implicit solvent models are important for many biomolecular simulations. The polarity of aqueous solvent is essential and qualitatively captured by continuum electrostatics methods like Generalized Born (GB). However, GB does not account for the solvent-induced interactions between exposed hydrophobic sidechains or solute-solvent dispersion interactions. These "nonpolar" effects are often modeled through surface area (SA) energy terms, which lack realism, create mathematical singularities, and have a many-body character. We have explored an alternate, Lazaridis-Karplus (LK) gaussian energy density for nonpolar effects and a dispersion (DI) energy term proposed earlier, associated with GB electrostatics. We parameterized several combinations of GB, SA, LK, and DI energy terms, to reproduce 62 small molecule solvation free energies, 387 protein stability changes due to point mutations, and the structures of 8 protein loops. With optimized parameters, the models all gave similar results, with GBLK and GBDILK giving no performance loss compared to GBSA, and mean errors of 1.7 kcal/mol for the stability changes and 2 Å deviations for the loop conformations. The optimized GBLK model gave poor results in MD of the Trpcage mini-protein, but parameters optimized specifically for MD performed well for Trpcage and three other small proteins. Overall, the LK and DI nonpolar terms are valid alternatives to SA treatments for a range of applications. © 2017 Wiley Periodicals, Inc.

Monte Carlo simulations of proteins at constant pH with generalized Born solvent, flexible sidechains, and an effective dielectric boundary.

Titratable residues determine the acid/base behavior of proteins, strongly influencing their function; in addition, proton binding is a valuable reporter on electrostatic interactions. We describe a method for pK(a) calculations, using constant-pH Monte Carlo (MC) simulations to explore the space of sidechain conformations and protonation states, with an efficient and accurate generalized Born model (GB) for the solvent effects. To overcome the many-body dependency of the GB model, we use a "Native Environment" approximation, whose accuracy is shown to be good. It allows the precalculation and storage of interactions between all sidechain pairs, a strategy borrowed from computational protein design, which makes the MC simulations themselves very fast. The method is tested for 12 proteins and 167 titratable sidechains. It gives an rms error of 1.1 pH units, similar to the trivial "Null" model. The only adjustable parameter is the protein dielectric constant. The best accuracy is achieved for values between 4 and 8, a range that is physically plausible for a protein interior. For sidechains with large pKa shifts, ≥2, the rms error is 1.6, compared to 2.5 with the Null model and 1.5 with the empirical PROPKA method.

Computational protein design: the Proteus software and selected applications.

We describe an automated procedure for protein design, implemented in a flexible software package, called Proteus. System setup and calculation of an energy matrix are done with the XPLOR modeling program and its sophisticated command language, supporting several force fields and solvent models. A second program provides algorithms to search sequence space. It allows a decomposition of the system into groups, which can be combined in different ways in the energy function, for both positive and negative design. The whole procedure can be controlled by editing 2-4 scripts. Two applications consider the tyrosyl-tRNA synthetase enzyme and its successful redesign to bind both O-methyl-tyrosine and D-tyrosine. For the latter, we present Monte Carlo simulations where the D-tyrosine concentration is gradually increased, displacing L-tyrosine from the binding pocket and yielding the binding free energy difference, in good agreement with experiment. Complete redesign of the Crk SH3 domain is presented. The top 10000 sequences are all assigned to the correct fold by the SUPERFAMILY library of Hidden Markov Models. Finally, we report the acid/base behavior of the SNase protein. Sidechain protonation is treated as a form of mutation; it is then straightforward to perform constant-pH Monte Carlo simulations, which yield good agreement with experiment. Overall, the software can be used for a wide range of application, producing not only native-like sequences but also thermodynamic properties with errors that appear comparable to other current software packages.

A computational protein design protocol for optimization of the SARS-CoV-2 receptor-binding-motif affinity for human ACE2.

The present protocol describes the computational design of the SARS-CoV-2 receptor binding motif (RBD) to identify mutations that can potentially improve binding affinity for the human ACE2 (hACE2) receptor. We focus on four positions located at the interface with the hACE2 receptor in the RBD:hACE2 complex. We conduct the design with a high-throughput computational protein design (CPD) program, Proteus, incorporating an adaptive Monte Carlo (MC) protocol that promotes the selection of sequences with good binding affinities. For complete details on the use and execution of this protocol, please refer to Polydorides and Archontis (2021).

Computational Design of Peptides with Improved Recognition of the Focal Adhesion Kinase FAT Domain.

We describe a two-stage computational protein design (CPD) methodology for the design of peptides binding to the FAT domain of the protein focal adhesion kinase. The first stage involves high-throughput CPD calculations with the Proteus software. The energies of the folded state are described by a physics-based energy function and of the unfolded peptides by a knowledge-based model that reproduces aminoacid compositions consistent with a helicity scale. The obtained sequences are filtered in terms of the affinity and the stability of the complex. In the second stage, design sequences are further evaluated by all-atom molecular dynamics simulations and binding free energy calculations with a molecular mechanics/implicit solvent free energy function.

Computational protein design with a generalized Born solvent model: application to Asparaginyl-tRNA synthetase.

Computational Protein Design (CPD) is a promising method for high throughput protein and ligand mutagenesis. Recently, we developed a CPD method that used a polar-hydrogen energy function for protein interactions and a Coulomb/Accessible Surface Area (CASA) model for solvent effects. We applied this method to engineer aspartyl-adenylate (AspAMP) specificity into Asparaginyl-tRNA synthetase (AsnRS), whose substrate is asparaginyl-adenylate (AsnAMP). Here, we implement a more accurate function, with an all-atom energy for protein interactions and a residue-pairwise generalized Born model for solvent effects. As a first test, we compute aminoacid affinities for several point mutants of Aspartyl-tRNA synthetase (AspRS) and Tyrosyl-tRNA synthetase and stability changes for three helical peptides and compare with experiment. As a second test, we readdress the problem of AsnRS aminoacid engineering. We compare three design criteria, which optimize the folding free-energy, the absolute AspAMP affinity, and the relative (AspAMP-AsnAMP) affinity. The sequences and conformations are improved with respect to our previous, polar-hydrogen/CASA study: For several designed complexes, the AspAMP carboxylate forms three interactions with a conserved arginine and a designed lysine, as in the active site of the AspRS:AspAMP complex. The conformations and interactions are well maintained in molecular dynamics simulations and the sequences have an inverted specificity, favoring AspAMP over AsnAMP. The method is not fully successful, since experimental measurements with the seven most promising sequences show that they do not catalyze at a detectable level the adenylation of Asp (or Asn) with ATP. This may be due to weak AspAMP binding and/or disruption of transition-state stabilization.

Backbone charge transport in double-stranded DNA.

Understanding charge transport in DNA molecules is a long-standing problem of fundamental importance across disciplines1,2. It is also of great technological interest due to DNA's ability to form versatile and complex programmable structures. Charge transport in DNA-based junctions has been reported using a wide variety of set-ups2-4, but experiments so far have yielded seemingly contradictory results that range from insulating5-8 or semiconducting9,10 to metallic-like behaviour11. As a result, the intrinsic charge transport mechanism in molecular junction set-ups is not well understood, which is mainly due to the lack of techniques to form reproducible and stable contacts with individual long DNA molecules. Here we report charge-transport measurements through single 30-nm-long double-stranded DNA (dsDNA) molecules with an experimental set-up that enables us to address individual molecules repeatedly and to measure the current-voltage characteristics from 5 K up to room temperature. Strikingly, we observed very high currents of tens of nanoamperes, which flowed through both homogeneous and non-homogeneous base-pair sequences. The currents are fairly temperature independent in the range 5-60 K and show a power-law decrease with temperature above 60 K, which is reminiscent of charge transport in organic crystals. Moreover, we show that the presence of even a single discontinuity ('nick') in both strands that compose the dsDNA leads to complete suppression of the current, which suggests that the backbones mediate the long-distance conduction in dsDNA, contrary to the common wisdom in DNA electronics2-4.

Proteus and the Design of Ligand Binding Sites.

This chapter describes the organization and use of Proteus, a multitool computational suite for the optimization of protein and ligand conformations and sequences, and the calculation of pK α shifts and relative binding affinities. The software offers the use of several molecular mechanics force fields and solvent models, including two generalized Born variants, and a large range of scoring functions, which can combine protein stability, ligand affinity, and ligand specificity terms, for positive and negative design. We present in detail the steps for structure preparation, system setup, construction of the interaction energy matrix, protein sequence and structure optimizations, pK α calculations, and ligand titration calculations. We discuss illustrative examples, including the chemical/structural optimization of a complex between the MHC class II protein HLA-DQ8 and the vinculin epitope, and the chemical optimization of the compstatin analog Ac-Val4Trp/His9Ala, which regulates the function of protein C3 of the complement system.

Crossreactivity to vinculin and microbes provides a molecular basis for HLA-based protection against rheumatoid arthritis.

The HLA locus is the strongest risk factor for anti-citrullinated protein antibody (ACPA)(+) rheumatoid arthritis (RA). Despite considerable efforts in the last 35 years, this association is poorly understood. Here we identify (citrullinated) vinculin, present in the joints of ACPA(+) RA patients, as an autoantigen targeted by ACPA and CD4(+) T cells. These T cells recognize an epitope with the core sequence DERAA, which is also found in many microbes and in protective HLA-DRB1*13 molecules, presented by predisposing HLA-DQ molecules. Moreover, these T cells crossreact with vinculin-derived and microbial-derived DERAA epitopes. Intriguingly, DERAA-directed T cells are not detected in HLA-DRB1*13(+) donors, indicating that the DERAA epitope from HLA-DRB1*13 mediates (thymic) tolerance in these donors and explaining the protective effects associated with HLA-DRB1*13. Together our data indicate the involvement of pathogen-induced DERAA-directed T cells in the HLA-RA association and provide a molecular basis for the contribution of protective/predisposing HLA alleles.

Predicting the acid/base behavior of proteins: a constant-pH Monte Carlo approach with generalized born solvent.

The acid/base properties of proteins are essential in biochemistry, and proton binding is a valuable reporter on electrostatic interactions. We propose a constant-pH Monte Carlo strategy to compute protonation free energies and pK(a)'s. The solvent is described implicitly, through a generalized Born model. The electronic polarizability and backbone motions of the protein are included through the protein dielectric constant. Side chain motions are described explicitly, by the Monte Carlo scheme. An efficient computational algorithm is described, which allows us to treat the fluctuating shape of the protein/solvent boundary in a way that is numerically exact (within the GB framework); this contrasts with several previous constant-pH approaches. For a test set of six proteins and 78 titratable groups, the model performs well, with an rms error of 1.2 pH units. While this is slightly greater than a simple Null model (rms error of 1.1) and a fully empirical model (rms error of 0.9), it is obtained using physically meaningful model parameters, including a low protein dielectric of four. Importantly, similar performance is obtained for side chains with large and small pK(a) shifts (relative to a standard model compound). The titration curve slopes and the conformations sampled are reasonable. Several directions to improve the method further are discussed.