A new Lamarckian genetic algorithm for flexible ligand‐receptor docking

We present a Lamarckian genetic algorithm (LGA) variant for flexible ligand‐receptor docking which allows to handle a large number of degrees of freedom. Our hybrid method combines a multi‐deme LGA with a recently published gradient‐based method for local optimization of molecular complexes. We compared the performance of our new hybrid method to two non gradient‐based search heuristics on the Astex diverse set for flexible ligand‐receptor docking. Our results show that the novel approach is clearly superior to other LGAs employing a stochastic optimization method. The new algorithm features a shorter run time and gives substantially better results, especially with increasing complexity of the ligands. Thus, it may be used to dock ligands with many rotatable bonds with high efficiency. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

A consensus line search algorithm for molecular potential energy functions

Force field based energy minimization of molecular structures is a central task in computational chemistry and biology. Solving this problem usually requires efficient local minimization techniques, i.e., iterative two‐step methods that search first for a descent direction and then try to estimate the step width. The second step, the so called line search, typically uses polynomial interpolation schemes to estimate the next trial step. However, dependent on local properties of the objective function alternative schemes may be more appropriate especially if the objective function shows singularities or exponential behavior. As the choice of the best interpolation scheme cannot be made a priori, we propose a new consensus line search approach that performs several different interpolation schemes at each step and then decides which one is the most reliable at the current position. Although a naive consensus approach would lead to severe performance impacts, our method does not require additional evaluations of the energy function, imposing only negligible computational overhead. Additionally, our method can be easily adapted to the local behavior of other objective functions by incorporating suitable interpolation schemes or omitting non‐fitting schemes. The performance of our consensus line search approach has been evaluated and compared to established standard line search algorithms by minimizing the structures of a large set of molecules using different force fields. The proposed algorithm shows better performance in almost all test cases, i.e., it reduces the number of iterations and function and gradient evaluations, leading to significantly reduced run times. © 2008 Wiley Periodicals, Inc. J Comput Chem, 2009     

Global optimization using bad derivatives: Derivative-free method for molecular energy minimization

A general method designed to isolate the global minimum of a multidimensional objective function with multiple minima is presented. The algorithm exploits an integral “coarse-graining” transformation of the objective function, U, into a smoothed function with few minima. When the coarse-graining is defined over a cubic neighborhood of length scale ϵ, the exact gradient of the smoothed function, 𝒰_ϵ, is a simple three-point finite difference of U. When ϵ is very large, the gradient of 𝒰_ϵ appears to be a “bad derivative” of U. Because the gradient of 𝒰_ϵ is a simple function of U, minimization on the smoothed surface requires no explicit calculation or differentiation of 𝒰_ϵ. The minimization method is “derivative-free” and may be applied to optimization problems involving functions that are not smooth or differentiable. Generalization to functions in high-dimensional space is straightforward. In the context of molecular conformational optimization, the method may be used to minimize the potential energy or, preferably, to maximize the Boltzmann probability function. The algorithm is applied to conformational optimization of a model potential, Lennard–Jones atomic clusters, and a tetrapeptide. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1445–1455, 1998     

Conformational search using a molecular dynamics–minimization procedure: Applications to clusters of coulombic charges,Lennard–Jones particles,and waters

A new conformational search method, molecular dynamics–minimization (MDM), is proposed, which combines a molecular dynamics sampling strategy with energy minimizations in the search for low-energy molecular structures. This new method is applied to the search for low energy configurations of clusters of coulombic charges on a unit sphere, Lennard–Jones clusters, and water clusters. The MDM method is shown to be efficient in finding the lowest energy configurations of these clusters. A closer comparison of MDM with alternative conformational search methods on Lennard–Jones clusters shows that, although MDM is not as efficient as the Monte Carlo–minimization method in locating the global energy minima, it is more efficient than the diffusion equation method and the method of minimization from randomly generated structures. Given the versatility of the molecular dynamics sampling strategy in comparison to Monte Carlo in treating molecular complexes or molecules in explicit solution, one anticipates that the MDM method could be profitably applied to conformational search problems where the number of degrees of freedom is much greater. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 60–70, 1998     

Exploring hierarchical refinement techniques for induced fit docking with protein and ligand flexibility

We present a series of molecular‐mechanics‐based protein refinement methods, including two novel ones, applied as part of an induced fit docking procedure. The methods used include minimization; protein and ligand sidechain prediction; a hierarchical ligand placement procedure similar to a‐priori protein loop predictions; and a minimized Monte Carlo approach using normal mode analysis as a move step. The results clearly indicate the importance of a proper opening of the active site backbone, which might not be accomplished when the ligand degrees of freedom are prioritized. The most accurate method consisted of the minimized Monte Carlo procedure designed to open the active site followed by a hierarchical optimization of the sidechain packing around a mobile flexible ligand. The methods have been used on a series of 88 protein‐ligand complexes including both cross‐docking and apo‐docking members resulting in complex conformations determined to within 2.0 Å heavy‐atom RMSD in 75% of cases where the protein backbone rearrangement upon binding is less than 1.0 Å α‐carbon RMSD. We also demonstrate that physics‐based all‐atom potentials can be more accurate than docking‐style potentials when complexes are sufficiently refined. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Geometry optimization of molecular clusters and complexes using scaled internal coordinates

Scaled internal coordinates are introduced for use in the geometry optimization of systems composed of multiple fragments, such as solvated molecules, clusters, and biomolecular complexes. The new coordinates are related to bond lengths, bond angles and torsion angles by geometry-dependent scaling factors. The scaling factors serve to expedite the optimization of complexes containing outlying fragments, without hindering the optimization of the intramolecular degrees of freedom. Trial calculations indicate that, at asymptotic separations, the scaling factors improve the rate of convergence by a factor of 4 to 5.     

An analysis of the structural and energetic properties of deoxyribose by potential energy methods

We discuss the three fundamental issues of a computational approach in structure prediction by potential energy minimization, and analyze them for the nucleic acid component deoxyribose. Predicting the conformation of deoxyribose is important not only because of the molecule's central conformational role in the nucleotide backbone, but also because energetic and geometric discrepancies from experimental data have exposed some underlying uncertainties in potential energy calculations. The three fundamental issues examined here are: (i) choice of coordinate system to represent the molecular conformation; (ii) construction of the potential energy function; and (iii) choice of the minimization technique. For our study, we use the following combination. First, the molecular conformation is represented in cartesian coordinate space with the full set of degrees of freedom. This provides an opportunity for comparison with the pseudorotation approximation. Second, the potential energy function is constructed so that all the interactions other than the nonbonded terms are represented by polynomials of the coordinate variables. Third, two powerful Newton methods that are globally and quadratically convergent are implemented: Gill and Murray's Modified Newton method and a Truncated Newton method, specifically developed for potential energy minimization. These strategies have produced the two experimentally-observed structures of deoxyribose with geometric data (bond angles and dihedral angles) in very good agreement with experiment. More generally, the application of these modeling and minimization techniques to potential energy investigations is promising. The use of cartesian variables and polynomial representation of bond length, bond angle and torsional potentials promotes efficient second-derivative computation and, hence, application of Newton methods. The truncated Newton, in particular, is ideally suited for potential energy minimization not only because the storage and computational requirements of Newton methods are made manageable, but also because it contains an important algorithmic adaptive feature: the minimization search is diverted from regions where the function is nonconvex and is directed quickly toward physically interesting regions.     

Conformational search of peptides and proteins: Monte Carlo minimization with an adaptive bias method applied to the heptapeptide deltorphin

The energy function of a protein consists of a tremendous number of minima. Locating the global energy minimum (GEM) structure, which corresponds approximately to the native structure, is a severe problem in global optimization. Recently we have proposed a conformational search technique based on the Monte Carlo minimization (MCM) method of Li and Scheraga, where trial dihedral angles are not selected at random within the range [-180 degrees,180 degrees ] (as with MCM) but with biased probabilities depending on the increased structure-energy correlations as the GEM is approached during the search. This method, called the Monte Carlo minimization with an adaptive bias (MCMAB), was applied initially to the pentapeptide Leu-enkephalin. Here we study its properties further by applying it to the larger peptide with bulky side chains, deltorphin (H-Tyr-D-Met-Phe-His-Leu-Met-Asp-NH(2)). We find that on average the number of energy minimizations required by MCMAB to locate the GEM for the first time is smaller by a factor of approximately three than the number required by MCM-in accord with results obtained for Leu-enkephalin.     

Self-healing umbrella sampling: a non-equilibrium approach for quantitative free energy calculations

We propose a new approach for the umbrella sampling method in molecular dynamics simulations of complex systems. An accelerated sampling of the slow degrees of freedom is achieved by generating a single self-adaptive trajectory that tends to span uniformly the reaction coordinate using a time dependent bias potential derived from the preceding history of the system. To show the convergent behavior and the efficiency of the method, we present the free energy surface of alanine dipeptide in water as a function of the backbone dihedral angles.     

Direct minimization of the energy functional in the LCAO-MO density matrix formalism

A previously proposed method of energy minimization is developed for MC SCF wavefunctions formed by all-pair excitations for a closed-shell system. The orbital coefficients are optimized by a gradient approach using a suitable orthogonal transformation of the atomic basis, while optimum CI coefficients are determined solving the usual secular problem for the lowest eigenvalue, after each optimization of the orbitals. Applications to LiH and NH₃ molecules show that the method is numerically well stable, and is capable of accounting for a large part of the correlation energy giving results which compare well with those of the conventional CI method.     

Prodock: Software package for protein modeling and docking

A new software package, Prodock , for protein modeling and flexible docking is presented. The protein system is described in internal coordinates with an arbitrary level of flexiblity for the proteins or ligands. The protein is represented by an all-atom model with the Ecepp /3 or Amber IV force field, depending on whether the ligand is a peptidic molecule or not. Prodock is based on a new residue data dictionary that makes the programming easier and the definition of molecular flexibility more straigthforward. Two versions of the dictionary have been constructed for the Ecepp /3 and Amber IV geometry, respectively. The global optimization of the energy function is carried out with the scaled collective variable Monte Carlo method plus energy minimization. The incorporation of a local minimization during the conformational sampling has been shown to be very important for distinguishing low-energy nonnative conformations from native structures. To make the Monte Carlo minimization method efficient for docking, a new grid-based energy evaluation technique using Bezier splines has been incorporated. This article includes some techniques and simulation tools that significantly improve the efficiency of flexible docking simulations, in particular forward/backward polypeptide chain generation. A comparative study to illustrate the advantage of using quaternions over Euler angles for the rigid-body rotational variables is presented in this paper. Several applications of the program Prodock are also discussed. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 412–427, 1999     

Optimization and dynamics of protein-protein complexes using B-splines

A moving-grid approach for optimization and dynamics of protein-protein complexes is introduced, which utilizes cubic B-spline interpolation for rapid energy and force evaluation. The method allows for the efficient use of full electrostatic potentials joined smoothly to multipoles at long distance so that multiprotein simulation is possible. Using a recently published benchmark of 58 protein complexes, we examine the performance and quality of the grid approximation, refining cocrystallized complexes to within 0.68 A RMSD of interface atoms, close to the optimum 0.63 A produced by the underlying MMFF94 force field. We quantify the theoretical statistical advantage of using minimization in a stochastic search in the case of two rigid bodies, and contrast it with the underlying cost of conjugate gradient minimization using B-splines. The volumes of conjugate gradient minimization basins of attraction in cocrystallized systems are generally orders of magnitude larger than well volumes based on energy thresholds needed to discriminate native from nonnative states; nonetheless, computational cost is significant. Molecular dynamics using B-splines is doubly efficient due to the combined advantages of rapid force evaluation and large simulation step sizes. Large basins localized around the native state and other possible binding sites are identifiable during simulations of protein-protein motion. In addition to providing increased modeling detail, B-splines offer new algorithmic possibilities that should be valuable in refining docking candidates and studying global complex behavior.     

Radial‐template approach for accurate density representation in computational quantum theory

A method for improved representation of electronic charge and spin densities for molecular and solid state systems is presented, based upon extensions of least squares fits to quantum mechanical “true” densities using basis functions of limited support. Attention is given to optimization of radial degrees of fit freedom, and the design of fit functions permitting rapid analytic manipulation and calculation of properties, such as Coulomb potentials. The method is demonstrated for covalent CO and for a large metal‐organic crystalline structure. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

An efficient molecular docking using conformational space annealing

Molecular docking falls into the general category of global optimization problems because its main purpose is to find the most stable complex consisting of a receptor and its ligand. Conformational space annealing (CSA), a powerful global optimization method, is incorporated with the Tinker molecular modeling package to perform molecular docking simulations of six receptor-ligand complexes (3PTB, 1ULB, 2CPP, 1STP, 3CPA, and 1PPH) from the Protein Data Bank. In parallel, Monte Carlo with the minimization (MCM) method is also incorporated into the Tinker package for comparison. The energy function, consisting of electrostatic interactions, van der Waals interactions, and torsional energy terms, is calculated using the AMBER94 all-atom empirical force field. Rigid docking simulations for all six complexes and flexible docking simulations for three complexes (1STP, 3CPA, and 1PPH) are carried out using the CSA and the MCM methods. The simulation results show that the docking procedures using the CSA method generally find the most stable complexes as well as the native-like complexes more efficiently and accurately than those using the MCM, demonstrating that CSA is a promising search method for molecular docking problems.     

Direct minimization of the energy functional in LCAO-MO density matrix formalism

A general theory is presented for the optimization of the coefficients of orbitals and configuration interaction expansion in the case of multiconfiguration wavefunctions containing all single excitations. The orbital coefficients are optimized by suitable orthogonal transformations of the atomic basis; the Cl coefficients are determined solving the usual secular problem. The energy minimization is performed directly by a gradient approach. The method works both for ground and excited states and no convergence difficulties are met. Computational examples are given for H₂O and H₂S molecules.     

Conical intersections in solution: formulation, algorithm, and implementation with combined quantum mechanics/molecular mechanics method

The significance of conical intersections in photophysics, photochemistry, and photodissociation of polyatomic molecules in gas phase has been demonstrated by numerous experimental and theoretical studies. Optimization of conical intersections of small- and medium-size molecules in gas phase has currently become a routine optimization process, as it has been implemented in many electronic structure packages. However, optimization of conical intersections of small- and medium-size molecules in solution or macromolecules remains inefficient, even poorly defined, due to large number of degrees of freedom and costly evaluations of gradient difference and nonadiabatic coupling vectors. In this work, based on the sequential quantum mechanics and molecular mechanics (QM/MM) and QM/MM-minimum free energy path methods, we have designed two conical intersection optimization methods for small- and medium-size molecules in solution or macromolecules. The first one is sequential QM conical intersection optimization and MM minimization for potential energy surfaces; the second one is sequential QM conical intersection optimization and MM sampling for potential of mean force surfaces, i.e., free energy surfaces. In such methods, the region where electronic structures change remarkably is placed into the QM subsystem, while the rest of the system is placed into the MM subsystem; thus, dimensionalities of gradient difference and nonadiabatic coupling vectors are decreased due to the relatively small QM subsystem. Furthermore, in comparison with the concurrent optimization scheme, sequential QM conical intersection optimization and MM minimization or sampling reduce the number of evaluations of gradient difference and nonadiabatic coupling vectors because these vectors need to be calculated only when the QM subsystem moves, independent of the MM minimization or sampling. Taken together, costly evaluations of gradient difference and nonadiabatic coupling vectors in solution or macromolecules can be reduced significantly. Test optimizations of conical intersections of cyclopropanone and acetaldehyde in aqueous solution have been carried out successfully.     

Controlling an Electron‐Transfer Reaction at a Metal Surface by Manipulating Reactant Motion and Orientation

The loss or gain of vibrational energy in collisions of an NO molecule with the surface of a gold single crystal proceeds by electron transfer. With the advent of new optical pumping and orientation methods, we can now control all molecular degrees of freedom important to this electron‐transfer‐mediated process, providing the most detailed look yet into the inner workings of an electron‐transfer reaction and showing how to control its outcome. We find the probability of electron transfer increases with increasing translational and vibrational energy as well as with proper orientation of the reactant. However, as the vibrational energy increases, translational excitation becomes unimportant and proper orientation becomes less critical. One can understand the interplay of all three control parameters from simple model potentials.     

Controlled-advancement rigid-body optimization of nanosystems

In this study, we propose a novel optimization algorithm, with application to the refinement of molecular complexes. Particularly, we consider optimization problem as the calculation of quasi-static trajectories of rigid bodies influenced by the inverse-inertia-weighted energy gradient and introduce the concept of advancement region that guarantees displacement of a molecule strictly within a relevant region of conformational space. The advancement region helps to avoid typical energy minimization pitfalls, thus, the algorithm is suitable to work with arbitrary energy functions and arbitrary types of molecular complexes without necessary tuning of its hyper-parameters. Our method, called controlled-advancement rigid-body optimization of nanosystems (Carbon), is particularly useful for the large-scale molecular refinement, as for example, the putative binding candidates obtained with protein–protein docking pipelines. Implementation of Carbon with user-friendly interface is available in the SAMSON platform for molecular modeling at https://www.samson-connect.net . © 2019 Wiley Periodicals, Inc.     

A kinematic view of loop closure

We consider the problem of loop closure, i.e., of finding the ensemble of possible backbone structures of a chain segment of a protein molecule that is geometrically consistent with preceding and following parts of the chain whose structures are given. We reduce this problem of determining the loop conformations of six torsions to finding the real roots of a 16th degree polynomial in one variable, based on the robotics literature on the kinematics of the equivalent rotator linkage in the most general case of oblique rotators. We provide a simple intuitive view and derivation of the polynomial for the case in which each of the three pair of torsional axes has a common point. Our method generalizes previous work on analytical loop closure in that the torsion angles need not be consecutive, and any rigid intervening segments are allowed between the free torsions. Our approach also allows for a small degree of flexibility in the bond angles and the peptide torsion angles; this substantially enlarges the space of solvable configurations as is demonstrated by an application of the method to the modeling of cyclic pentapeptides. We give further applications to two important problems. First, we show that this analytical loop closure algorithm can be efficiently combined with an existing loop-construction algorithm to sample loops longer than three residues. Second, we show that Monte Carlo minimization is made severalfold more efficient by employing the local moves generated by the loop closure algorithm, when applied to the global minimization of an eight-residue loop. Our loop closure algorithm is freely available at http://dillgroup. ucsf.edu/loop_closure/.