"Polarizable force field has been successfully used in molecular modeling for years, especially in biological and protein simulations. In this research thesis, development of a new polarizable force field ―POSSIM (POlarizable Simulations with Second order Interaction Model) involving electrostatic polarization is described and parameters for several protein residues were produced. In this research thesis, the POSSIM force field was extended to the side chains of the following residues: lysine, glutamic acid, prontonated hisidine, phenylalanine and tryptophan. This work involved producing parameters for methyl ammonium, acetate ion, imidazolium cation, benzene and pyrrole molecules. The parameters fitting procedure starts from the molecular complex with dipolar electrostatic probes of a many-body system to produce polarizabilities, compute the energies, then charges and Lennard-Jones parameters are produced by fitting to gas-phase dimerization calculations, followed by the torsional parameters fitting and end up with the pure liquid simulations. In all the cases, three-body energies, dimerization energies and distances agree well to the accurate quantum mechanical results. The final parameters obtained assured the error of less than 2% in the heat of vaporization and average volume results compared with the available experimental data. Unlike the quantum mechanical calculations, the polarizable force field computations require a relatively small amount of computational resources. Moreover, compared to fixed-charges empirical force fields, polarizable force fields are much more accurate in a number of energy calculations. In the following chapters, the results obtained with this particular polarizable force field are discussed."

Accurate prediction of protein-ligand binding affinity is essential to computational drug discovery. Current approaches are limited by the accuracy of the underlying potential energy model that describes atomic interactions. A more rigorous physical model is critical for evaluating molecular interactions to chemical accuracy. The objective of this thesis research is to develop a polarizable force field with an accurate representation of electrostatic interactions, and apply this model to protein-ligand recognition and to ultimately solve practical problems in computer aided drug discovery. By calculating the hydration free energies of a series of organic small molecules, an optimal protocol is established to develop the electrostatic parameters from quantum mechanics calculations. Next, the systematical development and parameterization procedure of AMOEBA protein force field is presented. The derived force field has gone through extensive validations in both gas phase and condensed phase. The last part of the thesis involves the application of AMOEBA to study protein-ligand interactions. The binding free energies of benzamidine analogs to trypsin using molecular dynamics alchemical perturbation are calculated with encouraging accuracy. AMOEBA is also used to study the thermodynamic effect of constraining and hydrophobicity on binding energetics between phosphotyrosine(pY)-containing tripeptides and the SH2 domain of growth receptor binding protein 2 (Grb2). The underlying mechanism of an "entropic paradox" associated with ligand preorganization is explored. / text

Computational chemists have successfully simulated many systems by applying the principles of quantum mechanics, while approximate molecular mechanical models have seen great utility in problems of biochemical interest. In recent years, a number of methods have been developed to combine the advantages of both techniques. In this study the so-called QM/MM method is developed and applied to the determination of the free energy of a simple Menshutkin S_N2 chemical reaction. This is an extremely demanding process, well beyond the computational capacity of an average workstation, and thus a Beowulf-class Linux cluster is constructed to perform the calculations, and tested for a variety of computational chemistry applications. A number of methods for improving the QM/MM approach are considered in this work. The Fluctuating Charge, or FlucQ, polarizable molecular mechanics force field is implemented in a flexible manner within the CHARMM package and tested for a variety of systems, including the S_N2 test case. Several drawbacks of the original method are addressed and overcome. Both molecular dynamics and Monte Carlo techniques are used within the QM/MM framework to investigate the S_N2 reaction, and the two methods are compared. Techniques are developed and tested to increase the efficiency of QM/MC calculations to the point where they become competitive with QM/MD. Extremely expensive QM treatments are shown to be required to obtain accurate energies for the Menshutkin reaction. A method is developed and tested, and compared with the traditional ONIOM technique, for dramatically reducing the computational time required to use these treatments for QM/MC simulations, paving the way for fully ab initio high basis set QM/MM simulation.

Water models are of great importance for different fields of studies such as fluid mechanics, nano materials, and biomolecule simulations. In this dissertation, we focus on the water models applied in atomistic simulations, including those of biomolecules such as proteins and DNA. Despite water's simple structure and countless studies carried out over the decades, the best water models are still far from perfect. Water models are normally divided into two types--explicit model and implicit model. Here my research is mainly focused on explicit models. In explicit water models, fixed charge n-point models are most widely used in atomistic simulations, but have known accuracy drawbacks. Increasing the number of point charges, as well as adding electronic polarizability, are two common strategies for accuracy improvements. Both strategies come at considerable computational cost, which weighs heavily against modest possible accuracy improvements in practical simulations. With a careful comparison between the two strategies, results show that adding polarizability is a more favorable path to take. Optimal point charge approximation (OPCA) method is then applied along with a novel global optimization process, leading to a new polarizable water model OPC3-pol that can reproduce bulk liquid properties of water accurately and run at a speed comparable to 3- and 4-point non-polarizable water models. For practical use, OPC3-pol works with existing non-polarizable AMBER force fields for simulations of globular protein or DNA. In addition, for intrinsically disordered protein simulations, OPC3-pol fixes the over-compactness problem of the previous generation non-polarizable water models. / Doctor of Philosophy / With the rapid advancements of computer technologies, computer simulation has become increasingly popular in biochemistry research fields. Simulations of microscopic substances that are vital for living creatures such as proteins and DNAs have brought us more and more insights into their structures and functions. Because of the fact that almost all the microscopic substances are immersed in water no matter they are in a human body, a plant, or in bacteria, accurately simulating water is crucial for the success of such simulations. My research is focused on developing accurate and fast water models that can be used by researchers in their biochemical simulations. One particular challenge is that water in nature is very flexible and properties of water can change drastically when its surroundings change. Many classical water models cannot correctly mimic this flexibility, and some more advanced water models that are able to mimic it can cost several times more computing resources. Our latest water model OPC3-pol, benefited from a new design, accurately mimics the flexibility and runs as fast as a traditional rigid water model.

Membrane proteins are critical to many biological processes, including molecular transport, signal transduction, and cellular interactions. Through the use of molecular dynamics (MD) simulations, we are able to model this environment at an atomistic scale. However, traditionally used nonpolarizable force fields (FF) are thought to model the unique dielectric gradient posed by the lipid environment with a limited accuracy due to the mean field approximation of charge. Advancements in polarizable FFs and computing efficiency has enabled the explicit modeling of polarization responses and charge distribution, enabling a deeper understanding of the electrostatics driving these processes. Through the use of the Drude FF, we study three specific model systems to understand where explicit polarization is important in describing membranes and membrane proteins. These studies sought to answer the questions: (1) How does explicit electronic polarization impact small molecule permeation and localization preference?, (2) What electrostatic interactions underlie membrane protein secondary structure?, and (3) How do conformational changes propagate between microswitches in G-Protein Coupled Receptors? In this work, we show small molecule dipole moments changing as a function of localization in the bilayer. Additionally, we show differences in the free energy surfaces of permeation for aromatic, polar, and negatively charged species reliant upon force field used. For secondary structure, we showed key interactions which aided to stabilize model helices in bilayers. Finally, we showed potential inductive effects of key microswitch residues underlying prototypical G-Protein coupled receptor activation. This dissertation has helped to show the importance of including explicit polarization in membrane protein systems, especially when considering interactions at the interface and modeling species with charge. This work enables a refined view of the electrostatics occurring in membranes and membrane protein systems, and in the future, can be used as a basis for methodologies in computer aided drug design efforts. / Doctor of Philosophy / Deepening our understandings of membranes and membrane proteins enable better informed and more efficient drug design. In order to do this, biological processes can be simulated through molecular dynamics (MD) simulations. MD simulations use mathematical models known as force fields (FF) to represent the physics of biological systems at an atomistic scale. This enables the study of key interactions which can be leveraged for drug discovery efforts. However, traditional FF neglect electronic structural changes which are crucial for accurately describing the membrane environment and the influence it has on surrounding and embedded molecules. Using enhanced FFs, known as polarizable FFs, we can model this response and gain an entirely new perspective on membranes and membrane proteins. This work helps to define when these FFs are most important to be used when studying membranes and membrane proteins, and in the future, serve as a basis for further simulations in drug discovery efforts.

One of the central tasks for biomolecular modeling is to develop accurate and computationally cheap methods. In this dissertation, we present the development of a brand new polarizable force field—Polarizable Simulations with Second order Interaction Model (POSSIM) involving electrostatic polarization. The POSSIM framework combines accuracy of a polarizable force field and computational efficiency of the second-order approximation of the full-scale induced point dipole polarization formalism. POSSIM force field has been extended to include parameters for small molecules serving as models for peptide and protein side-chains. Parameters have been fitted to permit reproducing many-body energies, gas-phase dimerization energies and geometries and liquid-phase heats of vaporization and densities. Quantum mechanical and experimental data have been used as the target for the fitting. The resulting parameters can be used for simulations of the parameterized molecules themselves or their analogues. In addition to this, these force field parameters have been employed in further development of the POSSIM fast polarizable force field for proteins. The POSSIM framework has been expanded to include a complete polarizable force field for proteins. Most of the parameter fitting was done to high-level quantum mechanical data. Conformational geometries and energies for dipeptides have been reproduced within average errors of ca. 0.5 kcal/mol for energies of the conformers (for the electrostatically neutral residues) and 9.7º for key dihedral angles. We have also validated this force field by simulating an elastin-like polypeptide GVG(VPGVG)3 in aqueous solution. Elastin-like peptides with the (VPGVG)n motif are known to exhibit anomalous behavior of their radius of gyration that increases when temperature is lowered (the so called inverse temperature transition). We have simulated the system with the OPLS-AA and POSSIM force fields and demonstrated that our newly developed polarizable POSSIM parameters permit to capture the experimentally observed decrease of the radius of gyration with increasing temperature, while the fixed-charges OPLS-AA ones do not. Furthermore, our fitting of the force field parameters for the peptides and proteins has been streamlined compared with the previous generation of the complete polarizable force field and relied more on transferability of parameters for non-bonded interactions (including the electrostatic component). The resulting deviations from the quantum mechanical data are similar to those achieved with the previous generation, thus the technique is robust and the parameters are transferable. At the same time, the number of parameters used in this work was noticeably smaller than that of the previous generation of our complete polarizable force field for proteins, thus the transferability of this set can be expected to be greater and the danger of force field fitting artifacts is lower. Therefore, we believe that this force field can be successfully applied in a wide variety of applications to proteins and protein-ligand complexes.

"The accurate, fast and low cost computational tools are indispensable for studying the structure and dynamics of biological macromolecules in aqueous solution. The goal of this thesis is development and validation of continuum Fuzzy-Border (FB) solvation model to work with the Polarizable Simulations Second-order Interaction Model (POSSIM) force field for proteins developed by Professor G A Kaminski. The implicit FB model has advantages over the popularly used Poisson Boltzmann (PB) solvation model. The FB continuum model attenuates the noise and convergence issues commonly present in numerical treatments of the PB model by employing fixed position cubic grid to compute interactions. It also uses either second or first-order approximation for the solvent polarization which is similar to the second-order explicit polarization applied in POSSIM force field. The FB model was first developed and parameterized with nonpolarizable OPLS-AA force field for small molecules which are not only important in themselves but also building blocks of proteins and peptide side chains. The hydration parameters are fitted to reproduce the experimental or quantum mechanical hydration energies of the molecules with the overall average unsigned error of ca. 0.076kcal/mol. It was further validated by computing the absolute pKa values of 11 substituted phenols with the average unsigned error of 0.41pH units in comparison with the quantum mechanical error of 0.38pH units for this set of molecules. There was a good transferability of hydration parameters and the results were produced only with fitting of the specific atoms to the hydration energy and pKa targets. This clearly demonstrates the numerical and physical basis of the model is good enough and with proper fitting can reproduce the acidity constants for other systems as well. After the successful development of FB model with the fixed charges OPLS-AA force field, it was expanded to permit simulations with Polarizable Simulations Second-order Interaction Model (POSSIM) force field. The hydration parameters of the small molecules representing analogues of protein side chains were fitted to their solvation energies at 298.15K with an average error of ca.0.136kcal/mol. Second, the resulting parameters were used to reproduce the pKa values of the reference systems and the carboxylic (Asp7, Glu10, Glu19, Asp27 and Glu43) and basic residues (Lys13, Lys29, Lys34, His52 and Lys55) of the turkey ovomucoid third domain (OMTKY3) protein. The overall average unsigned error in the pKa values of the acid residues was found to be 0.37pH units and the basic residues was 0.38 pH units compared to 0.58pH units and 0.72 pH units calculated previously using polarizable force field (PFF) and Poisson Boltzmann formalism (PBF) continuum solvation model. These results are produced with fitting of specific atoms of the reference systems and carboxylic and basic residues of the OMTKY3 protein. Since FB model has produced improved pKa shifts of carboxylic residues and basic protein residues in OMTKY3 protein compared to PBF/PFF, it suggests the methodology of first-order FB continuum solvation model works well in such calculations. In this study the importance of explicit treatment of the electrostatic polarization in calculating pKa of both acid and basic protein residues is also emphasized. Moreover, the presented results demonstrate not only the consistently good degree of accuracy of protein pKa calculations with the second-degree POSSIM approximation of the polarizable calculations and the first-order approximation used in the Fuzzy-Border model for the continuum solvation energy, but also a high degree of transferability of both the POSSIM and continuum solvent Fuzzy Border parameters. Therefore, the FB model of solvation combined with the POSSIM force field can be successfully applied to study the protein and protein-ligand systems in water. "

Genetic sequences are being collected at an ever increasing rate due to rapid cost reductions; however, experimental approaches to determine the structure and function of the protein(s) each gene codes are not keeping pace. Therefore, computational methods to augment experimental structures with comparative (i.e. homology) models using physics-based methods for building residues, loops and domains are needed to thread new sequences onto homologous structures. In addition, even experimental structure determination relies on analogous first principles structure refinement and prediction algorithms to place structural elements that are not defined by the data alone. Computational methods developed to find the global free energy minimum of an amino acid sequence (i.e. the protein folding problem) are increasingly successful, but limitations in accuracy and efficiency remain. Optimization efforts have focused on subsets of systems and environments by utilizing potential energy functions ranging from fixed charged force fields (Fiser, Do, & Sali, 2000; Jacobson et al., 2004), statistical or knowledge based potentials (Das & Baker, 2008) and/or potentials incorporating experimental data (Brunger, 2007; Trabuco, Villa, Mitra, Frank, & Schulten, 2008). Although these methods are widely used, limitations include 1) a target function global minimum that does not correspond to the actual free energy minimum and/or 2) search protocols that are inefficient or not deterministic due to rough energy landscapes characterized by large energy barriers between multiple minima. Our Global Optimization Using Metadynamics and a Polarizable Force Field (GONDOLA) approach tackles the first limitation by incorporating experimental data (i.e. from X-ray crystallography, CryoEM or NMR experiments) into a hybrid target function that also includes information from a polarizable molecular mechanics force field (Lopes, Roux, & MacKerell, 2009; Ponder & Case, 2003). The second limitation is overcome by driving the sampling of conformational space by adding a time-dependent bias to the objective function, which pushes the search toward unexplored regions (Alessandro Barducci, Bonomi, & Parrinello, 2011; Zheng, Chen, & Yang, 2008). The GONDOLA approach incorporates additional efficiency constructs for search space exploration that include Monte Carlo moves and fine grained minimization. Furthermore, the dimensionality of the search is reduced by fixing atomic coordinates of known structural regions while atoms of interest explore new coordinate positions. The overall approach can be used for optimization of side-chains (i.e. set side-chain atoms active while constraining backbone atoms), residues (i.e. side-chain atoms and backbone atoms active), ligand binding pose (i.e. set atoms along binding interface active), protein loops (i.e. set atoms connecting two terminating residues active) or even entire protein domains or complexes. Here we focus on using the GONDOLA general free energy driven optimization strategy to elucidate the structural details of missing protein loops, which are often missing from experimental structures due to conformational heterogeneity and/or limitations in the resolution of the data. We first show that the correlation between experimental data and AMOEBA (i.e. a polarizable force field) structural minima is stronger than that for OPLS-AA (i.e. a fixed charge force field). This suggests that the higher order multipoles and polarization of the AMOEBA force field more accurately represented the true crystalline environment than the simpler OPLS-AA model. Thus, scoring and optimization of loops with AMOEBA is more accurate than with OPLS-AA, albeit at a slightly increased computational cost. Next, missing PDZ domain protein loops and protein loops from a loop decoy data set were optimized for 5 ns using the GONDOLA approach (i.e. under the AMOEBA polarizable force field) as well as a commonly used global optimization procedure (i.e. simulated annealing under the OPLS-AA fixed charge force field). The GONDOLA procedure was shown to provide more accurate structures in terms of both experimental metrics (i.e. lower Rfree values) and structural metrics (i.e. using the MolProbity structure validation tool). In terms of Rfree, only one out of seven simulated annealing results was better than the Gondola global optimization. Similarly, one simulated anneal loop had a better MolProbity score, but none of the simulated annealing loops were better in both categories. On average, GONDOLA achieved an Rfree value 19.48 and simulated annealing saw an average Rfree value of 19.63, and the average MolProbity scores were 1.56 for GONDOLA and 1.75 for simulated annealing. In addition to providing more accurate predictions, GONDOLA was shown to converge much faster than the simulated annealing protocol. Ten separate 5 ns optimizations of the 4 residue loop missing from one of the PDZ domains were conducted. Five were done using GONDOLA and five with the simulated annealing protocol. The fastest four converging results belonged to the GONDOLA approach. Thus, this work demonstrates that GONDOLA is well-suited to refine or predict the coordinates of missing residues and loops because it is both more accurate and converges more rapidly.

Specific recognition of ligands including metal ions by proteins is the key of many crucial biological functions and systems. Accurate prediction of the binding strength not only sheds light on the mechanism of the molecular recognition but also provides the most important prerequisite of drug discovery. Computational modeling of molecular binding has gained a great deal of attentions in the last few decades since the advancement of computer power and availability of high-resolution crystal structures. However there still exist two major challenges in the field of molecular modeling, i.e. sampling issue and accuracy of the models. In this work, I have dedicated to tackling these two problems with a noval polarizable force field which is believed to produce more accurate description of molecular interactions than classic non-polarizable force fields. We first developed the model for divalent cations Mg²⁺ and Ca²⁺, deriving the parameters from quantum mechanics. To understand the hydration thermodynamics of these ions we have performed molecular dynamics simulations using our AMOEBA force field. Both the water structures around ions and the solvation free energies were in great accordance with experiment data. We have also simulated and calculated the binding free energies of a series of benzamidine-like inhibitors to trypsin using explicit solvent approach by free energy perturbation. The calculated binding free energies are well within the accuracy of experimental measurement and the direction of change is predicted correctly in all cases. Finally, we computed the hydration free energies of a few organic molecules and automated the calculation procedure. / text

L’angiogenèse étant un processus limité dans des conditions physiologiques et un processus clé dans la croissance tumorale, elle est devenue une cible thérapeutique prometteuse. La neuropiline-1 (Np1) est un corécepteur du VEGF, qui est le facteur pro-angiogénique le mieux décrit jusqu’à présent. Dans cette thèse, nous nous intéressons à la modélisation du ligand 47, une molécule active expérimentalement, en vue de l’amarrer dans la neuropiline-1. En raison de sa flexibilité conformationnelle, ce ligand pourrait en effet adopter une conformation étendue, comme la tuftsine, un ligand naturel de Np1, ou une conformation repliée. Une étape préliminaire essentielle, avant l'exploration de la complexation du lig-47 a Np1, est l'étude de sa flexibilité conformationnelle. Il est en effet impératif de s'assurer que les interactions intra-moleculaires (conformationnelles) dans le lig-47 sont calculées a une precision comparable a celles de ses interactions intermoléculaires avec Np1. Un écueil important tient au caractère conjugue du lig-47. Les quatre fragments constitutifs de cette molecule sont tous aromatiques ou conjugues, et sont tous connectes par des atomes insaturés. [...] Nous avons, dans une première étape, construit la molécule en quatre fragments : le benzimidazole, le méthylbenzène, le benzodioxane et le carboxythiourée (CTU). Nous avons approché ces quatre fragments par une molécule d’eau afin de calibrer les rayons de Van de Waals effectifs impliqués dans les contributions d’énergies électrostatique et de répulsion. L’anisotropie ayant été assurée, nous avons cherché à reproduire la conjugaison, par la calibration des barrières de torsion V0 primaires (n=1) et binaires (n=2). Nous obtenons des accords très satisfaisants entre les courbes conformationnelles obtenues avec SIBFA et celles obtenues par des calculs de chimie quantique (QC). Nous avons ensuite effectué des minimisations de l'énergie SIBFA en partant des minima des six courbes conformationnelles. Afin d'évaluer la transférabilité de la méthode, nous avons comparé les stabilités relatives de ces minima par des calculs QC ponctuels. Or les différences d'énergie séparant le minimum « global » des minima locaux se sont avérées sous-estimées par rapport aux calculs QC. Ces résultats nous ont amenés à envisager des façons différentes de représenter le fragment CTU. La possibilité la plus évidente consiste à le séparer en deux sous-fragments amide et thioamide. Les effets de la conjugaison et de la transférabilite des multipôles et polarisabilités sont ainsi perdus mais pourraient être compensés par la prise en compte explicite de l’énergie de polarisation des fragments amide et thioamide. Avec cette approche, la recalibration des rayons effectifs a permis de préserver des accords convenables avec les calculs quantiques pour l'approche des atomes du CTU par une molécule d'eau sonde. Les courbes conformationnelles reproduisent de près les courbes QC avec une recalibration minimale. Les minima de ces courbes ont été a nouveau minimisés en SIBFA, conduisant a des structures néanmoins très proches des minima correspondants de l'approche précédente avec un CTU construit d'un seul tenant. Mais à présent, les différences d'énergie séparant le minimum global des minima locaux sont très voisine de celles trouvées en QC. De plus, l'évolution des courbes conformationnelles en fonction de la structure considérée s'est avérée régie par l'énergie de polarisation. Par ailleurs, nous avons obtenu des résultats satisfaisants lors de l’amarrage de la tuftsine dans Np1. Ces résultats s'avèrent suffisamment probants pour permettre d'envisager à présent une étude détaillée des modes d'interaction du lig-47 avec Np-1. / Angiogenesis is a limited process in physiological conditions and a key process in tumor growth. Hence, it has become a promising therapeutic target. The neuropilin-1 (Np1) is a co-receptor for VEGF, which is today the best known pro-angiogenic factor. This manuscript deals with the molecular modelling of the ligand 47 (lig47), an experimentally active molecule, in order to dock it into Np1. Due to its conformational flexibility, this ligand could adopt and extended conformation such as tuftsin, a natural ligand, does in its complex with Np1, or a folded conformation. An essential preliminary step, before the exploration of the complexation of lig47 with Np1 is the study of its conformational flexibility. Indeed it must be ensured that the intramolecular (conformational) interactions are calculated with precision compared with the calculations of the interaction energies with Np1. An important issue comes from the polyconjugaison of the lig47. [...] We probed these four fragments with a water molecule in order to calibrate the effective Van der Waals radii implicated in electrostatic and repulsion contributions. Once the anisotropy was reproduced, we looked for reproducing the effect of the conjugaison on torsional barriers. We hence calibrated primary (n=1) and binary (n=2) barriers. We obtained very satisfactory agreements between the conformational curves obtained with SIBFA and those obtained with quantum chemistry (QC) calculations. Then we performed energy minimizations of SIBFA energy of the minima of the ix conformational curves. In order to evaluate the transferability of the method, we compared the relative stabilities of these minima with single-point QC calculations. But the differences of energy between the global minimum and local minima were underestimated in comparison with QC calculations. These results lead us to consider different ways of representing CTU. The more evident way consists in separating it in two fragments amide and thioamide. The effects of the conjugaison and the transferability of the multipoles and the polarisabilities are lost but could be compensated by the explicit consideration of the polarization energy between the amide and thioamide fragments. With this approach, the recalibration of the effective radii permitted to preserve good agreements with the QC calculations when probing the CTU by a water molecule. The conformational curves reproduce the QC curves after a minimal recalibration. The minima of these curves were then re-minimized with SIBFA, leading to structures close to those obtained with the first representation. But now, the differences of energy between the global minimum and the local minima are very close to those obtained in QC. Moreover, the evolution of the conformational curves as function of the number of the structure is ruled by the polarization energy. Otherwise, we obtained satisfactory results when we docked the tuftsin in Np1. These results will allow us to consider a detailed study of the interaction between lig47 et Np1.