Accurate ligand-protein docking in CASP15 using the ClusPro LigTBM server.

In the ligand prediction category of CASP15, the challenge was to predict the positions and conformations of small molecules binding to proteins that were provided as amino acid sequences or as models generated by the AlphaFold2 program. For most targets, we used our template-based ligand docking program ClusPro ligTBM, also implemented as a public server available at https://ligtbm.cluspro.org/. Since many targets had multiple chains and a number of ligands, several templates, and some manual interventions were required. In a few cases, no templates were found, and we had to use direct docking using the Glide program. Nevertheless, ligTBM was shown to be a very useful tool, and by any ranking criteria, our group was ranked among the top five best-performing teams. In fact, all the best groups used template-based docking methods. Thus, it appears that the AlphaFold2-generated models, despite the high accuracy of the predicted backbone, have local differences from the x-ray structure that make the use of direct docking methods more challenging. The results of CASP15 confirm that this limitation can be frequently overcome by homology-based docking.

Accurate de novo design of hyperstable constrained peptides.

Naturally occurring, pharmacologically active peptides constrained with covalent crosslinks generally have shapes that have evolved to fit precisely into binding pockets on their targets. Such peptides can have excellent pharmaceutical properties, combining the stability and tissue penetration of small-molecule drugs with the specificity of much larger protein therapeutics. The ability to design constrained peptides with precisely specified tertiary structures would enable the design of shape-complementary inhibitors of arbitrary targets. Here we describe the development of computational methods for accurate de novo design of conformationally restricted peptides, and the use of these methods to design 18-47 residue, disulfide-crosslinked peptides, a subset of which are heterochiral and/or N-C backbone-cyclized. Both genetically encodable and non-canonical peptides are exceptionally stable to thermal and chemical denaturation, and 12 experimentally determined X-ray and NMR structures are nearly identical to the computational design models. The computational design methods and stable scaffolds presented here provide the basis for development of a new generation of peptide-based drugs.

BioMThermDB 1.0: Thermophysical Database of Proteins in Solutions.

We present here a freely available web-based database, called BioMThermDB 1.0, of thermophysical and dynamic properties of various proteins and their aqueous solutions. It contains the hydrodynamic radius, electrophoretic mobility, zeta potential, self-diffusion coefficient, solution viscosity, and cloud-point temperature, as well as the conditions for those determinations and details of the experimental method. It can facilitate the meta-analysis and visualization of data, can enable comparisons, and may be useful for comparing theoretical model predictions with experiments.

Protein secondary structure motifs: A kinematic construction.

The kinematic geometry of protein backbone structures, constrained by either single or multiple hydrogen bonds (H-bonds), possibly in a periodic array, is discussed. These structures include regular secondary structure elements α-helices and ß-sheets but also include other short H-bond stabilized irregular structural elements like ß-turns. The work here shows that the variations observed in such structures have simple geometrical correlations consistent with constrained motion kinematics. A new classification of the ideal helices is given, in terms of the parameter α, the angle at a Cα atom to its two neighboring Cα 's along the helix, and shown how it can be generalized to include nonideal helices. Specifically, we derive an analytical expression of the backbone dihedrals, (ϕ, ψ), in terms of the parameter α subject to the constraint that the peptide planes are parallel to the helical axis. Helices constructed in this way exhibit near-vertical alignment of the C = O and N - H units and are the canonical objects of this study. These expressions are easily modifiable to include perturbations of parameters relevant to nonplanar peptide units and noncanonical angles. The addition of a second parameter, ε0 , inclination of successive peptide planes along a helix with respect to the helical axis leads to a generalization of the previous expression and provides an efficient parametrization of such structures in terms of coordinates consistent with H-bond parameters. An analogs parametrization of ß-turns, using inverse kinematic methods, is also given. Besides offering a unifying viewpoint, our results may find useful applications to protein and peptide design.

Kinematic Reconstruction of Cyclic Peptides and Protein Backbones from Partial Data.

We present an algorithm, QBKR (Quaternary Backbone Kinematic Reconstruction), a fast analytical method for an all-atom backbone reconstruction of proteins and linear or cyclic peptide chains from Cα coordinate traces. Unlike previous analytical methods for deriving all-atom representations from coarse-grained models that rely on canonical geometry with planar peptides in the trans conformation, our de novo kinematic model incorporates noncanonical, cis-trans, geometry naturally. Perturbations to this geometry can be effected with ease in our formulation, for example, to account for a continuous change from cis to trans geometry. A simple optimization of a spring-based objective function is employed for Cα-Cα distance variations that extend beyond the cis-trans limit. The kinematic construction produces a linked chain of peptide units, Cα-C-N-Cα, hinged at the Cα atoms spanning all possible planar and nonplanar peptide conformations. We have combined our method with a ring closure algorithm for the case of ring peptides and missing loops in a protein structure. Here, the reconstruction proceeding from both the N and C termini of the protein backbone (or in both directions from a starting position for rings) requires freedom in the position of one Cα atom (a capstone) to achieve a successful loop or ring closure. A salient feature of our reconstruction method is the ability to enrich conformational ensembles to produce alternative feasible conformations in which H-bond forming C-O or N-H pairs in the backbone can reverse orientations, thus addressing a well-known shortcoming in Cα-based RMSD structure comparison, wherein very close structures may lead to significantly different overall H-bond behavior. We apply the fixed Cα-based design to the reverse reconstruction from noisy Cryo-EM data, a posteriori to the optimization. Our method can be applied to speed up the process of an all-atom description from voluminous experimental data or subpar electron density maps.

Sampling and refinement protocols for template-based macrocycle docking: 2018 D3R Grand Challenge 4.

We describe a new template-based method for docking flexible ligands such as macrocycles to proteins. It combines Monte-Carlo energy minimization on the manifold, a fast manifold search method, with BRIKARD for complex flexible ligand searching, and with the MELD accelerator of Replica-Exchange Molecular Dynamics simulations for atomistic degrees of freedom. Here we test the method in the Drug Design Data Resource blind Grand Challenge competition. This method was among the best performers in the competition, giving sub-angstrom prediction quality for the majority of the targets.

RMSD and Symmetry.

A common approach for comparing the structures of biomolecules or solid bodies is to translate and rotate one structure with respect to the other to minimize the pointwise root-mean-square deviation (RMSD). We present a new, robust numerical algorithm that computes the RMSD between two molecules or all the mutual RMSDs of a list of molecules and, if desired, the corresponding rotation matrix in a minimal number of operations as compared to previous algorithms. The RMSD gradient can also be computed. We address the problem of symmetry, both in alignment (possible alternative alignments due to indistinguishable atoms) as well as geometry. In the latter case, it is possible to have degenerate superposition. A necessary condition is optimal superimposability to one's mirror image. Double (respectively, triple) degeneracy results in a one- (respectively, two)-parameter family of rotations leaving the superposition invariant. The software, frmsd, is freely available at http://www.ams.stonybrook.edu/~coutsias/codes/frmsd.tgz. © 2019 Wiley Periodicals, Inc.

A human transcription factor in search mode.

Transcription factors (TF) can change shape to bind and recognize DNA, shifting the energy landscape from a weak binding, rapid search mode to a higher affinity recognition mode. However, the mechanism(s) driving this conformational change remains unresolved and in most cases high-resolution structures of the non-specific complexes are unavailable. Here, we investigate the conformational switch of the human mitochondrial transcription termination factor MTERF1, which has a modular, superhelical topology complementary to DNA. Our goal was to characterize the details of the non-specific search mode to complement the crystal structure of the specific binding complex, providing a basis for understanding the recognition mechanism. In the specific complex, MTERF1 binds a significantly distorted and unwound DNA structure, exhibiting a protein conformation incompatible with binding to B-form DNA. In contrast, our simulations of apo MTERF1 revealed significant flexibility, sampling structures with superhelical pitch and radius complementary to the major groove of B-DNA. Docking these structures to B-DNA followed by unrestrained MD simulations led to a stable complex in which MTERF1 was observed to undergo spontaneous diffusion on the DNA. Overall, the data support an MTERF1-DNA binding and recognition mechanism driven by intrinsic dynamics of the MTERF1 superhelical topology.

LAGUERRE-INTERSECTION METHOD FOR IMPLICIT SOLVATION.

Explicit solvent molecular dynamics simulations of a macromolecule are slow as the number of solvent atoms considered typically increases by order of magnitude. Implicit methods introduce surface-dependent corrections to the force field, gaining speed at the expense of accuracy. Properties such as molecular interface surfaces, volumes and cavities are captured by Laguerre tessellations of macromolecules. However, Laguerre cells of exterior atoms tend to be overly large or unbounded. Our method, the inclusion-exclusion based Laguerre-Intersection method, caps cells in a physically accurate manner by considering the intersection of the space-filling diagram with the Laguerre tessellation. We optimize an adjustable parameter, the weight, to ensure the areas and volumes of capped cells exposed to solvent are as close as possible, on average, to those computed from equilibrated explicit solvent simulations. The contact planes are radical planes, meaning that as the solvent weight is varied, interior cells remain constant. We test the consistency of our model using a high-quality trajectory of HIV-protease, a dimer with flexible loops and open-close transitions. We also compare our results with interval-arithmetic Gauss-Bonnet based method. Optimal solvent parameters quickly converge, which we use to illustrate the increased fidelity of the Laguerre-Intersection method over two recently proposed methods as compared to the explicit model.

Characterization of Biomolecular Helices and Their Complementarity Using Geometric Analysis.

A general method is presented to characterize the helical properties of potentially irregular helices, such as those found in protein secondary and tertiary structures and nucleic acids. The method was validated using artificial helices with varying numbers of points, points per helical turn, pitch, and radius. The sensitivity of the method was validated by applying increasing amounts of random perturbation to the coordinates of these helices; 399 360 helices in total were evaluated. In addition, the helical parameters of protein secondary structure elements and nucleic acid helices were analyzed. Generally, at least seven points were required to recapitulate the parameters of a helix using our method. The method can also be used to calculate the helical parameters of nucleic acid-binding proteins, like TALE, enabling direct analysis of their helix complementarity to sequence-dependent DNA distortions.

Flexibility of Bricard's linkages and other structures via resultants and computer algebra.

Flexibility of structures is extremely important for chemistry and robotics. Following our earlier work, we study flexibility using polynomial equations, resultants, and a symbolic algorithm of our creation that analyzes the resultant. We show that the software solves a classic arrangement of quadrilaterals in the plane due to Bricard. We fill in several gaps in Bricard's work and discover new flexible arrangements that he was apparently unaware of. This provides strong evidence for the maturity of the software, and is a wonderful example of mathematical discovery via computer assisted experiment.

Constraint methods that accelerate free-energy simulations of biomolecules.

Atomistic molecular dynamics simulations of biomolecules are critical for generating narratives about biological mechanisms. The power of atomistic simulations is that these are physics-based methods that satisfy Boltzmann's law, so they can be used to compute populations, dynamics, and mechanisms. But physical simulations are computationally intensive and do not scale well to the sizes of many important biomolecules. One way to speed up physical simulations is by coarse-graining the potential function. Another way is to harness structural knowledge, often by imposing spring-like restraints. But harnessing external knowledge in physical simulations is problematic because knowledge, data, or hunches have errors, noise, and combinatoric uncertainties. Here, we review recent principled methods for imposing restraints to speed up physics-based molecular simulations that promise to scale to larger biomolecules and motions.

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

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

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

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

The FALC-Loop web server for protein loop modeling.

The FALC-Loop web server provides an online interface for protein loop modeling by employing an ab initio loop modeling method called FALC (fragment assembly and analytical loop closure). The server may be used to construct loop regions in homology modeling, to refine unreliable loop regions in experimental structures or to model segments of designed sequences. The FALC method is computationally less expensive than typical ab initio methods because the conformational search space is effectively reduced by the use of fragments derived from a structure database. The analytical loop closure algorithm allows efficient search for loop conformations that fit into the protein framework starting from the fragment-assembled structures. The FALC method shows prediction accuracy comparable to other state-of-the-art loop modeling methods. Top-ranked model structures can be visualized on the web server, and an ensemble of loop structures can be downloaded for further analysis. The web server can be freely accessed at http://falc-loop.seoklab.org/.

NUMERICALLY STABLE SOLUTION TO THE 6R PROBLEM OF INVERSE KINEMATICS.

In this paper, a stable and accurate algorithm to compute all solutions of the inverse kinematics problem of a 6 revolute manipulator chain is presented. A system of equations is constructed based on the fundamental closure conditions, leading to a closed algebraic system of 20 equations involving 16 quantities, composed of trigonometric functions of five among the six unknown joint angles. Two among these five are stably eliminated using singular value decomposition (SVD) avoiding the need to consider special cases. The resulting system of equations involving three unknowns is solved by conversion to a generalized eigenvalue problem. The remaining three unknown angles are obtained using the previously computed pseudoinverse. In this formulation we exploit the inherently complex form of the system reducing it to 10 complex equations in 9 quantities, which substantially accelerates the SVD computation. The method's robustness is demonstrated through a comparison to current methods and several examples including known problematic cases where some axis or link lengths vanish, or some joint angles are 180 degrees, as well as cases where multiple eigenvalues arise.

Crustwater: Modeling Hydrophobic Solvation.

We describe Crustwater, a statistical mechanical model of nonpolar solvation in water. It treats bulk water using the Cage Water model and introduces a crust, i.e., a solvation shell of coordinated partially structured waters. Crustwater is analytical and fast to compute. We compute here solvation vs temperature over the liquid range, and vs pressure and solute size. Its thermal predictions are as accurate as much more costly explicit models such as TIP4P/2005. This modeling gives new insights into the hydrophobic effect: (1) that oil-water insolubility in cold water is due to solute-water (SW) translational entropy and not water-water (WW) orientations, even while hot water is dominated by WW cage breaking, and (2) that a size transition at the Angstrom scale, not the nanometer scale, takes place as previously predicted.

Protein loop modeling by using fragment assembly and analytical loop closure.

Protein loops are often involved in important biological functions such as molecular recognition, signal transduction, or enzymatic action. The three dimensional structures of loops can provide essential information for understanding molecular mechanisms behind protein functions. In this article, we develop a novel method for protein loop modeling, where the loop conformations are generated by fragment assembly and analytical loop closure. The fragment assembly method reduces the conformational space drastically, and the analytical loop closure method finds the geometrically consistent loop conformations efficiently. We also derive an analytic formula for the gradient of any analytical function of dihedral angles in the space of closed loops. The gradient can be used to optimize various restraints derived from experiments or databases, for example restraints for preferential interactions between specific residues or for preferred backbone angles. We demonstrate that the current loop modeling method outperforms previous methods that employ residue-based torsion angle maps or different loop closure strategies when tested on two sets of loop targets of lengths ranging from 4 to 12.

Topology of cyclo-octane energy landscape.

Understanding energy landscapes is a major challenge in chemistry and biology. Although a wide variety of methods have been invented and applied to this problem, very little is understood about the actual mathematical structures underlying such landscapes. Perhaps the most general assumption is the idea that energy landscapes are low-dimensional manifolds embedded in high-dimensional Euclidean space. While this is a very mild assumption, we have discovered an example of an energy landscape which is nonmanifold, demonstrating previously unknown mathematical complexity. The example occurs in the energy landscape of cyclo-octane, which was found to have the structure of a reducible algebraic variety, composed of the union of a sphere and a Klein bottle, intersecting in two rings.

Iterative assembly of helical proteins by optimal hydrophobic packing.

We present a method for the computer-based iterative assembly of native-like tertiary structures of helical proteins from alpha-helical fragments. For any pair of helices, our method, called MATCHSTIX, first generates an ensemble of possible relative orientations of the helices with various ways to form hydrophobic contacts between them. Those conformations having steric clashes, or a large radius of gyration of hydrophobic residues, or with helices too far separated to be connected by the intervening linking region, are discarded. Then, we attempt to connect the two helical fragments by using a robotics-based loop-closure algorithm. When loop closure is feasible, the algorithm generates an ensemble of viable interconnecting loops. After energy minimization and clustering, we use a representative set of conformations for further assembly with the remaining helices, adding one helix at a time. To efficiently sample the conformational space, the order of assembly generally proceeds from the pair of helices connected by the shortest loop, followed by joining one of its adjacent helices, always proceeding with the shorter connecting loop. We tested MATCHSTIX on 28 helical proteins each containing up to 5 helices and found it to heavily sample native-like conformations. The average rmsd of the best conformations for the 17 helix-bundle proteins that have 2 or 3 helices is less than 2 A; errors increase somewhat for proteins containing more helices. Native-like states are even more densely sampled when disulfide bonds are known and imposed as restraints. We conclude that, at least for helical proteins, if the secondary structures are known, this rapid rigid-body maximization of hydrophobic interactions can lead to small ensembles of highly native-like structures. It may be useful for protein structure prediction.