Keywords: Computer simulations to engineer PDZ-peptide recognition Computational Protein Design, Molecular dynamics, Monte Carlo simulation, Polarizable force field, Free energy calculation, PDZ Dessin de protéine, dynamique moléculaire, simulation Monte Carlo, champ

).

Basdevant, Tâp Ha Duong et Georgios Archontis pour avoir accepté d'examiner mon travail de thèse. Je voudrais adresser un remerciement particulier à mon directeur de thèse, Thomas Simonson. Merci d'avoir accepté ma candidature il y a trois ans et d'avoir suivi l'avancement de mon travail jour après jour avec une rigueur scientifique extraordinaire. Merci pour les nombreuses discussions, les conseils et les enseignements. Merci enfin pour le temps consacré à lire et commenter ce manuscrit.

Je remercie particulièrement Ernesto J. Fuentes pour son importante collaboration dans le projet PDZ, Alex MacKerell et Benoît Roux pour leur grande disponibilité et soutien avec le développement du champ de force Drude.

Je remercie tout l'équipage du laboratoire BIOC, surtout le directeur Yves Méchulam pour son soutien pendant le processus de sélection de l'école doctorale. Merci à Thomas G. et David pour le support technique et scientifique pendant ces trois ans. Merci à Alexey pour les discussions et tous les pauses café, à Savvas et Karen pour m'avoir aidé à faire les premiers pas dans le monde du CPD avec Proteus. Merci aux camarades Nicolas P. et Vaitea, pour votre soutien pendant les jours de travail et après les nombreux bières partagé ensemble. Merci à Maxime et Yann. Je remercie aussi Rojo, Alexandrine, Marie-Pierre et à tous ceux qui ont participé temporairement aux activités de l'équipe de bioinformatique. Merci à Paula, et bon courage pour ta thèse. Je remercie Clara, Pierre, Marc D. pour les pauses déjeuner. Un grand merci à Mélanie et Caroline pour toute votre aide. Je remercie aussi tous les autres membres du labo: Pierre-Damien, Emma, Marc G., Alexis, Titine, Mimi, Gabrielle, Michou, Nathalie et Guillaume. Merci aux doctorants et postdoctorants Giuliano, Nicolas M., Amlam, Ditipriya, Irène, Ramy, Nhan et Angelina.

Merci aux amis, particulièrement aux ritals pour tous les soirées parisiennes passées ensemble et aussi pour les futures. Ovunque saremo fra tre anni, avremo sempre Parigi. Merci Chiara, pour tout ce que tu as partagé avec moi jusqu'ici. Merci à ma belle famille, à Roberto et Ivana, pas besoin d'expliquer pourquoi. E infine grazie a Mamma, Papà e Stella, perché senza di voi non sarei mai arrivato fin qui.

All the best. i

Table of contents

List of figures vii

List of tables ix

Introduction

Protein-protein interactions (PPIs) regulate complex functional networks in cells. They are often mediated by specialized domains such as SH3 and PDZ domains. The specific recognition between such domains controls important biological functions. The understanding of these phenomena is extremely difficult and challenging, even when the interaction partners are relatively small protein domains. Indeed, PPI networks exploit a large number of binding events to transfer information through communication pathways, governing signalling inside and outside cells. Disrupting or altering the equilibrium between PPIs plays an important role in several diseases: one example is the observed altered expression of important signalling proteins found in diseased tissues. The inibition of targeted PPIs is a recognized strategy for the development of new drugs [START_REF] Böhm | Protein-ligand interactions: From molecular recognition to drug design[END_REF]).

Understanding and engineering PPIs can be addressed by molecular modeling approaches, including free energy simulations and Computational Protein Design. Computational Protein Design (CPD) allows one to perform large screenings of many protein variants in a limited amount of time. Several protein mutants can be studied starting from a native structure, optimizing a preferred quantity such as stability or binding. CPD applications increased dramatically in the last ten years, but there is a strong need of continuous improvement of both theoretical and computational models to obtain reliable and predictive results that match experiments as closely as possible.

In the present thesis we focus on PDZ domains, which are among the most widespread PPI domains. The human genome has more than 400 PDZ domains belonging to more than 200 different proteins. They are also present in bacteria, yeasts, plants, insects and vertebrates (from the SMART database -Schultz [2000] ; [START_REF] Letunic | 20 years of the SMART protein domain annotation resource[END_REF]). PDZs are small domains, composed of roughly 90 amino acids that specifically recognize the 5 -10 C-terminal amino acids of their binding partners. They are often able to bind the corresponding peptides in isolation. The experimental screening of binding of peptide and protein variants is the most common strategy to better characterize their interactions [START_REF] Songyang | Recognition of unique carboxyl-terminal motifs by distinct PDZ domains[END_REF] ; Wiedemann et al. [2004] ; Stiffler et al. [2007]).

Introduction PDZ domains share a common secondary structure, composed of 5-6 β strands and two α helices, as shown in figure 1. In their typical bound configuration, C-terminal residues of the binding partner assume an elongated β strand conformation in the binding pocket between the strand β 2 and helix α 2 .

Figure 1 -PDZ domain structure. Left: with labels indicating α helices and β strands from Nter to Cter residues, colored in green and blue. The peptide ligand is colored in orange and interacts with residues in the binding pocket situated between β 2 and α 2 . Right: with labels indicating ligand positions, using negative numbering: position P 0 is the most buried and corresponds to the peptide Cter residue.

Several PDZs of different amino acid sequences show similar affinities for the same ligand [START_REF] Ernst | Coevolution of PDZ domain-ligand interactions analyzed by high-throughput phage display and deep sequencing[END_REF]). This suggests that important hotspot residues for their PPIs are situated in specific domain regions. Crystal structures and combinatorial screening showed that strands β 2 , β 3 as well as helix α 2 are important for the PDZ:peptide binding [fig. 1]. PDZs can be classified according to their specificity for different peptide ligands. We number the C-ter residues of the peptide ligand in inverse order where position P 0 is the most buried C-ter residue, P -1 is the second-last and so on. Four PDZ classes are defined according to the composition of the last 3 C-terminal positions of the peptide ligand, P 0,-1,-2 . The first class binds the pattern S/T -X -Φ (Serine/Threonine followed by an amino acid X and then a hydrophobic residue). The second class binds Φ -X -Φ (any residue between two hydrophobic amino acids). The third binds D/E -X -Φ and X -Ψ -D/E where Ψ is an aromatic residue. Another PDZ classification uses 16 patterns based on the last seven Cter residues (Tonikian et al. [2008]). In this study the authors analyzed binding specificity on a large scale. A huge combinatorial peptide library was used to produce more than 3000 peptide ligands able to bind to 54 PDZ domains from the human genome and 28 from the Caenorhabditis elegans worm genome. PDZs showed capability to bind many different peptides, with patterns which are conserved between worm and human. All these experimental strategies need expression and purification of PDZ domains: the procedure can be quite expensive especially when considering several mutants. Also the structural characterization of variants can be quite difficult and time-demanding. For this reason experiments have been coupled with computational studies, often involving molecular dynamics (MD) simulations which are also able to give insights about thermodynamic properties of the systems of interest. [START_REF] Basdevant | Thermodynamic basis for promiscuity and selectivity in protein-protein interactions: PDZ domains, a case study[END_REF] performed MD simulation of 12 different PDZ:peptide systems and evaluated binding contributions. Comparing the simulated complexes, they found that non-polar interactions are predominant for specificity, as expected from the observed promiscuity. The authors also found the importance of dynamical behavior and peptide reorganization upon binding. Other studies can be found in Tian et al. Among the known PDZs, we will focus on the domain belonging to the Tiam1 protein. It is known for modulating signaling activity for cell proliferation and migration, whose dysregulation increases metastatis growth in cancer. Indeed, both overexpression and subexpression of Tiam1 are connected to several tumors. For example, overexpression was found in sick tissues [START_REF] Minard | Tiam1 regulates cell adhesion, migration and apoptosis in colon tumor cells[END_REF]), suggesting that it can be used as a marker [START_REF] Zhao | Overexpression of t lymphoma invasion and metastasis 1 predict renal cell carcinoma metastasis and overall patient survival[END_REF] ; [START_REF] Ding | Overexpression of tiam1 in hepatocellular carcinomas predicts poor prognosis of HCC patients[END_REF]). Starting from now, we will use "Tiam1" to indicate the PDZ domain of the Tiam1 protein. It is a class 2 PDZ domain, meaning that it has a preference to bind Φ -X -Φ peptides (any residue between two hydrophobic amino acids at positions P 0,-2 ). Four X-ray structures of Tiam1 complexes where the domain is bound to different ligands can be found in the PDB with access codes 3KZD, 3KZE, 4GVC and 4GVD [START_REF] Shepherd | The tiam1 PDZ domain couples to syndecan1 and promotes cell-matrix adhesion[END_REF]). Structures are also known for three quadruple mutants (4NXP, 4NXQ and 4NNR, Liu et al. [2013]) and one NMR structure has also been solved for the unbound, Apo form of Tiam1 (2D8I, [START_REF] Qin | Solution structure of the PDZ domain of t-cell lymphoma invasion and metastasis 1 varian[END_REF]). Tiam1 is composed of 94 residues. An experimental screening of an 8-residue peptide library [START_REF] Shepherd | Distinct ligand specificity of the tiam1 and tiam2 PDZ domains[END_REF]) showed its capability to bind 70 different peptides. For selected ligand variants, their relative binding free energies were measured. A 2.2 kcal/mol interval was observed between the best and worst binding free energies.

Introduction

In the first part of this manuscript, we begin with a general introduction about Computational Protein Design (CPD), with more detailed information regarding the Proteus software [START_REF] Simonson | Computational protein design: The proteus software and selected applications[END_REF]) developed at École Polytechnique. Next, we apply the method to design peptide sequences that bind the Tiam1 PDZ domain. Results were published (Villa et al. [2018b]) and the paper can be found attached to chapter two.

A more detailed and quantitative study of Tiam1:peptide binding will be described in the second part of this thesis. We performed polarizable free energy simulations to study relative binding free energies for charge mutations in the binding pocket, using the Drude polarizable force field. Results were compared to those obtained with two additive force fields (Amber and Charmm). General information about polarizable models with a particular interest on Drude will be given in chapter three. In chapter four, we describe more technical aspects related to Drude free energy simulations using the a dual topology approach. Our results were published in a larger study [START_REF] Panel | Accurate PDZ/peptide binding specificity with additive and polarizable free energy simulations[END_REF]).

In the third part of the thesis, we describe the development of a polarizable Drude model for methyl phosphate and phosphotyrosine using established Drude parametrization procedures. Methyl phosphate parameters were validated computing Mg 2+ standard binding free energies which showed good agreement with experiments. Results were recently published (Villa et al. [2018a]).

Finally, we also performed a CPD study of acid:base equilibrium constants for a collection of proteins (more information can be found in [START_REF] Villa | Comparing pairwiseadditive and many-body generalized born models for acid/base calculations and protein design[END_REF]). We compared two approximations used to compute electrostatic interactions in Proteus, which both rely on a Generalized Born (GB) implicit solvent model. The more sophisticated model could be used in future design work.

Chapter 1

Computational design of PDZ-peptide binding

The protein design problem can be defined as looking for the amino acid sequences that fit a specific structure according to a scoring function which depends on the system energy. Random mutagenesis is used to produce a set of candidates that are subjected to selective pressure. Since CPD must be efficient, combinatorial approaches are preferred. A simplified and discretized conformational space, a simple scoring function and a fast exploration algorithm are key ingredients of all CPD software.

The starting point of a typical CPD simulation is the three-dimensional structure of a given protein. The wildtype sequence associated to this native conformation can then be altered by modification of side chain atoms belonging to one or more residues. This is done efficiently if the conformational space is discretized, usually keeping the protein backbone fixed. An energy function is used to score different sequences generated by the sampling algorithm. The conceptual difference between energy function and sequence score is important: the first is the formula used to compute interactions between groups of atoms in a specific configuration. The latter depends on the system energy but can vary according to the quantity to optimize. For example, optimizing protein stability one uses a scoring function which is based on the unfolding free energy difference upon mutation. On the other hand, to optimize ligand binding one has a scoring function which is based on the change in binding affinity. To optimize protonation states, one has to use a scoring function based on protonation free energy differences, and so on.

In this chapter we start by describing the most popular energy functions in CPD, how the conformational space is usually discretized and some widely-used exploration routines. Then we will give more detailed information about the Proteus software and how sequences are scored according to different quantities, with a special interest in ligand binding.

Chapter 1. Computational design of PDZ-peptide binding 1 Energy functions for CPD CPD energy functions are based on molecular mechanics plus an implicit solvent model. The first is borrowed from an existing force field for molecular simulation (for example Amber, Charmm) while the latter is often based on simple dielectric shielding. One widely-used implicit solvent model is Generalized Born (GB) combined with solvent accessible surface area (SA), also called GB/SA.

Several CPD softwares use optimized, statistical potentials, where the physical-based energy function terms are weighted or scaled by empirical parameters in order to fit experimental data. However in the present discussion we will focus on physical potentials, in the spirit of our inhouse CPD software Proteus. We start by describing molecular-mechanics for internal degrees of freedom, then nonbonded interaction terms. Bonded energy term Interactions between bonded atoms (up to 1-2, 1-3, 1-4 couples) are computed using a sum of several contributions

E bonded = E bond + E angl + E U rey-Bradley + E dihedral + E improper (1.1)
Most of the terms share the same, harmonic potential functional form (bonds, angles, Urey-Bradley, Impropers) while the energy associated to dihedral angles is usually periodic

E bond = k b (b -b 0 ) 2 E angl = k α (α -α 0 ) 2 E U -B = k u (u -u 0 ) 2 E impr = k ω (ω -ω 0 ) 2 (1.
2)

E dih = k φ [1 + cos (nφ -δ)]
In general, bond-angle-dihedral terms are applied to all 1-2, 1-3, 1-4 couples [fig. 1.1], while Urey-Bradley and Impropers are used in more special cases [fig. 1.2]. Urey-Bradley controls angle bending and is applied between a few 1-3 couples which form an angle α with a central Bonded parameters (like bond distance b 0 and force constant k b ) usually depend on atom types. Their values are optimized in order to reproduce quantum-mechanical (QM) optimized geometries or experimental structures, as well as QM vibrational spectra and QM energy scans. In most of the existing force fields, bonded parameters are quite transferable between different molecules when they refer to the same chemical group.

Nonbonded energy term Nonbonded interactions are not computed for 1-2 and 1-3 pairs (their interactions are completely included in internal terms). 1-4 couples can be excluded or scaled. The interaction between all other not excluded pairs of atoms is then sum of two terms

E nonbonded = E LJ + E elec (1.3)
Van der Waals interactions between atoms of type i and j at distance r ij are expressed in the attractive part of the Lennard-Jones potential energy

E LJ (r ij ) = 4ǫ ij   σ ij r ij 12 - σ ij r ij 6   (1.4)
where ǫ ij and σ ij are force field constants for a couple of atom types. The first term expresses the strong nuclear repulsion that occurs at short distances because of the Pauli principle.

The force field parameters ǫ ij and σ ij are optimized targeting QM interaction energies and experimental consended-phase data. As for the corresponding bonded interaction terms, LJ parameters are usually transferable between different molecules.

Partial atomic charges define the magnitude of Coulomb electrostatic interactions in a dielectric ǫ between atoms i and j at distance r ij E elec (r ij ) = q i q j ǫr ij (1.5)

Chapter 1. Computational design of PDZ-peptide binding

In contrast to bonded and LJ parameters, atomic charges are less transferable. They depend on the specific molecule and are usually parametrized targeting its QM dipole moment, as well as interactions between water and model compound acceptors/donors.

Implicit solvent models

Implicit solvent models are based on the decomposition of the system solvation free energy ∆G solv by integration over the solvent degrees of freedom. In general, the total ∆G solv is split into two parts: the first counts electrostatic interactions between charged atoms and the second counts the nonpolar, non-electrostatic interactions:

∆G solv = ∆G solv elec + ∆G solv apol (1.6)
The solvation process is decomposed in 3 steps: 1) formation of a solvent cavity which depends on the solute volume/shape, without solute-solvent interactions 2) turn on solutesolvent vdW dispersion interactions, 3) turn on electrostatic interactions. Steps 1 and 2 define ∆G solv apol while step 3 defines ∆G solv elec . We now describe how popular implicit models used in CPD express the different terms of equation 1.6.

Surface-accessible Surface Area

In SASA models, the non-polar solvation free energy is written as a linear combination of contributions from solute atoms:

∆G solv surf = i σ i S i (1.7)
For each atom i the solvent accessible surface S i [fig. 1.3] is multiplied by a coefficient σ i . It should be noted that equation 1.7 is based on a strong assumption of additivity.

probe w

Energy functions for CPD

CASA model In the Coulombic Accessible Surface Area model, the total (polar+apolar) solvation free energy is expressed using the simple function

∆G solv = 1 ǫ -1 E elec + ∆G solv surf (1.8)
The apolar contribution to the solvation free energy is obtained using the surface energy term (SASA model), while the electrostatic contribution is obtained scaling Coulomb interactions. This is one of the first implicit solvent models used in CPD. Despite its simplicity, successful applications can be found in the literature, for example in the work of Dahiyat & Mayo [1997].

Generalized Born model To describe continuum electrostatics in the context of GB, we follow the analytical treatment of [START_REF] Schaefer | A comprehensive analytical treatment of continuum electrostatics[END_REF]. Their integral formula for the system electrostatic energy

E el = 1 8πǫ w w D2 (x)d 3 x + 1 8πǫ p p D2 (x)d 3 x (1.9)
is based on the assumption that the system is split into two regions (solvent and solute) characterized by two distinct dielectric constants ǫ w and ǫ p . Only solute atoms i contribute to the additive dielectric displacement D = N i=1 Di .

q i q j p w b i b j

Figure 1.4 -Interaction between solute charges Left: potential at point M generated by a charge q buried in the solute interior. Right: Solvation radii of two charges q i and q j , interacting in the solute ǫ p and screened by the solvent ǫ w .

One can write It is possible to manipulate the D2 integrals of equation 1.10. The full calculation is shown in [START_REF] Schaefer | A comprehensive analytical treatment of continuum electrostatics[END_REF]. One can split diagonal D2

E el = 1 8πǫ w R 3 D2 (x)d 3 x + τ 8π p D2 (
i and off-diagonal Di Dj elements and then sum up results in self-energy and interaction-energy terms. Separating the standard, wellknown Coulomb interaction between charges in the solute environment and the contribution from the solvent one obtains E el = i<j q i q j ǫ p r ij + ∆E solv (1.12)

∆E solv = i ∆E self i + i<j ∆E int ij (1.13)
Using the approximation of Still et al. [1990], the solvent contribution is computed using

∆E int ij = - τ q i q j r 2 ij + b i b j exp -r 2 ij /4b i b j 1/2 (1.14)
where b i is called solvation radius of the solute charge i. It has a physical interpretation, as the measure of the distance between charge i and the solute surface. It is an analytical function of all other solute atoms positions b i = b i (r 1 ,r 2 ,...,r n ) (1.15)

1. Energy functions for CPD

.5 shows the limiting behaviour of this interaction.

The accuracy of all GB models relies on the accuracy of solvation radii calculation. In several applications, they are obtained from the definition of the self-energy potential

∆E self i def = -τ q 2 i 2b i (1.16)
Using this definition, solvation radii are obtained from the calculation of associated selfenergies. To calculate this quantity, it is necessary to approximate the integral of the squared dielectric displacement using a continuos density function P k (x) to describe the solute volume

∆E self i = -τ q 2 i 2ǫ w R i + τ q 2 i k =i E self ik = -τ q 2 i 2ǫ w R i + k =i τ 8π R 3 D 2 i P k (x)dx (1.17)
where R i is an atomic radius, part of the force field parameters.

We now describe two GB variants that differ in the method used to approximate equation 1.17. In GB/ACE, the integral is calculated using a density function P k (x) where volume occupied by charges is expressed with gaussians

E self ik = 1 ω ik exp(-r 2 ik σ 2 ik ) + V k 8π r 3 ik r 4 ik + µ 4 ik (1.18)
ω ik and µ ik are functions of the volume V k , of the atomic radii R k = (3V k /4π) and a smoothing parameter λ. Another approximation of equation 1.17 is GB/HCT, where

4E self ik = 1 L ik - 1
U ik + r ik 4 1 U 2 ik - 1
L 2 ik + 1 2r ik ln L ik U ik + R 2 ik 4r ik 1 L 2 ik - 1
U 2 ik (1.19)
and

L ik =            1 r ik + R k ≤ R i R i r ik -R k ≤ R k ≤ r ik + R k r ij -R k R i ≤ R k ≤ r ik -R k U ik =      1 r ik + R k ≤ R i r ik -R k R i ≤ r ik + R k (1.20)
Onufriev proposed a slightly modified version of equation 1.16 to compute b i . The atomic radius R i of each atom is reduced by a constant and the self-energy is scaled by a factor λ, Chapter 1. Computational design of PDZ-peptide binding then the solvation radius is reduced by

δ b i =   (R i -ρ 0 ) -1 -λ k =i E self ik   -1 -δ (1.21)
In the Proteus GB implementation, λ = 1.4, ρ 0 = 0.09 Å and δ = 0.15 Å.

LK model

In the Lazaridis-Karplus model [START_REF] Lazaridis | Effective energy function for proteins in solution[END_REF]), the free energy of solvation is a sum of atomic terms

∆G solv i ∆G solv = i ∆G solv i (1.22)
where the sum runs over atoms and

∆G solv i = G ref i - j V j f i (r ij )dV ≃ G ref i - j =i f i (r ij )V j (1.23)
The "reference" solvation free energies G ref i measure the change in solvation free energy for a process where the atom i is transferred from a fully-solvated, unfolded state to a partially solvated or buried folded state. Each integral is approximated by volumes V j occupied by the other atoms j. The f i terms are gaussian energy density functions

4πr 2 ij f i (x ij ) = α i exp(-x 2 ij ) x ij = r ij -R i λ i (1.24)
where R i is the Van der Waals radius of atom i, r ij is the distance between i and j and λ i is a parameter called gaussian correlation length such that when i is fully solvated, the total

∆G solv i equals zero.
Dispersion term In most popular implicit solvation models, the apolar contribution to the solvation free energy in the three steps process [eqn. 1.6] has a cavity-formation term calculated with the accessible surface area model. Moreover, the attractive part of the LJ potential [eqn. 1.4] can be integrated for each solute atom, in order to obtain total solute-solvent vdW dispersion

∆G solv disp = i U vdW i = i -ρ w V solv 4ǫ iw σ 6 iw (r -r i ) 6 dV (1.25)
In the formula, ǫ iw and σ iw are the nonbonded parameters of eqn. 1.4 between solute atom i and the water oxygen. The integral is performed over the solvent volume, considering the number density of water molecules ρ w as a constant. It can be split in two parts

U vdW i = - 4π 3 ρ w 4ǫ iw σ 6 iw 1 R 3 i - 3 4π solute V i 4ǫ iw σ 6 iw (r -r i ) 6 dV (1.26)
where V i is the solute volume excluding atom i. Onufriev and coworkes proposed an approximation of the integral: the solution consists in two terms

I disp i = 3 4π V i dV (r -r i ) 6 ≃ I main i + I neck i (1.27)
where main and neck are respectively contributions relative to the volume region occupied by atomic spheres and the volume region between atomic spheres. The two integrals are obtained analytically (Aguilar et al. [2010]).

Structural models for CPD

Discretization of the protein conformational space Using an implicit solvent model, the cost to compute interactions is aggressively reduced. Protein degrees of freedom must also be simplified in order to explore large ensembles of protein variants in an acceptable computing time. Most CPD implementations use a fixed backbone while side chains assume a discrete set of conformations, or rotamers. These are defined by a rotamer library, which can be backbone dependent [START_REF] Dunbrack | Backbone-dependent rotamer library for proteins application to side-chain prediction[END_REF]) or not [START_REF] Tuffery | A new approach to the rapid determination of protein side chain conformations[END_REF]). The use of rotamer libraries is justified by the fact that amino acid side chains in proteins assume a limited number of configurations, as confirmed by several studies (Finkelstein & Ptitsyn [1977] ; Janin et al. [1978] ; [START_REF] Ponder | Tertiary templates for proteins[END_REF]). CPD models can also use native rotamers, extending existing libraries by introduction of side chain configurations extracted from the crystallographic structure of the system of interest.

Pairwise-additive energy function

To explore vast ensembles of structures and sequences, CPD software needs to rapidly compute interactions between couples of atoms in a discretized conformational space. In many situations it is convenient to follow a combinatorial approach. The protocol is based on a pre-calculated energy matrix, where pairwise-additive interactions between groups of atoms (in general between couples of residues) are stored and can be accessed during the sequence/structure exploration (Dahiyat & Mayo [1997] ; [START_REF] Gaillard | Pairwise decomposition of an MMGBSA energy function for computational protein design[END_REF]). Additional information will be given below, with more details about our in-house software Proteus.

We describe now the general case, where such an energy function would be

E tot = i E ii + i<j E ij (1.28)
The sum runs over residues i and j which occupy one of their possible rotamers. Using equation 1.28, a problem arises when one wants to define a pairwise-additive energy function based on molecular-mechanics, where the solvent is described implicitly with one of the previously Chapter 1. Computational design of PDZ-peptide binding ). Another problem arises considering the calculation of the solvent accessible surface (important when using CASA or GBSA). Indeed, considering a couple of residues i and j the intersection between their surfaces can also include their intersection with a third residue, as well as all other residues in the system [fig. Heuristic methods Heuristic methods are useful to determine a group of low energy configurations starting from one or more random states. A popular method was introduced by Wernisch et al. [2000]. They proposed a simple sampling method, where starting from a random protein sequence, a single position was picked and and its rotamer optimized. The procedure was repeated for many iterations, until a local minimum was found. The method successfully reproduced side chain conformations for surface residues and stability changes for mutations applied in the protein core or at its surface. However, this method generally produces only a few variants. For this reason it is not particularly adapted to perform high-throughput design, where one wants to generate a distribution of many protein variants with a single simulation.

Monte-Carlo methods Monte Carlo (MC) sampling can be used to obtain a set of protein sequences generated from a stochastic process. One popular Monte Carlo method suitable for protein design is based on the Metropolis algorithm. With MC, it is possible to generate a distribution of protein sequences that are distribuited according to a particular probability density function. The system energy is used to accept or refuse new system configurations that populate the desired distribution. For MC, convergence is assured only for a very long simulation, and the sampling can be stuck in local minima. However, several advanced sampling techniques can be employed to avoid too long simulations and to jump free energy barriers. More technical information will be given below, with a detailed description about implementation in our software Proteus.

The main advance of Monte Carlo is that is usually simple to implement and can be easily adapted to many different problems. Sampling Boltzmann probability, is also possible to extract statistical-mechanics properties (for example free energy differences) which can be easily related to experimental data or to results of molecular dynamics simulations.

CPD softwares

Before introducing our in-house software Proteus, we briefly describe two CPD programs used in the lab during the last few years.

Rosetta is a collection of programs suitable for molecular modelling developed by a large community of researchers (RosettaCommons is an organization that counts 150 developers around the world). RosettaDesign is a utility used for CPD of protein stability; other tools like RosettaDock or RosettaAbInitio are used to predict conformations of protein-ligand complexes or de novo protein structure prediction. Several Rosetta energy functions are inspired Chapter 1. Computational design of PDZ-peptide binding by molecular mechanics but make extensive use of statistical terms used to fit experimental data. The first versions were based on a fixed backbone and a backbone-dependent rotamer library, while models with limited backbone flexibility were introduced more recently. Despite its remarkable success in the protein design field, the Rosetta energies are usually expressed in unphysical, rosetta units. Even if some effort has been made to translate them to kcal/mol (Alford et al. [2017]), results are still difficult to interpret, especially when compared to those obtained with experiments or state-of-the-art molecular simulations. Moreover, for the design procedure the developers preferred fast algorithms (for example Monte Carlo with Simulated Annealing) which are able to generate a few variants in a limited computer time, rather than generating many protein variants with a single run. The user interested in generating an ensemble of sequences often needs to run several independent Rosetta simulations and then discard repeated variants. A typical Proteus simulation is composed of three steps: 1) model construction, 2) precalculation of residue-residue interactions using a pairwise-additive energy function and then 3) fast exploration of structure-sequence space [fig. 1.7]. The first steps are performed starting from an input structure, given in PDB format. The initial model usually needs to be refined by adding missing hydrogens and performing few energy minimization steps. The system energy is expressed with an energy function which uses classical molecular mechanics plus an implicit solvent based on GB/SA as described later.

The Proteus CPD framework

A given protein composed of n amino acids has sequence S = {t 1 ,t 2 ,...,t n } where t i is the amino acid type at site i. In a Proteus simulation, selected residues i called active positions are able to mutate. They can vary their type t i within a mutation space T i which is defined a priori by the user. The only exception is for glycines and prolines, usually not mutable and modelled as part of the fixed protein backbone. Considering all active positions with their mutation space, the protein is then able to explore a sequence space which is defined by the ensemble of all possible protein sequences. This is a usually a very large number: as an example, for a mutation space of 20 possible types for N active positions one has 20 N possible sequences.

Folded and unfolded states

For each visited sequence we consider its folded and unfolded states, such that for two sequences s and s ′ the associated folding energy difference is

∆E = ∆E f (s → s ′ ) -∆E uf (s → s ′ ) (1.29)
The backbone conformation is the one of the initial model, given in input. We refer to it as native backbone. Side chains can assume a limited number of configurations belonging to a library of rotamers described by different torsional angles. By default Proteus uses the rotamer library from Tuffery (1995Tuffery ( or 2003)), but additional rotamers can easily be included.

C β atoms of each site are fixed together with the backbone, while the rest of the side chain can vary. For a given amino acid type t i (for example ASP or GLU shown in figure 1.9) rotamers are defined by a group of different torsional angles χ that belong to the rotamer space R i . Active positions can assume rotamers of different amino acid types, while the other positions, called inactive, can vary their conformation within rotamers of their fixed, native amino acid type. Another capability of Proteus is to include ligand rotamers, specified by the user, or multiple backbone conformations for desired protein segments. These backbone rotamers are usually obtained from molecular dynamics simulation snapshots of the wildtype protein. Folded state interactions The system energy in the folded state is described by a pairwiseadditive function

E f (s,{r i }) = n i E ii (t i ,r i ) + n i<j E ij (t i ,r i ,t j ,r j ) (1.30)
the sum runs over all the residues i = 1, ..., n. Diagonal terms include the interactions between atoms belonging to side chain i (of type t i and rotamer r i ) with itself and with the fixed backbone. Off-diagonal terms contain the interactions between atoms belonging to two different side chains i and j.

The Proteus CPD framework

Unfolded state interactions On the other hand, the unfolded state energy is a sum of diagonal terms:

E uf (S) = i E uf i (t i ) (1.31)
It corresponds to a simple physical model of non-interacting sites. Each E uf i (t i ) term is called the unfolded energy of amino acid type t i at position i. Unfolded energy values depend on the specific problem of interest, where they assume a different physical meaning. For example, to study protein stability, unfolded energies of different amino acid types are calculated using fully solvent-exposed amino acid model compounds. Moreover, the unfolded energy of an amino acid type can be the same for all protein positions i = 1, ..., n or it can vary for different sites. This flexibility allows to define particular models for the unfolded state. For example, one can use a set of unfolded energies for exposed positions and a different set for buried positions. Using equation 1.29, this choice allows to better represent the energetic cost in transferring side chains from pure solvent to different dielectric environments.

Proteus Energy Function

Interactions are computed using molecular-mechanics and an implicit solvent

E = E M M + E solv (1.32)
where

E M M = E bond + E angl + E dih + E impr + E vdw + E Coul (1.33)
Interactions between solute atoms are computed borrowing parameters from the Amber ff99SB force field

E vdw (i,j) + E Coul (i,j) = 4ǫ ij   σ ij r ij 12 - σ ij r ij 6   + q i q j ǫ p r ij (1.34)
where ǫ p is the solute (protein) dielectric constant.

With the pairwise-additive energy function, interactions between pairs of residues (and between residues and the protein backbone) in presence of the implicit solvent are computed before exploring the sequence-structure space, for each couple of possible amino acid types and each couple of rotamers. Values are stored in an energy matrix [fig. 1.10]. This calculation can be accomplished in parallel on a small cluster or on a multi-core desktop machine. During the exploration, energy updates are calculated by picking the necessary elements from the pre-computed matrix.

Pairwise residue GB interaction

We consider briefly a method to compute the many-body GB interaction accurately yet efficiently, so that the NEA can be avoided. Using the Still formula to approximate the GB interaction between a pair of charges, we compute the interaction between a couple of residues R and R ′ as

∆E int RR ′ = i∈R j∈R ′ τ q i q j r 2 ij + B R B R ′ exp -r 2 ij /4B R B R ′ 1/2 (1.35)
The right side of the equation has residue solvation radii B R and B R ′ instead of the usual atomic radii b i , b j . Radii B R can be calculated as a harmonic average over b i ∈ R,

∆E self R = i∈R ∆E self i = - τ 2 i∈R q 2 i b i def = - τ 2 ( i∈R q 2 i ) B R (1.36)
This formulation is called Residue-GB. As in atomic-GB, the interaction between two group of charges depends on the configuration of the whole solute through the residue solvation radii which are average of the atomic radii b i = b i (r 1 ,...,r N ).

Even if residue-GB is not pairwise in the interaction term, an alternative functional form was proposed by Archontis and Simonson [2005] and can be used to extrapolate the pairwise "information" about a couple of residues. Indeed, equation 1.35 can be fitted for a large range 4. The Proteus CPD framework of B values using the formula

∆E int RR ′ ({r i ,q i }; B = B R B R ′ ) ≃ c 1 ({r i ,q i }) + c 2 ({r i ,q i })B + +c 3 ({r i ,q i })B 2 + c 4 ({r i ,q i })B -1/2 + c 5 ({r i ,q i })B -3/2
For all pairs of residues in all posible rotamers whose coordinates and charges {r i ,q i } are known, 5 coefficients c i ({r i ,q i }) can be extracted and then inserted in the matrix. Then ∆E int can be directly computed during the MC simulation, keeping updated residues solvation radii from pairwise self-energies. More information can be found in [START_REF] Villa | Comparing pairwiseadditive and many-body generalized born models for acid/base calculations and protein design[END_REF]).

Monte Carlo exploration

Figure 1.11 -Amino acid mutation in MC. Two different types are represented using black and white circles. For each mutation in the folded state, the inverse mutation is performed in the unfolded protein, as described in the text. Without knowing the system partition function, we can compute a priori the relative probability to visit points in the configurational and sequence space. We generate a Markov chain of states using the well-known Metropolis acceptance method [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF]). In principle, the algorithm allows to sample any kind of probability density function. In our case, we look for states that are populated according to a Boltzmann distribution. Starting from an initial state (defined by a protein sequence and the set of occupied rotamers) we perform N iterations to generate N states of the system. At each iteration we perform trial moves to change conformation. These moves are accepted or rejected according to acceptance rules that exploit the detailed balance property of long, well-behaved Markov chains.

One possible Proteus trial move is a mutation. At a given simulation timestep, we modify the side chain at a random amino acid position i in the folded protein, mutating it to a random amino acid type t i → t ′ i picked from mutation space T i . A random rotamer is then selected from the rotamer space R i of the extracted type. At the same time, the reverse mutation is performed in an unfolded copy of the system, t which counts both folded and unfolded copies is calculated using equation 1.29. Another trial move changes the rotamer of a randomly picked residue r i → r ′ i without modification of its amino acid type.

Mutation and rotamer moves can be combined together to obtain a trial move that involves multiple positions in a single attempted transition. Other possible moves are the so-called mutmut or rot-rot moves where we change two side chains at two different positions. To increase the efficiency of these moves, we define lists of neighbours that are pre-computed before the simulation looking at pairwise interactions between residues. If two residues interact with an energy that is lower than a threshold value, they cannot be selected together in one of these "double" moves.

With the transitions and energy changes defined above, the desider Boltzmann distribution of states is determined by accepting or rejecting moves. Since moves are selected randomly they have an associaged "generation" probability α(o → n) where o and n indicates the old (current) and the new, generated states. We call acc(o → n) the corresponding move acceptance probability. The overall probability associated with the move is the product π = α(o → n)acc(o → n). At thermal equilibrium, the average number of leaving moves from a state o must be equal to the number of moves going from all other states n to the state o. In other words, the flux of configurations in one direction is canceled by the flux in the opposite direction. This is the detailed balance condition, which holds in the limit of a very long simulation. On average, the probability of finding the system in the two configurations corresponds to the equilibrium populations N (o) and N (n) which obey Boltzmann distributions. The detailed balance reads

N (o)π(o → n) = N (n)π(n → o) (1.37) so N (n) N (o) = π(o → n) π(n → o) = α(o → n)acc(o → n) α(n → o)acc(n → o) (1.38) Since N (n) N (o) = exp [-β∆E(o → n)] (1.39)
the acceptance probability must obey

acc(o → n) acc(n → o) = exp [-β∆E(o → n)] α(n → o) α(o → n) (1.40)
This is verified for the Metropolis scheme

acc(o → n) = exp [-β∆E(o → n)] α(n → o) α(o → n) if ∆E(o → n) > 0; 1 otherwise (1.41)
The main goal of CPD is to produce one or more protein variants with an associated sequence score in order to optimize a desired property. For example, one can use a scoring function based on the protein unfolding free energy in order to study protein stability. In general, one is not directly interested in the quantitative estimation of free energy differences but in producing a relatively small set of protein variants, which will be eventually studied with more sophisticated methods. For this reason, several scoring functions do not directly represent a physical quantity and are expressed in arbitrary units. However, the CPD model along with the Monte Carlo sampling described above can be used to estimate equilibrium thermodynamic quantities like binding free energy differences or protonation probabilities at constant pH. Proteus is particularly suitable for this kind of calculations: the fact that it is based on well-established physical models allows to define sampling methods which target different equilibrium properties.

Constant-activity Monte Carlo and ligand titration A special case is the design of ligand binding. We consider the simple case of a protein of fixed sequence S wich can bind two different ligands, represented by one mutating site of mutation space L and L ′ . ′ with protein residues are pre-calculated for all possible couples of protein rotamers. In this special case, the bound state is represented using the folded state energy function [eqn. 1.30] while the unbound state is represented using the unfolded state energy function [eqn. 1.31]. The ligand unfolded energies (in this case we call them unbound energies) have a logarithmic dependence on the concentration

E ub (X) = kT ln [X] + e ub (X) (1.42)
where X is equal to L or L ′ and kT is the thermal energy. The second term e ub (X) is the energy of the unbound ligand calculated ahead of time using the molecular mechanics energy function, including the implicit solvent. During Monte Carlo, mutations are performed swapping the ligand type between L and L ′ .

Figure 1.13 -Ligand mutation where two different ligands are represented with black and white circles. For a ligand mutation in the bound state, the inverse mutation is performed in the unbound state.

For two ligand states X and Y, moves are accepted according to the energy difference

∆E bound (X → Y ) -∆E ub (X → Y ) where ∆E ub (X → Y ) = kT ln [Y ] [X] + e ub (Y ) -e ub (X) (1.43)
Choosing [X], [Y ] ahead of time, we then perform a "constant activity" MC. Constant activity MC can be used to estimate binding free energy differences ∆∆G 0 (X → Y ). This method has been recently applied to compare the binding affinity between tyrosyl-tRNA synthetase (TyrRS) and several substrates [START_REF] Druart | Protein:ligand binding free energies: A stringent test for computational protein design[END_REF]). It consists in performing a series of simulations, each one at a different ∆E ref (X → Y ). From equation 1.43, this corresponds to varying the ratio of ligand concentrations [Y ] / [X]. From each simulation, one extracts the MC populations of X and Y , then standard binding free energy difference is computed from the titration midpoint [fig. 1.14] where X and Y are equally populated Constant-pH Monte Carlo As in constant-activity MC, the titration protocol can be applied to simulate protonation-deprotonation events at constant pH (Polydorides & Simonson [2013]). This allows to calculate pK a values of side chains of proteins, a quantity of great interest as it is connected to their stability, binding and many other biochemical properties.

∆∆G 0 bind (X → Y ) = kT ln [Y ] [X] mid = ∆E ref mid -∆e ub (X → Y ) (1.
Several precise and well-established methods to obtain protein pK a s are based on explicit and implicit solvent free energy calculations [START_REF] Simonson | Proton binding to proteins: pKaCalculations with explicit and implicit solvent models[END_REF]). However, CPD allows to screen larger proteins in a limited amount of time.

Following the Proteus protocol, titrable protein amino acid side chains are able to mutate in their mutation space. An energy matrix is computed, where pairwise interactions are obtained using protonated and deprotonated rotamers at all titrable positions. We indicate their generic deprotonated and protonated forms X and X + . For each couple of types X and X + , their relative population in solution (model compounds) can be controlled using the unfolded energy difference

∆E uf (X → X + ) = 2.303kT (pH -pK model a ) + e ub (X + ) -e ub (X) (1.45)
where pK model a is a tabulated quantity, known from experiments. The standard protonation free energy ∆G 0 depends on the dissociation constant at equilibrium through ∆G 0 = kT lnK a = -2.303kT pK a . In equation 1.45, the difference pH -pK a is used to compute the protonation free energy at a different pH value. The pK a shift ∆pK a is then obtained from the ratio of protonated/deprotonated populations in the protein extracted from a MC simulation:

∆∆G 0 = ∆G 0 protein -∆G 0 model = -2.303kT ∆pK a (1.46)
Chapter 2

High throughput design of protein-ligand binding

The computational design of ligand binding is the search for preferred sequences using a scoring function based on the binding free energy. However, most of CPD methods usually identify few candidates, looking for low energy states close to the global minimum energy configuration (GMEC). They often use empirical energy functions, without giving binding free energies to be compared with experimental affinities. As a consequence, most of CPD programs cannot be used for the more challenging high-throughput design, where one wants to optimize the binding between a receptor and many ligand variants in a short cpu time, obtaining formally exact binding free energies.

Here we present a general high-throughput design method that allows to overcome these limitations. It was implemented in our in-house program Proteus, but can be easily adapted to other CPD programs. Our protocol adaptively flattens the energy landscape in sequence space, and then extracts formally exact binding free energy differences from biased Monte Carlo (MC) simulations. Our test system was the Tiam1 PDZ domain (Tiam1), which binds to its natural peptide ligand Syndecan1 (Sdc1) composed of 8 amino acids. We used the method to design peptides that bind to Tiam1 and could serve as inhibitors of its activity. We targeted four of the last five C-terminal amino acids of Sdc1. Considering 18 possible amino acid types at the four "active" positions (we allow all amino types except Gly and Pro), there are 18 4 =104976 competiting ligands. Our protocol is efficient, since we obtained binding free energies for almost 75000 peptide variants using just one CPU hour. Method, simulations and results are presented in the attached paper (Villa et al. [2018b]).

Selected peptide variants were finally studied using a semi-empirical PB/LIE method, previously optimized to predict binding specificity of the Tiam1 PDZ domain. It was recently developed in our laboratory (Panel et al. [2017]). Details about PB/LIE calculations and results are given below, after the attached paper. At the end of this chapter, we also describe some theoretical aspects of the new design protocol, and how biased MC sampling techniques can be used to design sequences either by binding affinity or by specificity. For the high throughput design of protein:peptide binding, one must explore a vast space of amino acid sequences in search of low binding free energies. This complex problem is usually addressed with either simple heuristic scoring or expensive sequence enumeration schemes. Far more efficient than enumeration is a recent Monte Carlo approach that adaptively flattens the energy landscape in sequence space of the unbound peptide and provides formally exact binding free energy differences. The method allows the binding free energy to be used directly as the design criterion. We propose several improvements that allow still more efficient sampling and can address larger design problems. They include the use of Replica Exchange Monte Carlo and landscape flattening for both the unbound and bound peptides. We used the method to design peptides that bind to the PDZ domain of the Tiam1 signaling protein and could serve as inhibitors of its activity. Four peptide positions were allowed to mutate freely. Almost 75 000 peptide variants were processed in two simulations of 10 9 steps each that used 1 CPU hour on a desktop machine. 96% of the theoretical sequence space was sampled. The relative binding free energies agreed qualitatively with values from experiment. The sampled sequences agreed qualitatively with an experimental library of Tiam1-binding peptides. The main assumption limiting accuracy is the fixed backbone approximation, which could be alleviated in future work by using increased computational resources and multi-backbone designs. Published by AIP Publishing. https://doi.org/10.1063/1.5022249

I. INTRODUCTION

Computational protein design (CPD) is an emerging technique for engineering and understanding protein structure and function. 1-4 CPD explores a large space of amino acid sequences and conformations to identify protein variants that have predefined properties, such as ligand binding. Conformational space is usually defined by a discrete library of side chain rotamers and a small set of allowed backbone conformations. The energy function usually combines physical and empirical terms; 5-7 the solvent and the unfolded state of the protein are described implicitly. The space of sequences and conformations is often explored with a simulated annealing Monte Carlo (MC) approach. 1,8,9 Other approaches search for the lowest energy state, the "Global Minimum Energy Conformation" or GMEC, using combinatorial optimization. 10-14 To design ligand binding, the quantity to optimize is the binding free energy, which is not a property of the GMEC but a Boltzmann average over many states. This is a complex problem, and most CPD applications have simply tried to identify one or a few low energy states. Methods have been developed to systematically enumerate low energy states within a certain interval above the GMEC and to compute the configuration integral of the molecule(s) up to a predefined accuracy. They are referred to as "partition function" a) Electronic mail: thomas.simonson@polytechnique.fr methods. Enumeration methods are applicable to problems with limited numbers of degrees of freedom. They exploit the discrete nature of sequence and rotamer space and use efficient methods from computer science to explore them. Thus, cost function network and dead end elimination approaches have been developed that provide deterministic guarantees on configuration integral accuracy. 15-17 The binding free energy is then obtained as the ratio of bound and unbound configuration integrals. Related methods have been developed to enumerate secondary structures of simplified RNA molecules 18,19 and transmembrane proteins. 20 Of course, partition functions are not needed to obtain binding free energies. Molecular dynamics (MD) and Monte Carlo (MC) simulation methods have been developed for decades that focus directly on partition function ratios, instead. One explores states that are "intermediate" between the bound and unbound states, through importance sampling, using a bias energy to increase their statistical weights. 21-23 To compute binding free energy differences between two ligands, one alchemically transforms one ligand into the other. 24-26 These methods carefully avoid calculating partition functions yet provide excellent precision when comparing pairs of ligands if enough sampling is performed. They can be accelerated by employing an implicit model of the solvent. 27 MD and MC simulations have also been applied to more complex problems. For the competitive binding of many ligands to the same protein, 28-30 an adaptive, "lambda dynamics" method was proposed. The ligand chemical potentials are adaptively adjusted during MD until all possible complexes are similarly populated. Related adaptive MD methods have been used to sample conformational basins of proteins or peptides. An energy penalty is added to basins that have already been explored, leading to a gradual flattening of the conformational free energy surface. Metadynamics, 31-34 Wang-Landau methods, 35 and multicanonical simulations 36-38 are three main examples.

In protein design, the space to explore is especially complex since it includes amino acid sequence variations at multiple positions in a protein or peptide. Thus, adaptive simulations are an attractive perspective. A Wang-Landau MC approach was applied to RNA folding 39 and to protein backbone conformations. 40 A simpler, constant-activity MC protocol was used to study the competition between pairs of ligands binding to the tyrosyl-tRNA synthetase enzyme; 41,42 the relative binding free energies obtained for four small ligands were in qualitative agreement with experiment. Finally, an adaptive, Wang-Landau MC approach was applied to the binding of peptides to several PDZ protein domains. 43 A simulation of the unbound peptide was performed with several positions allowed to mutate freely; a bias potential was adaptively constructed so that all available sequences were equally populated. The same bias potential was then employed in a simulation of the protein:peptide complex. In this simulation, sequences are then populated according to their relative binding free energies, an elegant solution to the free energy design problem.

Here, we propose three improvements to the adaptive Wang-Landau MC approach 43 and apply it to a PDZ:peptide system of biological interest. The first improvement is the use of Replica Exchange MC (REMC) for improved sampling. 44,45 This allows us to explore a much larger mutation space than previously. Below, we simulate four mutating positions and over 10 5 competing ligand sequences using one CPU hour. The second improvement is the use of separate adaptive simulations and bias potentials for the bound and unbound states. By adaptively flattening sequence space for the protein:peptide complex, we make the method more robust and allow improved sampling of weak-binding peptides. Third, we have implemented the method within the established, Proteus CPD software. 41,46 This gives us access to improved, Generalized Born (GB) implicit solvent models, 47-49 efficient exploration using a precalculated energy matrix, a user-friendly interface, and additional Proteus features not applied here, such as multi-backbone design. 40 We consider the PDZ domain of the Tiam1 protein. The peptide ligand corresponds to the eight C-terminal residues of the Syndecan-1 protein. Tiam1 interacts functionally with Syndecan-1 by binding to its C-terminus through its own PDZ domain, helping to regulate downstream signaling pathways in several cell types. 50,51 From now on, we refer to the peptide as Sdc1 and to the PDZ domain as Tiam1. Our goal is to design peptides analogous to Sdc1 that can bind and inhibit Tiam1 either in vitro or in a cellular context. Such peptides could modulate Tiam1 function and have applications as reagents in cell biology or as inhibitors to downregulate Tiam1 in tumor cells. We redesign the Sdc1 peptide by allowing four of its last five residues to mutate (out of eight). This corresponds to a pool of 18 4 = 104 976 competing ligands (we allow all amino types except Gly and Pro, 18 in all). The protein and peptide backbone are held fixed, while side chains mutate and/or explore rotamers. The fixed backbone allows precalculation of an energy matrix, which makes the MC exploration very efficient. 41,52 The ligand pool is simulated twice, in the bound and unbound states, respectively. The ligand populations are gradually equalized by an adaptive Wang-Landau bias potential. They can then be reranked based on their binding free energies. Thus, our protocol directly designs peptides for their binding free energies.

Computed binding free energies were obtained for almost 75 000 peptide variants using one simulation of the peptide and one of the PDZ:peptide complex, each included one billion MC steps per REMC replica and required 1/2 h of CPU time on a desktop computer. The free energies were in fair agreement with experiment for seven variants, with a mean unsigned error of 0.8 kcal/mol (excluding one outlier). A sequence logo based on the MC sequences was qualitatively similar to the one based on a combinatorial library of peptide sequences selected experimentally for Tiam1 binding. 53 In future applications, longer simulations would allow larger numbers of mutating positions to be explored and more accurate variants of the Proteus solvent model to be used. We expect that using multiple conformations for the protein and/or peptide backbone, with the help of a hybrid MC method, 40 would also give improved accuracy. The ability to design sequences directly for binding affinity, without the need for expensive partition function calculations, should facilitate a number of biotechnology applications.

II. COMPUTATIONAL METHODS

A. MC exploration of side chain rotamers and types

We consider a polypeptide of L amino acids. Below, it will be either a folded protein:peptide complex or an unbound, extended peptide. We assume that each amino acid i can adopt a few different types s, s ′ , . . ., leading to different possible sequences S. The polypeptide backbone has a fixed geometry. The side chains can each explore a few discrete conformations r, r ′ , . . . called rotamers (around 10 per side chain type s). The energy of any particular conformation depends on the sequence and the particular set of rotamers. We perform a Monte Carlo exploration, whose goal is to generate a Markov chain of states, 55-57 such that the states are populated according to a Boltzmann distribution. The simulation system explicitly includes one copy of the polypeptide, whose sequence and rotamers can fluctuate. One possible MC move is to change a rotamer r i at one particular position; the energy change is ∆E on = E(. . .s i , r ′ i . . .) ☞ E(. . .s i , r i . . .), where the subscripts o and n refer to the old and new rotamer states. Another possible elementary move is a mutation: we modify the side chain type s i → s ′ i at a chosen position i in the folded protein, assigning a particular rotamer r ′ i to the new side chain. Let α(o → n) be the probability to select a move between two states o and n; let acc(o → n) be the probability to accept it. The overall probability of the move is π

(o → n) = α(o → n) acc(o → n).
The Metropolis-Hastings scheme 55-57 chooses the acceptance probability as

acc(o → n) = min(1, exp(-β∆E on )) α(o → n) α(n → o) , ( 1
)
where β is the usual inverse temperature. This, along with the assumption of detailed balance 45,56,57 leads to Boltzmann statistics,

N(n) N(o) = e -β∆E on , (2)
where N(o) is the equilibrium population of state o (respectively, n).

B. Adaptive bias energy

The adaptive simulation method is closely analogous to the Wang-Landau and metadynamics approaches. 32,35 For each position i among a selection of positions to be designed, we perform a long MC simulation where only i can mutate, and we gradually increase a bias potential until all the side chain types at i have roughly equal populations. The bias potential, E B i (r i , t), thus depends on time t and side chain type s i . In practice, an increment e b i (s i ; t) is added at a particular time t; then, we run the simulation for an interval T before adding a new bias increment. At the end of each interval T, the bias increment is only added to one side chain type s: the one populated at the end of the interval. The magnitude of the bias increment depends exponentially on the current bias value, which leads to a gradually decreasing increment. Specifically, the bias potential is defined as follows:

E b i (s; t) = n;nT <t e b i (s i (nT ); nT )δ s i (nT ),s , (3)
e b i (s(t); t) = e 0 exp -E b i (s(t); t)/E 0 , (4)
where e 0 and E 0 are constant energies and δ s,s ′ is the Kronecker delta. The above bias increment scheme is adapted from the well-tempered metadynamics methodology. 32-34 Once the single-position bias potentials E b i have been optimized separately, we are ready to perform simulations where all the selected positions can mutate simultaneously. We number these positions arbitrarily from 1 to p. We apply a bias that is the sum of the single-position biases,

E B (s 1 , s 2 , . . . , s p ; t) = p i=1 E b i (s i ; t). (5)

C. Relative binding free energy

The protein:peptide complex and the unbound peptide are simulated separately, each with its own bias potential. In the biased simulations, the populations of a particular sequence S in the bound and unbound states are denoted p b (S) and p ub (S). The binding free energy relative to a reference sequence S r , in the presence of the biases, has the form

∆∆G biased (S) = -kT ln p b (S) p b (S r ) + kT ln p ub (S) p ub (S r ) . ( 6
)
The relative binding free energy in the absence of bias has the form

∆∆G(S) = ∆∆G biased (S) -E B b (S) -E B ub (S) -E B b (S r ) + E B ub (S r ) . (7)
Subscripts b and ub indicate the bias energy in the bound or unbound state, respectively. In practice, Eq. ( 7) is applied to sequences that have been sufficiently sampled in both the bound and unbound states.

D. Protein design model and test system

The PDZ:Sdc1 complex is modeled using an experimental X-ray structure [Protein Data Bank (PDB) code 4GVD 58 ]. The backbone geometries of the protein and peptide are held fixed. Side chains explore a discrete rotamer library. 41,59 Selected peptide positions can mutate freely. The unbound peptide is modeled as a single, extended conformation, with a fixed backbone geometry. A molecular mechanics energy function is used along with a Generalized Born (GB) implicit solvent model. In principle, GB is a many-body model, such that the interaction energy between two amino acids I, J depends on the positions of the other amino acids. Here, a "Native Environment Approximation (NEA)" is made such that the I ☞ J interaction is computed assuming the rest of the system occupies its native sequence and structure. 41,47 With this approximation, all inter-residue interactions can be computed ahead of time and stored in an energy matrix. Later, they are used for sequence exploration. More details are given below.

We considered the X-ray structure of the Tiam1 PDZ domain bound to an octapeptide derived from the C-terminus of the Syndecan-1 (Sdc1) target protein (PDB entry 4GVD). We refer to the peptide as Sdc1. Amino acids in the peptide are numbered backwards from the C-terminus, as usual for PDZ ligands. The positions allowed to mutate here are positions ☞4 to ☞1. Position 0 is the C-terminal position. The N-terminus is acetylated and the C-terminus carboxylated. The peptide wildtype (WT) sequence is TKQEEFYA. Protonation states of histidines in Tiam1 were assigned to be neutral, based on visual inspection of hydrogen-bonding patterns in the 3D structure.

E. Effective energy function and energy matrix for MC

The energy calculations are done with a modified version of the Xplor program, 41,60 using the following effective energy function:

E = E bonds + E angles + E dihe + E impr + E vdw + E Coul + E solv . (8)
The first six terms in [START_REF] Hu | A gradient-directed Monte Carlo approach for protein design[END_REF] are taken from the Amber ff99SB molecular mechanics energy function. 61 The last term on the right, E solv , represents the contribution of a solvent. We use a "Generalized Born + Surface Area," or GBSA implicit solvent model, 62,63

E solv = E GB + E surf = 1 2 1 ǫ W - 1 ǫ P ij q i q j r 2 ij + a i a j e -r 2 ij /4a i a j -1/2 + i σ i A i . (9)
Here, ǫ W = 80 and ǫ P = 4 are the solvent and protein dielectric constants; r ij is the distance between atoms i, j and a i is the "solvation radius" of atom i. 62,63 A i is the exposed solvent accessible surface area of atom i; σ i is a parameter that reflects each atom's preference to be exposed or hidden from the solvent. The solute atoms are divided into 4 groups with the following σ i values [cal/(mol/Å 2 )]: unpolar (☞5), aromatic (☞12), polar (☞8), and ionic (☞9). Hydrogen atoms have a surface coefficient of zero. In the GB term, the solvation radius a i is a function of the coordinates of all the protein atoms. 62,63 We use a "Native Environment Approximation," where the solvation radii a i of a particular group are computed with the rest of the system having its native sequence and conformation. 41,47 The first stage in our CPD procedure 41,52,64 is to calculate the interaction energies between all side chain pairs and between side chains and the backbone, considering all allowed side chain types and the discrete rotamer library of Tuffery et al. (enhanced by allowing several orientations for SH and OH groups). 59,65 To alleviate the rotamer approximation, each side chain pair is initially positioned using two library rotamers; then, 15 steps of conjugate-gradient energy minimization are performed, where only the pair of interest can move and the energy is limited to the interaction within the pair, plus the pair's interaction with the backbone and solvent. 41,64 The final interaction energy is stored in an "energy matrix" for later use. 52 Thanks to the NEA, the GB interaction energy only depends on the types and rotamers of the pair of interest.

F. Replica exchange MC protocol

The MC moves included side chain mutations and rotamer changes at one or two positions. For two-position moves, the first position was chosen randomly and the second within its neighborhood, defined by an interaction energy threshold of 3 kcal/mol. Simulations were done with Replica Exchange MC, 45 with 4 replicas at thermal energies of kT = 0.6, 0.9,1.3, and 1.8 kcal/mol. A conformation swap between a random pair of replicas (at neighboring temperatures) was attempted every 2000 steps. In the adaptive MC calculations, the bias update period was T = 1000 steps. The bias increment amplitude was e 0 = 0.2 kcal/mol. The scaling energy E 0 was 29.4 kcal/mol. The simulations used to determine the single-position biases, E b i , lasted 10 8 steps (per replica). Once the biases were optimized, a production simulation was run with the total bias E B for 10 9 steps (per replica), both for the unbound and bound peptide. The positions allowed to mutate were positions ☞4 to ☞1 of the peptide (position 0 being the C-terminal position). Each position could adopt any side chain type except Gly or Pro, for a total of 18 possible types per position and 18 4 = 104 976 possible sequences overall. During the production simulation of the bound state, 101 159 sequences were sampled, and 74 963 (71%) were visited at least 1000 times during the PDZ:peptide simulation. For these, the relative binding free energies were estimated. The production simulations lasted about 1/2 h each on a desktop computer, thanks to the precomputed energy matrix and a shared memory, OpenMP parallelization.

III. RESULTS

The Sdc1 peptide was simulated in the unbound state, with its backbone fixed in the extended, β sheet geome-FIG. 1. The Tiam1-PDZ:Sdc1-peptide complex (cross-eyed stereo). Top: Wildtype complex, X-ray structure. Bottom: Complex with the TKQYTCTA peptide, which has the strongest computed binding affinity, representative MC structure. The peptide is yellow; its residues are labeled with their type.

try it has in the PDZ:peptide complex (Fig. 1). An adaptive bias potential E b i was constructed separately for each position i = ☞4, . . ., ☞1, such that when position i was the only position allowed to mutate, the distribution of side chain types was approximately flat. An overall bias potential was then defined as

E B (s , s 3 , s 2 , s 1 ) = i E b i (s i ).
The peptide was simulated for one billion steps with positions ☞4, . . ., ☞1 allowed to mutate, in the presence of the overall bias E B , using REMC with four replicas (see Sec. II). The C-terminal position 0 was not allowed to mutate, but kept its wildtype Ala side chain. All of the 18 4 = 104 976 possible sequences were sampled during the simulation. Table I reports the mean energy differences between each side chain type and the Ala type, taken as a reference, with the bias energy contribution removed. The values are averaged over the simulation. Values for positions ☞4, . . ., ☞1 are similar. These energies can be thought of as mutation energies for an extended, unfolded peptide chain. They contain contributions from intra-side chain interactions, side chain:backbone and side chain-side chain interactions, and side chain-solvent interactions (estimated with the GBSA implicit solvent model). Some of the energy differences are large and the corresponding types would not be sampled without using a bias potential E B to flatten the energy landscape.

A bias potential for the bound state was estimated in the same way, one position at a time (i = ☞4, . . ., ☞1). The complex was then simulated for one billion REMC steps with positions ☞4, . . ., ☞1 allowed to mutate, in the presence of the bias. 101 159 out of 104 976 possible sequences were sampled during the simulation; 74 963 (71%) were visited at least 1000 times. For these 74 963, we estimated the binding free energy difference ∆∆G, relative to Sdc1, by subtracting the mutation free energies in the bound and unbound states. The mutation free energies were estimated from the probability ratios between each sequence and the wildtype Sdc1 (see Sec. II). Populations were taken from the room temperature replica. Some of the sequences led to very unfavorable binding free energies, ∆∆G > 4 kcal/mol, because of steric clashes between a mutated peptide side chain and the PDZ protein.

Indeed, the protein and peptide backbones were held fixed in their X-ray conformation during the MC simulation and were unable to shift and make room for some of the larger side chain types. These large values can therefore be considered artefacts of the fixed backbone simulation. Figure 2 reports the ∆∆G values for the 50 000 sequences with the strongest binding, which span the range from ☞1.5 to 4 kcal/mol. For six sequences, reference values are available from experiment or, in two cases, from high-level, molecular dynamics (MD) free energy simulations. 66 The MD free energy simulations used explicit solvent and were found earlier to give good accuracy, within 1 kcal/mol of experiment. The reference values are also is TKQYTCTA, which differs from Sdc1 (TKQEEFYA) by four mutations. Figure 3 shows the MC sequences in the form of a sequence logo, where the sequences are Boltzmannweighted according to their MC binding affinities. An experimental logo is shown for comparison, corresponding to a combinatorial library of peptides, selected experimentally for Tiam1 binding. 53 Positions P 0 and P ☞2 are the most important for PDZ binding specificity. Position P 0 occupies three main types experimentally, C, A, and F, but was held fixed during the simulations. Of the four positions allowed to mutate, P ☞1 , P ☞3 , and P ☞4 are highly variable in both the MC and the experimental logos. Of the top ten MC types at these positions, 7 or 8 are present in the experimental logo, and vice versa, with somewhat different occupancies. Position P ☞2 is more conserved, both experimentally and in the simulations. Of the top four experimental types, Y, F, M, and T, all but T are in the top five MC types. While T has less than 1% occupancy in the MC sequences, the chemically similar types A, C, and S are highly populated. Overall, the MC logo is in reasonable agreement with the experimental one. The top scoring sequence TKQYTCTA is a candidate for enhanced Tiam1 binding.

IV. CONCLUDING DISCUSSION

To design ligand binding, the quantity to optimize is the binding free energy, which is a Boltzmann average over many states. Nearly all CPD applications to-date have simply tried to identify one or a few low energy states. In previous studies, we computed protein:ligand relative binding free energies within the CPD framework that were formally rigorous using constant-activity MC. 41,42 Here, we went further, designing sequences based on binding affinity in a high throughput manner. We used a Wang-Landau adaptive MC approach, proposed earlier for protein:peptide binding 43 and for backbone conformational free energies. 40 The method is analogous to the recent metadynamics and lambda dynamics methods. We modified the earlier protein:peptide binding method in three ways. First, we used a more powerful, Replica Exchange Monte Carlo sampling method. Second, we introduced two adaptive bias potentials instead of one: one for the unbound peptide and a new one for the protein:peptide complex. This leads to enhanced sampling of more weakly binding sequences. Third, we implemented the method in our Proteus software, which is a flexible, user-friendly, and freely available package for CPD that allows very fast Monte Carlo simulations. This will make the method easier to apply and to combine with various solvent models or force fields.

To perfectly flatten the energy landscape in amino acid sequence space, a complex bias energy would be needed. Here, we used a simple bias that is a sum of one-position terms, where each term depends only on the side chain type of a single amino acid. As a result, we achieve only a partial landscape flattening. This is not a difficulty since we can easily compute the relative binding free energy taking into account unequal type distributions; this is done by Eq. ( 6). Despite the imperfect landscape flattening, almost 75 000 peptide sequences were characterized in a single REMC simulation of the bound peptide and another of the unbound peptide, which required just 1/2 h each of CPU time, thanks to the precomputed energy matrix. For more complex applications, it would be straightforward to include selected side chain:side chain interactions in the bias potential to increase flattening.

Sequences were ranked by their binding free energies and those with large, positive ∆∆G values (about 1/3 of the sequences) were deemed artefacts of the fixed backbone approximation and discarded. For eight sequences, comparing to experimental or high quality MD free energy simulations gave a mean unsigned deviation of 0.8 kcal/mol, 1/3 of the experimental range and close to a simple Null model. The MC sequence logo was also in reasonable agreement with experiment. Among the predicted strong binders, the peptide variant E ☞3 D was characterized recently by extensive molecular dynamics simulations, 66 which predicted a weaker binding. This variant may thus represent a false positive of the design method (one out of eight variants tested). Additional predicted strong binders should be characterized experimentally to test for other false positives. Ideally, we would like to have experimental binding data for a full combinatorial library of peptide positions ☞5, . . ., ☞1, or an equivalent library for another PDZ domain. With decreasing costs for high throughput experiments, such data may become available in the future.

Among the discarded sequences, we expect there are some that would actually give strong binding, which can be considered false negatives of the fixed backbone design method. Similar problems were seen earlier 45 when designing four specificity positions within the Tiam1 PDZ domain. Depending on the backbone conformation of the peptide that was bound, large side chains were not always sampled at the designed protein positions. We expect that this artefact can be alleviated or eliminated by considering a few different peptide backbone conformations. The design calculations can either be done separately for each peptide backbone conformation or the conformations can be used together in a single, multibackbone MC simulation using a hybrid MC method proposed recently. 40 For a system of the same size as Tiam1, with the Proteus multi-backbone approach (which samples a Boltzmann distribution), we estimate it would take a few days to run 10 9 MC steps using a few 16-core machines. While much slower, the method remains feasible on a laboratory computer cluster.

Here, as a proof of principle, we used one CPU hour on a desktop machine to explore four mutating positions. In future applications, more resources could be used, allowing sampling to be increased by 2-3 orders of magnitude. This would allow us to explore larger numbers of mutating positions, with a multi-backbone model and a more rigorous Generalized Born model that eliminates the Native Environment Approximation. 49 The method can be applied to any protein:peptide or protein:ligand binding problem, allowing the binding free energy to be used directly as the design criterion. For the present, Tiam1:Sdc1 system, the top-scoring sequence TKQYTCTA (which is also the consensus sequence derived from the top 30 sequences) is a promising candidate for enhanced Tiam1 binding and experimental testing.

Chapter 2. High throughput design of protein-ligand binding

PB/LIE analysis

We computed binding free energies for 14 of the high scoring designed peptide sequences using PB/LIE. Sequences were extracted from the set of "top 30" designed binders (listed in Villa et al. [2018b], figure 3). We give first information about the PB/LIE method, from the published work of Panel et al. [2017].

Semi-empirical Free Energy Function

The binding free energy between Tiam1 and each designed peptide was estimated using the scoring function

∆G = α∆E vdW + β∆G P B + γ∆A + δ (2.1)
where α, β, γ are adjustable constants. The other terms ∆E vdW , ∆G elec and ∆A are differences between bound and unbound states, averaged over snapshots of MD simulations of the protein-peptide complex. More precisely, ∆E vdW is the average van der Walls interaction energy between protein and peptide. ∆G elec is the electrostatic free energy difference between bound and unbound states obtained with Poisson Boltzmann (PB) continuum electrostatics, and ∆A is the change in molecular surface upon binding. The weights α, β and γ were optimized (Panel et al. [2017]) to reproduce affinities for a target data set composed of 35 variants of the Tiam1:Sdc1 complex (including both peptide and protein mutations). Optimized parameters were α=0.02 (vdW), β=0.25 (elec) and γ=-4 cal/mol/Å 2 . After optimization the authours predicted experimental binding affinities for the target dataset, obtaining a mean unsigned error of 0.43 kcal/mol and Pearson correlation of 0.64 between experimental and calculate binding free energies.

Structural models and simulations setup

3D structures of the Tiam1:peptide complexes were obtained after reconstruction of the designed sequences with Proteus. For each peptide sequence under analysis, we extracted its rotamer configurations visited during the MC sampling in the bound state (hundreds for each sequence). Configurations were then sorted according to their Proteus energy (after subtracting the bias potential). The best rotamer configuration of each sequence was finally applied on the fixed backbone of the complex, and minimized with XPLOR. Histidines were set to be neutral, after inspection of hydrogen bonds in the 3D structures and pK a calculations performed with PropKa [START_REF] Bas | Very fast prediction and rationalization of pKa values for protein-ligand complexes[END_REF]).

PB/LIE analysis

Each reconstructed complex was then solvated in an octahedral box, filled with explicit solvent (TIP3P water). The system was then neutralized with sodium or chloride ions. After equilibration, production MD was run using NAMD 2.12. Simulations were done at room temperature and pressure with the Amber ff14SB force field [START_REF] Maier | ff14sb: Improving the accuracy of protein side chain and backbone parameters from ff99sb[END_REF]), using Langevin dynamics with a Langevin Piston Nosé-Hoover barostat (Feller et al. [1995]). Long range electrostatic interactions were computed using PME [START_REF] Darden | Computational biochemistry and biophysics, treatment of long-range forces and potential[END_REF]). Large energy fluctuations during dynamics were avoided using a flat-bottomed harmonic restraint of force constant 3.0 kcal/mol/Å 2 applied on the N-terminal peptide residue. Using such a restraint, the harmonic potential started to increase beyond a threshold distance of 3 Å between peptide and protein.

For each complex, production MD was run for 100 ns.

PB calculations

MD snapshots were extracted from trajectories (1000 frames for each complex). For each snapshot, water molecules were deleted and the system was described using continuum electrostatics. We used a protein dielectric constant ǫ p =8 and a solvent dielectric constant ǫ w =80. The electrostatic free energy difference was obtained taking the difference between the electrostatic free energy of the complex and of the peptide, with structures built from the same snapshot. The ∆G elec term was finally obtained averaging the electrostatic free energy differences over all the MD frames. PB calculations were run using CHARMM (Im et al. [1998]), solving the PB equation on a grid of 181 planes in each direction, with a 0.4 Å spacing between planes.

Results

Results of the PB/LIE analysis for the 14 selected peptides are resumed in table 2.1. We include, in the same table, the relative binding free energies obtained with Monte Carlo (from Villa et al. [2018b]). PB/LIE affinities were within 0.7 kcal/mol, while the ∆∆G obtained with MC were within 1.5 kcal/mol. We also note that with all MC ∆∆Gs were negative (they belong to the top 30 designed binders), while PB/LIE gived two variants (sequences ETMNA and EEMNA) that have essentially the same affinity of the WT, while all other PB/LIE variants have positive ∆∆G. It was recognized that the PB/LIE model has a small systematic error, overestimating the experimental binding affinities by 0.30 kcal/mol on average (Panel et al. [2017]). However, this cannot explain the observed discrepancies. Moreover, a real comparision is difficult because for the selected variants experimental data is unknown. It would be interestig to test some of these variants with more sophisticated, alchemical free energy Chapter 2. High throughput design of protein-ligand binding simulations [START_REF] Panel | Accurate PDZ/peptide binding specificity with additive and polarizable free energy simulations[END_REF]) which allow to extract free energies without any adjustable parameter.

Peptide

Bias convergence

The relative binding free energies in Villa et al. [2018b] were obtained by difference, simulating the ligand pool twice. Two separate bias potentials were constructed, running adaptive Monte Carlo in bound and unbound state. In a second round, the bias potentials were used in two separate production simulations. The final free energy differences ∆G bound and ∆G unbound were obtained removing the respective bias potentials from the flattened distributions.

The method is formally exact, and does not require perfect landscape flattening [fig. 2.1]. However, let us suppose that the bias converges perfectly, flattening the sequence space of the simulated system. In the long time limit, the free energy difference between any couple of 2. Bias convergence sequences goes to zero, and the bias potential converges to

E b (s,t → ∞) = - E 0 E 0 + kT F (s) + C(t → ∞) (2.2)
where kT =0.6 kcal/mol is the thermal energy, E 0 is the energy factor used to define the bias increment rule (equation 4 in the article) and C is a constant that does not depend on the sequence but grows as t → ∞. In equation 2.2, the quantity F(s) is the unbiased free energy of a sequence s at temperature T , that is completely compansated by the bias in the long time limit, up to a factor E 0 /(E 0 + kT ) and an arbitrary constant, as expected. A detailed derivation of equation 2.2 can be found in [START_REF] Barducci | Well-tempered metadynamics: A smoothly converging and tunable free-energy method[END_REF].

With an exponentially decreasing bias update rule (equation 4 in the article), the bias potential converges to the free energy F (s) of an ensemble at shifted temperature respect to the temperature T used in the adaptive MC. To compute free energy differences at temperature T , the bias potential must be multiplied by a factor that approaches 1 as kT is negligible respect to E 0 . For very large E 0 , one retrieves the result of standard metadynamics where 3 Perspectives: advanced sampling

E b (s) → -F (s) + C.
The bias potential defined in the article (equation 5) is made of "diagonal terms" or singleposition biases. As a consequence, it cannot fully compensate the different interactions between couples of residues of fluctuating amino acid types. To better flat the F (s) surface, one may consider a more complex (but still simple) bias potential, with a funcional form that includes double-position biases:

E b ij (s i (t),s j (t); t) = n;nT <t e b ij (s i (nT ),s j (nT ); nT )δ s i (t),s i (nT ) δ s j (t),s j (nT ) e b ij (s i (t),s j (t); t) = e 0 exp -E b ij (s i (t),s j (t); t)/E 0 (2.3)
Here e 0 and E 0 are adjustable parameters, that control the growth of the bias potential during the adaptive MC simulation. As for the single-position bias, the increment e b ij (s i ,s j ; t) is added at a particular MC step t only to one couple of side chains of type s i and s j . Such a doubleposition bias can be combined with the previously defined single-position bias (equation 3 and 4), to give:

E b (s 1 ,s 2 ,...,s p ; t) = p 1 E b i (s i ; t) + i<j E b ij (s i (t),s j (t); t) (2.4)
The Proteus implementation is under development, and it has been verified that this "combined" bias better flattens the sequence space in various situations. For example, it can be used to flat the free energy landscape of the unbound state of a protein:peptide complex. After convergence and supposing almost perfect flattening, the difference in bias potential for two couple of sequences can be used to directly obtain their free energy difference in the unbound state. Using equation 2.2, we write

∆G unbound (s → s r ) = E b (s r ,t ← ∞) -E b (s,t → ∞) E 0 + kT E 0 def = Ẽb unbound (s r ) -Ẽb unbound (s) (2.
5) The resulting, "scaled" bias potential Ẽb unbound (s) effectively flattens the free energy landscape of the unbound state at temperature T . It can be used to sample sequences of the bound state. The interesting result is that sequences sampled in this biased bound state will be populated according to their relative binding affinities:

∆∆G bind (s → s r ) = ∆G bound (s → s r )-∆G unbound (s → s r ) = ∆G bound (s → s r )-∆ Ẽb unbound (s → s r ) (2.
6) This represents a clear improvement respect to the former protocol described in Villa et al.

[2018b], since using the last equation one avoids to run the production MC simulations in the unbound state. Moreover, the design by affinity opens interesting potential applications for the engineering of PPI networks. Once a ligand pool has been flattened in the unbound state, the resulting bias can be applied to different complexes in separate simulations. This allows to screen binding affinities for many peptide variants binding to different receptors. In this case, one may also consider to compute an average Ẽb unbound (s) for many unbound state conformations, to better catch the peptide conformational changes upon binding.

Another interesting application is the design by specificity. It is based on the idea that flattening the sequence space for a given protein:peptide complex, the resulting potential Ẽb bound (s) can be used to sample sequences in the bound state of a different receptor. Indeed, as for the unbound state we can write

∆G bound (s → s r ) = E b (s r ,t ← ∞) -E b (s,t → ∞) E 0 + kT E 0 def = Ẽb bound (s r ) -Ẽb bound (s) (2.7)
Considering the relative binding free energies for a generic couple of receptors A and B

∆∆G A bind (s → s r ) = ∆G A bound (s → s r ) -∆G unbound (s → s r ) ∆∆G B bind (s → s r ) = ∆G B bound (s → s r ) -∆G unbound (s → s r ) (2.8)
and indicating with Ẽb bound,A (s) the bias potential that flattens the bound state for receptor A, one has

∆∆G A,B bind (s → s r ) = ∆∆G B bind (s → s r )-∆∆G A bind (s → s r ) = ∆G B bound (s → s r )-∆ Ẽb bound,A (s → s r )
(2.9) where the unbound state contributions cancel out. The interesting result is that sequences sampled in the bound state of B using Ẽb bound,A (s) will be populated according to their binding specificity for receptors B and A.

Chapter 3

Polarizable free energy simulations 1 Introduction

Molecular simulations need accurate energy functions to describe complex electrostatic interactions in crowded environments. Molecular dynamics (MD) needs also a simple, empirical representation of molecule's degrees of freedom. Most popular force fields like Amber (Maier [START_REF] Oostenbrink | A biomolecular force field based on the free enthalpy of hydration and solvation: the gromos force-field parameter sets 53a5 and 53a6[END_REF]) and OPLS [START_REF] Jorgensen | Development and testing of the opls all-atom force field on conformational energetics and properties of organic liquids[END_REF]) are carefully parametrized to describe both bonded and nonbonded interactions via a similar energy function. Their functional forms differ for few terms, usually related to internal degrees of freedom. They are parametrized following slightly different strategies, even if they use similar approaches to fit QM and experimental data. All these well-established force fields are able to reproduce both equilibrium states and thermal fluctuations of a large number of molecules, from small model compounds to big polymers. They are also able, in some situations, to fold small proteins in a reasonable amount of time (Lindorff-Larsen et al. [2011]). However, all these force fields share a common simplification in their description of electrostatic interactions.

et al. [2015]), Charmm (MacKerell et al. [1998]), Gromos
E nonbonded ({r}) = nonbonded pairs    ǫ min ij   R min ij r ij 12 -2 R min ij r ij 6   + q i q j ǫr ij    (3.1)
where {r} is the collection of distances r ij between nonbonded pairs of the system. The electrostatic interactions are computed as a sum over non-bonded pairs of the classical, Coulomb energy. Because of that, all these force fields are called additive. They have a limited represenation of the real molecular system because the charge of each atom i is condensed at one point in space r i and has a constant magnitude, neglecting induced polarization effects. This common feature makes all these force fields problematic when describing interactions between ionic groups in a heterogeneous environment, for example when charges are buried inside the Chapter 3. Polarizable free energy simulations protein core or are transferred from the bulk solvent to binding pockets. These force fields have also a limited accuracy in describing complex networks of salt bridges, when many ionic interactions fluctuate at thermodynamic equilibrium [START_REF] Debiec | Evaluating the strength of salt bridges: A comparison of current biomolecular force fields[END_REF]). In many applications this level of accuracy is acceptable, but there are cases where consequences can be dramatic.

Free energy calculations are among the most challenging applications in computational biochemistry. Despite the development of advanced sampling techniques and increasing computational power, there is still a great need of improved accuracy. Limitations are especially obvious when computing free energy differences for charging processes,

∆F (q 0 → q 1 ) q 0 = q 1 (3.2)
The precise calculation of electrostatic interactions for this kind of transormation is required when studying solvation, mutation free energies or ligand binding. For example, relative binding free energies for several point and double mutations in PDZ:peptide complexes were recently obtained from MD simulations, using a free energy perturbation protocol (Panel et al.

[2018]). Many peptide variants were studied using the Amber 99SB force field, computing their relative affinities and comparing them to experimental data. Good agreement was obtained for non-ionic mutations, giving a root mean square deviation of 0.37 kcal/mol, a mean unsigned error of 0.32 kcal/mol, a maximum error of 0.6 kcal/mol and a correlation coefficient of 0.84. The situation was quite different considering a set of ionic mutations, where the mutated side chains have a different charge. In this case errors were larger than 1 kcal/mol. This result is comparable to other studies in literature [START_REF] Deng | Calculation of standard binding free energies: Aromatic molecules in the t4 lysozyme l99a mutant[END_REF] ; [START_REF] Boyce | Predicting ligand binding affinity with alchemical free energy methods in a polar model binding site[END_REF]).

A new generation of molecular-mechanics force fields aims to overcome difficulties related to an additive, point-charge representation of electrostatic interactions. Because they include an empirical representation of electronic polarization, these polarizable force fields are more complex both in their functional form and in their parametrization strategies. Since they model the reponse of a molecule's electron charge distribution in different environments, their tuning requires extensive testing, carried out over several years of development. In the field of biomolecular simulation there are still many applications where polarizable MD has not yet been used. There are few examples of polarizable force fields that can be used routinelly to perform MD simulations at the 100 -1000 nanosecond timescale for free energy calculations, where one has to sample many microstates of systems composed of tons of thousands of atoms. The community recognises (Lemkul et al. [2016]) that there are actually three main models, implemented in different force fields, that have been shown to be more suitable and reliable for biomolecules: induced point-dipoles, fluctuating charges and the Drude oscillator model. The first is the one used in AMOEBA [START_REF] Shi | Polarizable atomic multipole-based amoeba force field for proteins[END_REF]), the second in the Charmm Fluctuating

Induced dipole model

Charge Force Field [START_REF] Patel | CHARMM fluctuating charge force field for proteins: I parameterization and application to bulk organic liquid simulations[END_REF] ; [START_REF] Patel | CHARMM fluctuating charge force field for proteins: II protein/solvent properties from molecular dynamics simulations using a nonadditive electrostatic model[END_REF]) and the latter in the Drude force field (Lamoureux & Roux [2003]). The present thesis focuses on Drude, even if a short description of AMOEBA and Charmm Fluctuating Charge force fields will be given. This chapter starts with a theoretical introduction about classical models used to study interaction between induced dipoles in molecular systems. The point-dipole model and key quantities will be presented. After a brief discussion of existing polarizable force fields suitable for molecular simulation, we will describe the Drude Force field. Several important concepts for the calculation of binding free energies will be discussed, with particular attention to charging processes. Since free energy differences are extracted from molecular simulations performed in infinite, periodic simulation boxes, several artefacts need to be taken into account.

2 Induced dipole model

Fields and dipoles

Modeling induced electronic polarization with an empirical, classical model dates back to 1972. The original work of Applequist et al. [1972] was based on the older study from [START_REF] Silberstein | Molecular refrectivity and atomic interaction[END_REF] that introduced the point-dipole model for polyatomic molecules. Atomic polarizabilities are assigned to each atom. Considering a system composed of N atomic dipoles μj located at positions rj , the electric field at position ri is

Ēi = - N j=1 T ij μj (3.3)
where T ij is the dipole field tensor for atomic centers i and j, which is a 3x3 matrix. If r is the distance between atoms i and j and where x = x i -x j , y = y i -y j and z = z i -z j , then

T ij = - 3 r 5      x 2 -1/3r 2 xy xz yx y 2 -1/3r 2 yz zx zy z 2 -1/3r 2      (3.4)
The induced dipole for atom i at position ri is computed using

μi = α i Ēi = α i   Ēext i - N j=1;j =i T ij μj   (3.5)
where Ēext i is an external electric field at position r i and α i is the polarizability tensor

α i =      α xx α xy α xz α yx α yy α yz α zx α zy α zz      (3.6)
For a molecule composed of N atoms in the presence of an applied electric field Ē, the induced molecular dipole moment μmol can be expressed as function of the molecular polarizability

α mol μmol = N i=1 μi =   N i=1 N j=1 B ij   Ē = α mol Ē α mol = N i=1 N j=1 B ij (3.7)
The matrix B ij contains atomic polarizabilities and the dipole field tensor. Its mathematical expression depends on the functional form of T ij and the interaction is parametrized by molecular polarizabilities. Indeed, equation 3.5 can be seen as a system of equations for each μi where i = 1, ..., N . In matrix form, it corresponds to

           α -1 1 T 12 ... T 1N T 21 α -1 2 ... T 2N ... ... T N 1 ... ... α -1 N                       μ1 ... ... ... μN            = Aμ =            Ē1 ... ... ... ĒN            = Ẽ (3.8) 2. Induced dipole model if B = A -1 then μ = B Ẽ, so that for each atom μi = N j=1 B ij Ēj (3.9)
This simple model proposed by Applequist makes good predictions of mean, molecular polarizabilities. However, this point-like description suffers from several problems. First, it produces a large anisotropy in molecular systems. Second, in some situations it leads to infinite polarizabilities. This problem arises when the polarizable atoms are too close. It can be avoided by modification of the dipole field tensor. Thole [1981] proposed an alternative interaction model in order to damp interactions between dipoles when they get closer than a threshold value. More precisely, he proposed to compute dipole-dipole interactions using a slightly different dipole field tensor, which is built by modification of charge distributions according to a shape function. This shape function, which can be linear or exponential, is parametrized by atomic polarizabilities and a general screening parameter τ :

ρ = τ 3 8π exp - τ r ij (α i α j ) 1/6 (3.10)
Using the threshold distance r * ij = τ (α i α j ) 1/6 the modified dipole filed tensor is then computed using

T (u = r ij /r * ij ) = - 3 r ij (r * ) 4      x 2 (1 -4 3u ) xy xz yx y 2 (1 -4 3u ) yz zx zy z 2 (1 -4 3u )      (3.11) when r ij < r * ij 2.

Electrostatic energy

The energy of a system of point dipoles is composed of several contributions. Its general prototype is

E elec ({r}) = E q-q + E q-µ + E µ-µ + E pol (3.12)
The first terms are charge-charge (q-q), charge-dipole (q-µ) and dipole-dipole (µ -µ) interactions. E pol is the work of polarization, also called "self polarization energy". Electronic degrees of freedom, parametrized by atomic polarizabilities α i , are represented using the smeared charge distribution. The dipole field-tensor T ij is the mathematical object used to calculate Chapter 3. Polarizable free energy simulations non-additive interactions (q-µ) and (µ -µ). It is function of the charge distributions of atoms i and j which depend, following equation 3.10, on the distance between centers i and j. A more explicit expression of equation 3.12 is

E elec ({r}) = i<j   q i q j r ij -q i α=x,y,z T α ij µ jα + α=x,y,z µ iα T α ij q j - α,β=x,y,z µ iα T αβ ij µ jβ   + U pol (3.13)
where T ij is the interaction tensor, such that

T α ij def = ∇ α T ij = - r ij,α r 3 ij T α,β ij def = ∇ α ∇ β T ij = 3r ij,α r ij,β -r 2 ij δ α,β r -5 ij (3.14)
We immediately see that, respect to the simple charge-charge Coulomb interaction, additional terms appear in the sum. These terms are parametrized by molecular polarizability and Thole factors, which control strength of induced dipoles and the shape of the charge distribution. The last term is the work of polarization, and can be derived considering the induced dipole μi of atom i in presence of electric field Ēi at point ri [START_REF] Böttcher | Theory of electric polarization[END_REF]). The energy of the induced dipole is

U = -μ i • Ēi + U pol (3.15)
At equilibrium, for small changes dμ i of the induced dipole moment, one has dU=0. So we can write

dU pol = -d -μ i • Ēi = Ēi • dμ i = μi α i • dμ i = 1 2α i d μ2 i (3.16)
The work of polarization is then obtained easily:

U pol = 1 2α i μi 0 d μ2 i = 1 2α i μ2 i (3.17)
For the system of dipoles μi , the total work of polarization or total "self-polarization energy" is finally

U pol = 1 2 i μi • μi α i (3.18)

Induced dipole force fields

AMOEBA (Ren & [START_REF] Bibliography Ren | Polarizable atomic multipole water model for molecular mechanics simulation[END_REF]) belongs to the class of force fields that make use of equation 3.13, as well as Amber ff02 (Cieplak et al. [2001]). AMOEBA includes also permanent, traceless quadrupoles, which are 3x3 tensors. So for each atom i it associates properties of the point-dipole model (charge q and dipole moment μi ) but also additional, quadrupole Q αβ components. All atoms are polarizable, including hydrogens. The direct use of the induceddipole representation is considered useful especially for its ability to reproduce anisotropy effects [START_REF] Baker | Polarizable force fields for molecular dynamics simulations of biomolecules[END_REF]). However, increased computational cost is demanded to compute interactions, as well as more-expensive force evaluations in MD. To give an idea, one can have a look at recent timings of all-atom, explicit solvent molecular dynamics simulations performed using AMOEBA, implemented on an optimized version of Tinker (Lagardère et al. [2018]).

Calculations were run on the Occigen supercomputer, hosted at the CINES centre of the French supercomputing agency GENCI. For plain MD of one solvated protein (dihydrofolate reductase, DHFR) composed of 23558 atoms the performances are of 7.2 ns/day on 960 cores. Recent calculations on a similar system performed in our laboratory using 224 cores using NAMD [START_REF] Phillips | Scalable molecular dynamics with NAMD[END_REF]) and an additive Amberff99SB force field show performances of ∼ 74 ns/day, over 40 times faster. Despite the strong effort in software development and the accuracy in the calculation of electrostatic interactions, long simulations at timescales of the order of 1 µs with AMOEBA have an enormous computational cost.

3 Point charge models

Fluctuating charge model

Given the high computational cost demanded by induced-dipole representations, a lot of effort has been made to build more empirical but less expensive energy functions suitable for MD.

In the fluctuating charge model, charges are considered as dynamical variables in an extended lagrangian framework. It is based on the principle of electronegativity equalization. Given a molecular system, charges fluctuate according to conformational changes that alter the local value of the electrostatic potential. Each atom has an associated electronegativity and within a molecule the system is constrained to keep the same electronegativity for each atom. In pratice, each charge has an associated mass M Q that is small compared to the real, atomic mass. In this way, electronic degrees of freedom (which evolve in the extended lagrangian) are able to relax quickly in reponse to conformational changes. The main advantage of the method relies is that the energy function used to compute interactions is unchanged compared to the standard, additive force fields. However, it suffer from a strong limitation: intramolecular electronic charge distributions are exchanged between atoms, so that charge transfer occurs mostly along bond lines. This is a strong limitation when one wants to model polarization in the direction perpendicular to the molecule plane.

3.2 Drude pseudo-particle model

q D q A E Figure 3
.1 -Drude system composed of a polarizable atom q A and its Drude particle q D , displaced according to an external field

In the Drude pseudo-particle model, the electronic cloud surrounding a polarizable, heavy atom is represented by a negative charge q D attached to the nucleus with a harmonic oscillator of force constant k D [fig. 3.1]. All atoms except hydrogens are by default polarizable. The harmonic spring has the same force constant for all atoms. The total charge Q of the atom is sum of the nucleus and Drude charges. Q is split between the nucleus A and its Drude particle D depending on the atomic polarizability α:

Q = q A + q D α = q 2 D k D M = m A + m D m D = 0.4 a.m.u. (3.19)
where the mass is expressed in atomic mass units. In general, each atom has an isotropic polarizability α xx = α yy = α zz . In the absence of an external field, each Drude particle oscillates around its reference nucleus. Otherwise it oscillates around a displaced position, defined by the vector

d = q D Ē k D (3.20)
which depends on the Drude particle charge, the force constant and the applied field. From the formula we see that for a given Ē, a bigger atomic polarizability α (a more negative charge q D ) produces a bigger displacement. The atom and Drude particle behave like a point charge Q of polarizability α, in the spirit of the point-dipole model.

Interactions between atoms and Drude oscillators that are separated by three or more bonds are computed using the standard, Coulomb energy function, while dipole-dipole interactions between 1-2 and 1-3 atom pairs are screened with the previously discussed method proposed by Thole. In the Drude force field, instead of using a general screening parameter τ for all atoms, interactions between dipoles i and j are screened by multiplying their interactions by

S(r ij ) = 1 -1 + τ ij r ij 2(α i α j ) 1/6 exp - τ ij r ij (α i α j ) 1/6 τ ij = τ i + τ j (3.21)
where τ i and τ j are screening (or Thole) parameters for atoms i and j. This formula corresponds to a smeared charge distribution

ρ = τ 3 ij 8π exp - τ ij r ij (α i α j ) 1/6 (3.22)
The potential function includes also a self-polarization energy

U self ( d) = 1 2 k D d2 (3.23)
The 2013 version of the Drude force field includes also anisotropic polarizabilities, where α xx = α yy = α zz . In this case, the self-polarization energy is calculated using the more general formula

U self ( d) = 1 2 k xx d 2 x + k yy d 2 y + k zz d 2 z (3.24)
Drude has also other auxiliary particles, called lone pairs [fig. 3.2], not present in the first release of the force field. Lone pairs represent electron pairs in valence shells of atoms that do not participate in bonding. Their role is to include a charge density that is delocalized with respect to their reference atoms, using virtual rigid interaction sites. In fact, lone particles are massless and their forces are transferred to their reference atoms during MD simulations. The introduction of lone pairs in Drude improved agreement with QM calculations, both for energetics and for optimized geometries (Harder et al. [2006]).

Drude polarizable water models

The first Drude water model was published in 2003 (Lamoureux et al. [2003]) and named SWM4 : Simple Water Model with 4 interaction sites. In a former parametrization, it was called SWM4-DP because the oxygen carried a positive Drude charge. Along with the usual oxygen and hydrogen atoms, a rigid particle called "LPW" was attached to the oxygen in order to better reproduce the gas-phase dipole moment and quadrupole moments. In analogy with the TIP4P water model, this lone charge is negative and its magnitude is twice the charge assigned to hydrogen atoms.

To be more consistent with the physical model, where negative Drude charges represent electronic degrees of freddom, a new model called SWM4-NDP [START_REF] Lamoureux | A polarizable model of water for molecular dynamics simulations of biomolecules[END_REF], Negative Drude Polarization) was introduced. This is the "default" Drude water in the current version of the force field. When compared to the non-polarizable and widely used TIP3P water, SWM4-NDP better reproduces several thermodynamic properties, especially in situations where a reponse to different environments is needed. In particular, the molecular dipole moment increases from the gas phase to a hydrogen-bonding environment. Another model, called SWM6-NDP (Yu et al. [2013]) adds two lone pairs out of the molecular plane. It reproduces water proporties even better than SWM4, but its simulation is more expensive. An important feature of polarizable water is not just related to average properties, but also to dynamics. Polarizable models reproduce water self-diffusion and viscosity better than fixed-charge models and this should give more realistic timescales in molecular dynamics simulations. However, more interaction sites need more parametrization, especially in charged enviroments. One

Point charge models

example is the case of divalent ions, in particular the interaction between SWM4-NDP and Mg 2+ (Lemkul & MacKerell [2016a]).

Simulating the Drude force field

To carry out consistent polarizable MD simulations, electronic degrees of freedom { d} need to correctly relax depending on the nuclei positions {r} in order to reproduce the response to variations of conformation-dependent electric fields. In the original Drude model electrons were represented by oscillating massless particles. For each oscillator there is a unique displacement d⋆ i such that the energy is minimized:

∂U ∂ di = ∂U self ({ d}) ∂ di + ∂U int ({ d,r}) ∂ di = 0 (3.25)
This is referred to as the Self Consistent Field condition, or SCF condition. The SCF condition can be achieved by minimizing the Drude positions at each MD step. This is the most straightforward solution to produce conistent trajectories, and is compatible with massless Drude particles. However, the procedure needs accurate implementation, is computationally expensive and can be affected by numerical errors.

The energy minimization can be avoided by following an alternative procedure. In the current Drude force field, Drude particles are not massless, so their motion can be integrated as the rest of the system. One "physical" temperature T is used to thermalize the real atoms while Drude oscillators are kept at a low temperature T D ≪ T using a dual Nosé-Hoover thermostat. They evolve almost adiabatically in response to nuclear conformational changes which modify the local field. The SCF condition is approximated, following Lamoureux and Roux, using the expansion

U ({r, d}) = U SCF ({r}) + i δ di ∂U ∂ di (d ⋆ i ) + 1 2 ij δ di ∂U ∂ di (d ⋆ i ) ∂U ∂ dj (d ⋆ j )δ dj + ... (3.26)
where d ⋆ i are the relaxed displacements for which the SCF condition is satisfied, and δ di = di -d ⋆ i .

MD implementation

Drude systems are preparated using CHARMM, which implements the routines necessary to transform a system described with an additive force field to a Drude system. Drude topology files have additional keywords: we give, as an example, the SWM4-NDP water topology As in the additive force fields, the residue is declared using the RESI keyword followed by its name and its charge. The optional keyword GROUP is used to indicate that the following atoms have a total charge equal to zero. Then the residue's atoms are listed, each one using the ATOM keyword followed by its name, its atom type and its charge. Polarizable atoms (in this case the oxygen OH2) have additional keywords, used to indicate the value of their polarizability ALPHA and their THOLE screening parameter. If THOLE is missing, a default value of 1.3 is used. Drude particles attached to polarizable atoms are not explicitly included in topology files, as they are automatically created by CHARMM. Their name is assigned using the letter D followed by the name of its reference atom. For example, OH2 carries a Drude particle of atom name "DOH2". CHARMM automatically assigns the default atom type "DRUD" to all oscillators, except in special cases where a different type is indicated using the optional keyword "TYPE". In SWM4, the oxygen Drude particle has such a special atom type, called "DOH2". This particular choice allow to specify special nonbonded interaction parameters (using NBFIX), explained later in the text. Virtual interaction sites of default atom type LP are indicated in the residue topology. Three atoms are used to define the particle position through the LONEPAIR line. In the example, coordinates of the OM lone pair particle are computed from positions of the OH2-H1-H2 atoms. Like the DOH2 oscillator, OM has its own (non-default) atom type called "LPDW". As in the additive water model TIP3P, SWM4 has an H1-H2 bond used to fix bond lengths and angles during dynamics.

Another Drude keyword called ANISOTROPY is used to specify magnitudes of diagonal elements α xx , α yy , α zz in the polarizability matrix. Since anisotropy depends on the molecule internal geometry, (x,y,z) are defined in a reference frame constructed using positions of three Using the equation A11+A22+A33=3, matrix elements α xx , α yy and α zz are set scaling AL-PHA according to the ratios between A11, A12 and A13.

After reading the Drude topology and parameter files with CHARMM, then the system is built using the CHARMM generate command generate @segname first @nterpatch last @cterpatch setup warn drude dmass 0.4

The sintax is similar to the one used to construct systems with additive force fields, the only difference is a drude keyword that creates vectors to manage Drude particles and LPs. The dmass m D = 0.4 a.m.u. is subtracted from polarizable atomic centers and assigned to their Drude particle. Charges are split between atomic centers q A and oscillators q D , using the ratio between polarizability and force constant, according to equation 3.19.

Once the coordinates for atomic centers are read in, Drude and Lone particle positions are assigned using the CHARMM commands coor sdrude and coor shake. Displacements of all oscillators are initially set to zero, requiring minimization before dynamics. This minimization is usually performed first on Drude particles with restrained atomic centers, then relaxing the whole system without restraints. Parameter files are very similar to the ones of the additive force field C36. Internal parameters are speficied in the BOND, ANGLE, DIHEDRAL sections and LJ parameters are specified after the NONBONDED keyword. As in C36, NBFIX is used to modify LJ interactions between two particular atom types i and j, using a special value of r min and e min . A new keyword called NBTHOLE (in the spirit of NBFIX) is used to specify τ ij for atom types i and j. It is usually introduced to balance particular interactions of charged groups, for example between nucleic acids and divalent ions.

The Drude force field is currently implemented in CHARMM, NAMD, OpenMM and GROMACS, even if it is only partially supported by the latter. NVT and NPT simulations can be run with CHARMM and NAMD exploiting the dual thermostat, while GROMACS completely supports just the NVT ensemble. Its implementation for the NPT ensemble is less efficient because it uses the older "SCF minimization" method to relax electronic degrees of freedom. Moreover, GROMACS does not support NBTHOLE parameters. In NAMD, pressure is controled using Langevin dynamics to control fluctuations in the barostat [START_REF] Martyna | Constant pressure molecular dynamics algorithms[END_REF] ; Feller et al. [1995]). Details about extended lagrangian implementation can be found in Jiang et al. [2010]. It is clear that to simulate a system using Drude, extra computing time is demanded. In terms of CPU time, the Drude force field is almost four times more expensive than a non-polarizable force field like Charmm C36. A factor of two is introduced by the many auxiliary particles. SWM4 has four interaction sites, all heavy atoms have an additional Drude particle and several atoms carry lone pairs. In addition, harmonic forces (with the relative self-polarization energies) need to be calculated at each timestep for all electronic degrees of freedom. Another factor of two comes from the dual thermostat: equations of motions need to be integrated with a small timestep, which should never be greater than 1 fs.

As an example, in table 3.1 we compare MD performances of one additive force field (C36) and the Drude force field. We considered a protein:peptide complex, the Tiam1 PDZ domain (89 amino acids) with its natural peptide binder Syndecan1 (8 amino acids). The system was solvated in an octahedral box containing 10139 water molecules. We used TIP3P water for C36 simulations, SWM4-NDP water for Drude simulations. As expected, the number of atoms of each system was proportional to the number of water molecules. TIP3P is composed of three atoms, where SWM4-NDP has four interaction sites plus a Drude particle on the oxygen atom. As a consequence, the C36 system was composed of roughly 30k atoms and the Drude system was composed of roughly 50k atoms (including Drude and lone particles). MD was done using NAMD 2.12. Long-range electrostatic interactions were computed with PME, LJ interactions were shifted to zero for separations between 10 and 12 Å. We used a timestep of 2 fs for C36 and of 1 fs for Drude. Simulations were run on the Occigen supercomputer of the CINES center of the French computer agency GENCI. For both simulations we used 224 cores, distribuited over 8 nodes equipped with 28 processors. As expected, C36 was four time faster than Drude. Using the extended langrangian which exploit the dual-thermostat protocol described above, polarizable simulations at timescales up to microseconds can be run using massive parallel architectures. However, some considerations about the stability of dynamics are necessary. Since the SCF condition is satisfied approximately, there are situations where atoms are in very polar environments and the Drude particle may oscillate too far from the corresponding nucleus. simulation where the Drude oscillator follows a normal regime, where the distance between q A and q D is always lower than 0.2 Å. Right: simulation where the Drude oscillator tends to cross the threshold value of 0.2 Å. A DrudeHardWall option is applied at frames shown in red. In a stable system, Drude particles usually oscillate around their reference nucleus at average distance of 0.1 Å [fig. 3.3, left]. However, large forces, arising from strong electrostatic interactions, sometimes lead to the so-called polarization catastrophe. In the worst case, the system will be so unstable that the dynamics will be broken and the MD engine will stop the simulation. The polarization catastrophe can be avoided using a smaller timestep to integrate equation of motions, decreasing it to 0.75 fs or less. However, because this represents a further increase of CPU time, an alternative solution called DrudeHardWall has been introduced. When a Drude particle exceeds a threshold distance respect to its center, by default equal to 0.2 Å, an hard-wall potential is applied [fig. 3.3, right] and the Drude particle is reprojected along its bond direction. The negative displacement is generated using the projection of its velocity along the bond line, flipped and rescaled according to the oscillator temperature T D . This should be sufficient to correct the trajectory and restore the correct regime. Since DrudeHardwall can be considered as a workaround, it should not be applied for large portions of trajectories, as in the example shown in figure 3.3. Because of its implementation, it leads to unrealistic frames and wrong energetics. The force field has been extensively tested and this problem should not occour too often. However, this is not always true for unofficial, in-house parametrized molecules. To avoid this situation, one may increase screening of electrostatic interactions by adjustment of Thole factors τ i of particular atoms.

Description

Considering a receptor R and a ligand L at equilibrium, their binding is governed by the equation

RL ⇆ R + L (3.27)
The binding constant (dimensions of a volume) is

K eq = ρ RL ρ R ρ L = [V ] (3.28)
It depends on the concentrations of the complex RL and of the two unbound species R and L. We note that the use of equation 3.28 needs a definition of bound and unbound states.

We want to connect the binding constant (which can be obtained from experiments) to the corresponding binding free energy, more precisely to the free energy that can be obtained from a computer simulation. We use X to indicate the different species in the solution: the ligand L, the receptor R or the complex RL. From the Gibbs free energy change

dG = V dP -SdT + X µ X dN X (3.29)
at the conditions of constant temperature and pressure (dT = 0 and dP = 0) we obtain the change in free energy upon binding of one molecule of R and one molecule of L, to form one molecule of the complex RL

∆G = µ RL -µ R -µ L (3.30)
µ X is the chemical potential of species X,

µ X = ∂G ∂N X N Y =X ,P,T = -kT ∂lnQ ∂N X N Y =X ,P,T (3.31)
For a non-ionic, dilute solution it has a logarithmic dependence on concentration µ X (P,T ) = µ 0 X (P,T ) + kT ln

ρ X ρ 0 X (3.32)
4. Standard and relative binding free energies µ 0 X (P,T ) depends on P and T and is the standard chemical potential. It is defined by an arbitrary standard state concentration ρ 0 X = ρ 0 for all chemical species X, so that

∆G = µ 0 RL -µ 0 R -µ 0 L + kT ρ RL ρ 0 -kT ρ R ρ 0 -kT ρ L ρ 0 = ∆G 0 + kT ln ρ 0 ρ RL ρ R ρ L (3.33)
We note that the right term in the formula depends on K eq . The formula holds for a hypothetical equilibrium state where concentrations are constrained to be constant; we call it reaction free energy. Without constraints, the system goes to a state where ∆G = 0 and so

∆G 0 = -kT ln ρ 0 ρ RL eq ρ L eq ρ R eq = -kT ln ρ 0 K eq (p,T ) (3.34)
The standard chemical potential of equation 3.32 can be computed from the ratio of two partition functions. Indeed, we can write

µ 0 X = -kT ln Z N,X (V N,X ) Z N,0 (V N,0 )V N,X ρ 0 + P 0 VX (3.35)
where Z N,X is the canonical partition function of a system containing one molecule of solute X solvated in N solvent molecules at volume V N,X . Z N,0 is the partition function of the same system without X at volume V N,0 . The two volumes are the equilibrium volumes of the two systems (with or without solute) when they are mantained at pressure P 0 . V = V N,X -V N,0 is the change in equilibrium volume when a molecule of solute X is added to the system with N solvent molecules, or partial molar volume.

In equation 3.34, the free energy has a logarithmic dependence on the binding constant K eq . We also note that, in the equation 3.28, K eq has dimensions of a volume and the logarithm argument in equation 3.34 is dimensionless. This trivial argument based on dimensional analysis (General [2010]) is foundamental to note that computer simulations and experiments may divide k eq by two different concentration factors. The choice of standard state concentration is arbitrary, it is usually set to be 1M. It is important, when comparing K eq obtained from an experiment and the one obtained from a simulation, that both equilibrium constants (and so free energies) refer to the same standard state. The relationship between the free energy difference ∆G 0 and the one measured in a simulation at ρ = ρ comp (which is usually run at a different concentration) is

∆G 0 = ∆G comp + kT ln ρ comp ρ 0 (3.36)
For ρ 0 = 1M , a cubic simulation box with one solute molecule has a volume V ≃ 1661 Å 3 , for L ≃ 11.84 Å. This is smaller than box size usually used to solvate proteins, and lead to correction terms that are significative: 1.9 kcal/mol for a box of L = 30 Å, to 3.9 kcal/mol for a box of L = 90 Å.

The above formulation is exact with classical mechanics and dilute solutions. Another requirement of [eqn. 3.34] is that the binding involves nonionic solutes, otherwise another correction must be applied. Indeed, the chemical potential of ions can be expressed analytically in the limit of infinite dilution (Hill, [1992])

µ X = µ 0 X (p,T ) + kT ln ρ X ρ 0 X + kT ln [γ X (ρ X ,p,T )] (3.37)
At infinite dilution and small ionic strength κ the activity coefficient γ can be approximated:

kT ln [γ X (ρ X ,p,T )] ≃ κ q 2 X 2ǫ w κ 2 = 4π kT ǫ w q 2 R ρ L + q 2 L ρ L (3.38)
Even if we know an analytical expression for standard binding free energies [eqn. 3.34] which has a single standard-state dependence, more often we are interested in computing the relative binding affinity between a receptor R and two different ligands L and L ′ .

∆∆G 0 (L → L ′ ) = ∆G 0 (RL ′ ) -∆G 0 (RL) = -kT ln Q RL ′ Q RL -ln Q L ′ Q L (3.39)
We note that the standard state dependence cancels out and the calculation corresponds to measuring ratios of configuration integrals in bound and unbound states.

Free energy perturbation

From the last equations we note that free energy differences are obtained from the ratio of configuration integrals. In general, one is interested in the free energy change between two states A and B

∆G(A → B) = -kT ln Q B Q A (3.40)
For a system of hamiltonian

H = T + U (3.41)

Free energy perturbation

the canonical partition function is defined as

Z(N,V,T ) = exp - H(x,p) kT dx N dp N (3.42)
It is proportional to the configuration integral

Q(N,V,T ) = exp - U (x) kT dx N (3.43)
In biological applications we are more often interested in the isothermal-isobaric ensemble where the configuration integral is

Q(N,P,T ) = exp - U (x) + P V kT dx N dV (3.44)
The exact calculation of two partition functions Q B and Q A from computer simulations is normally too expensive. However, one can consider states A and B which are close, connected by a small perturbation in the potential function

U ′ = U + ǫ (3.45)
The corresponding free energy difference

F ′ -F = -kT ln Z ′ Z = -kT ln exp -U (x) kT exp -ǫ(x) kT dx N exp -U (x) kT dx N (3.46)
can be obtained from

F ′ -F = -kT ln exp - ǫ(x) kT U (3.47)
For a perturbation ǫ the free energy difference F ′ -F can be measured as an average over an ensemble of configurations, collected from a canonical simulation carried out with potential U(x).

The formula [eqn. 3.47] holds for small perturbations, where ensembles generated with unperturbed and perturbed potentials overlap. In general, for perturbations larger than one kcal/mol the formula cannot be used directly. To solve this problem, a good strategy is to split the path between U and U + ǫ into several smaller steps ǫ = i ǫ i (called windows) and Chapter 3. Polarizable free energy simulations then obtain the total free energy change as a sum of smaller differences. A generalization of this step-wise procedure consists in the definition of the energy function

V (λ) = (1 -λ)V A + λV B (3.48)
where A and B are two distinct thermodynamic states of the same conformational space. The parameter λ can vary continuously between 0 and 1, connecting the two states. This situation is easily implemented in MD: when λ = 0 the interactions of B are turned off, while for λ = 1 the interactions of A are turned off. For intermediate values, the system is in a hybrid thermodynamic state between A and B. This kind of transformation with unphysical states is called alchemical, as it cannot be reproduced in the real world. Using the energy function [eqn. 3.48] one can split the path between any two states A and B into many, small windows characterized by a different value of λ and then obtain the free energy difference as sum of perturbations

∆F (A → B) = F B -F A = i ∆F (λ i → λ i+1 ) (3.49)

Thermodynamic integration

The derivative of F respect to λ can be computed:

∂F (λ) ∂λ = ∂U (x,λ) ∂λ exp -U (x,λ) kT dx N exp -U (x,λ) kT dx N = ∂U (x,λ) ∂λ U (x,λ) (3.50)
The conformations used to calculate the average value of ∂U/∂λ are collected from a simulation performed at fixed λ. Calculating the derivative for several values of the coupling parameter between 0 and 1, one can integrate the results and obtain the desired free energy difference

∆F (A → B) = F B -F A = 1 0 ∂F ∂λ dλ = N -1 i=0 λ i+1 λ i ∂F ∂λ dλ (3.51)
In general, several numerical integration schemes can be used to estimate the integral. For a simple, trapezoid approach This estimation can be inaccurate. A better interpolation scheme is the cubic spline, where the function ∂G/∂λ is interpolated using a polynomial [fig. 3.4], then integrated.

∆F (λ i → λ i+1 ) = 1 2 (λ i+1 -λ i ) ∂F ∂λ (λ i+1 ) + ∂F ∂λ (λ i ) (3.52)

Free energy perturbation

Bennet acceptance ratio

The BAR method was introduced in 1976 [START_REF] Bennett | Efficient estimation of free energy differences from monte carlo data[END_REF]). It estimates the ratio of configuration integrals Q A /Q B between two states defined by potential U A and U B , following a rule inspired by Monte Carlo simulations. A detailed formulation can be found in the original paper, here we summarize the main concepts. The idea is to consider a hypothetical MC simulation where the trial move does not change the system configuration but it switches the potential function from U A to U B . For this move, equilibrium requires that

M (U B -U A ) exp(-U A ) = M (U A -U B ) exp(-U B ) (3.53)
where M is the Metropolis acceptance probability M (x) = min(1,exp(-x)) that gives Boltzmann probability. Integrating both sides of the equation and multiplying each term by 1

Q A Q A Ω M (U B -U A ) exp(-U A )dq = Q B Q B Ω M (U A -U B ) exp(-U B )dq (3.54)
one obtains the ratio

Q A Q B = < M (U A -U B ) > U B < M (U B -U A ) > U A (3.55)
As in the simple free energy perturbation method, the ratio between configuration integrals is computed using ensemble-averages. The formula can be generalized for a function f that satisfies f (x) = f (-x) ; the free energy difference between states A and B is then estimated using

∆F AB = kT ln < f (U A -U B + C) > U B < f (U B -U A + C) > U A + C (3.56)
where f (x) = 1/(1 + exp(x/kT )) is the Fermi function. C is equal to

C = ln Q A N B Q B N A (3.57)
where N A and N B are the number of configurations for the two states identified with U A and U B . The constant C is determined with an iterative procedure, until the convergence criterion

f (U A -U B + C) = f (U B -U A ) + C is satisfied. It is easy to show that if N A = N B , after convergence one has ∆F AB = C (3.58)

BAR and TI with Drude

To obtain the free energy difference between two simulation windows with BAR, one has to compute ensemble-averages. System configurations are extracted from the trajectories of separate MD simulations, that were performed using two different potentials U (λ) and U (λ + δλ). Then two ensemble-averages are computed using equation 3.56

f (U (λ + δλ) -U (λ) + C) U (λ) f (U (λ) -U (λ + δλ) + C) U (λ+δλ) (3.59)
In practice, for each MD snapshot extracted from the trajectory obtained using U (λ), one has to compute the potential energy with U (λ) and also with U (λ + δλ). However, the Drude particle displacements of these snapshots {d ⋆ i (λ)} satisfy the SCF condition ∂U /∂d=0 for the potential U (λ) but not for the potential U (λ + δλ). For this reason, the application of BAR with Drude is not formally correct. Before to compute U (λ + δλ) for a shapshot of the U (λ) trajectory, one should minimize the Drude particles using U (λ + δλ). However, this solution is not practicable since it is too expensive.

The problem does not persist if free energy differences are obtained with TI, where one has to compute derivatives:

∂F (λ) ∂λ = ∂U (λ) ∂λ U (λ)
(3.60)

Free energy perturbation

In practice, the derivative ∂U/∂λ can be obtained by a finite-difference method, by taking the difference between energies computed at two slightly different λ values:

∂U ∂λ = U (λ + δλ) -U (λ) δλ (3.61)
Thus, with TI, we compute energies at, or very close to the λ value used to produce the trajectory. In this case, the relaxed displacements are supposed to be similar d ⋆ i (λ+δλ) ≃ d ⋆ i (λ) and the SCF condition is still well approximated. Indeed, evaluating the Taylor expansion of equation 3.26 for small δd i the first derivative is equal to zero, leaving a quadratic term that is negligible when

δd i = d ⋆ i (λ + δλ) -d ⋆ i (λ) ≃ O(δλ 2 ) (3.62)
Few comments can be made. First, we note that while BAR is not formally correct with Drude, it can be applied to the AMOEBA polarizable model, where the induced dipoles are automatically adjusted whenever the parameter λ is changed. Second, Drude simulations are run in the "Cold oscillators regime" where the oscillators are thermalized at a low temperature T D in order to approximate the SCF condition ∂U/∂ d = 0. When combined with a barostat, MD samples a (N,P,T,T D ) ensemble. Several studies describe methods to compute free energies for similar dual-temperature systems, where some degrees of freedom evolve adiabatically respect to others equilibrated at room temperature (Rosso et al. ). Even if their application for Drude is not straightforward, the problem could be studied mode in detail in the future.

Alchemical transformation

The calculation of binding free energies for the reaction [eqn. 3.27] has been extensively discussed in literature. The article from Boresch et al. [2003] includes an overview about techniques used to obtain absolute binding free energies, with a discussion about the standard state dependence. In figure 3.5 we show the thermodynamic cycle corresponding to the Double Decoupling method. It relies on the decomposition of the binding reaction

R s + L s ∆G --→ bind (RL) s (3.63) into two transformations L s ∆G --→ 1 L g (RL) s ∆G --→ 2 R s : L g ∆G --→ restr R s + L g (3.64)
the subscript s is used to indicate the solvated system. In the first reaction, the ligand is decoupled from the solvent and transferred to the gas phase (indicated using the subscript g).

The free energy cost for this transformation is called ∆G 1 . The second reaction is composed of two parts. First, the ligand is transferred to the gas phase under specitic restraints introduced in order to keep L bound to R, even when their interactions are turned off. The free energy cost for this reaction is called ∆G 2 . The endpoint is called R s :L g to indicate that R and L are bound just because of the introduced restraints. Then, the restraint potential is removed and the ligand is fully transferred to the gas phase. The free energy cost for this reaction is called ∆G restr . The the energy to apply the restraint in the bound state (RL) s is usually neglected. Finally, To compute the free energy necessary to remove the restraint between receptor and ligand, one can use a simulation that starts from the endpoint of reaction ( 2), where the restraint potential is gradually turned off. In practice, this is always complicated [START_REF] Boresch | Absolute binding free energies: a quantitative approach for their calculation[END_REF]). The most straightforward method to obtain this free energy difference is to compute it ana-6. Artefacts in charging free energy calculations lytically: for a cubic box, ∆G restr = -kT ln(Q restr /L 3 ) where V = L 3 is the volume of the simulation box. For example, for the simple case of an ion linked to the receptor with an harmonic spring one can estimate the configuration integral exploiting the spherical symmetry

∆G bind = ∆G 2 + ∆G restr -∆G 1 (3.
Q restr =

Artefacts in charging free energy calculations

When comparing simulation results with theory and experiments, several artefacts must be taken into account. In most cases, electrostatic interactions are computed on a periodic, infinite lattice using Particle Mesh Ewald (PME). As a result, alchemical free energies differ respect to the real quantities of interest depending on box size, charge distributions, the presence of counterions, and the nature of the solvent (Rocklin et al. [2013]). In general, one has to compute corrections such that the real free energy change is obtained from the alchemical free energy difference plus a term which can be computed analitically or numerically, depending on the situations. In few situations simulations need to be postprocessed, for example using continuum electrostatics to estimate several artificial contributions. Otherwise, a well-established theoretical background allows to estimate some of these corrections analytically.

PME with tinfoil boundary conditions: neutralizing gellium

We consider a simulation box and its infinite periodic images, where the total net charge of each box Q = 0. To avoid the (unphysical) situation where the infinite, periodic system has a divergent potential, the simulation box must be neutral. A possible strategy consists in introducing one (or more) counterions to balance the total net charge Q. For large boxes, the counterions will be sufficiently isolated, without affecting the simulation interacting with the system or perturbing each other's solvation. However, a more straightforward and practical solution is to introduce a uniform compensating charge density. For each charge q inserted in the system, a neutralizing charge -q is spread uniformly over the entire box volume. This charge density is called the gellium. Its charge density depends on the volume of the box: ρ = -q/V = -q/L 3 . It generates a uniform potential that now we call φ gel . Since ∇φ gel (x) = 0 for all points inside the box, the gellium does not introduce an electric field in the MD simulation or contribute to the forces. However, the role of the gellium is to shift the electrostatic potential φ in such a way that its average over the entire simulation box is zero:

φ(x) box = 0 (3.67)
To see this, we consider a charge q inserted in the system. It will give rise to infinite periodic images and an associated gellium, so that at a point r inside the box, the associated charge density is

ρ(r) = q ∞ i=0 (δ(r -r i ) -L -3 ) (3.68)
where r i is the position of the perturbing charge or one of its images. The energy change per box is

∆U = box drρ(r)φ(r) = box drq δ(r -r i ) -L -3 φ(r) = q(φ(r i ) -φ box ) (3.69)
The perturbing charge q thus experiences a shifted potential φ ′ = φ(r i )φ box . In effect, the gellium has shifted the mean box potential φ to zero; in other words, the system is electrically grounded.

The potential shift φ box depends on the nature of the solute and the volume fraction of the solvent in the box [START_REF] Lin | An overview of electrostatic free energy computations for solutions and proteins[END_REF]). When free energy calculations are performed exploiting thermodynamic cycles, different simulations can be run in different simulation boxes. Since different simulation boxes can be shifted by a different potential φ box , this effect must be taken into account. One example is the calculation of relative binding affinities illustrated 6. Artefacts in charging free energy calculations in figure 3.6, where different shifts can artificially overstabilize bound or unbound states. A practical example to illustrate how we can correct for this effect will be given in the next chapter.

Solvent polarization artefacts in a periodic system

Other artefacts are introduced by the the fact that electrostatic interactions are calculated in a periodic system. We consider the simple case of a spherical ion placed at the center of a periodic, cubic box of size L. Because of periodicity, the ion will interact with its periodic images and the neutralizing gellium. These interactions arising from the lattice summation are also called ion self-energy and are related to the box volume through a 1/L dependency. The self energy must be considered when, computing alchemical charging free energies ∆G(L), one wants to extrapolate the so called large box limit ∆G(L → ∞) where artefacts are removed. As pointed out by Hummer et al. [1996], the self contribution to the charging free energy of a small spherical ion of charge q is

∆G self = q 2 ξ 2L (3.70)
where ξ = -2.837297 is the Wigner constant for a cubic lattice. The term includes interactions between the ion, its images and the gellium. However, ∆G self is not sufficient to remove all the artefacts. As observed by Figueirido et al. [1997], the total charging free energy has a weaker dependence on the simulation box volume because of the solvent screening. For high dielectric solvents, the real deviation between ∆G(L) and ∆G(L → ∞) is smaller. It depends on the solvent dielectric constant and is proportional to q 2 /(ǫ W L). Another finite-size artefact can be removed considering another correction term which scales as L -3 (Hummer et al. [1997]), valid for spherical ions of radius R. The final correction formula that removes all these artefacts is

∆G(L → ∞) = ∆G(L) + q 2 |ξ| 2ǫ W L + 2πq 2 R 2 3L 3 ǫ w -1 ǫ w + O(R 2 L -3 ǫ -3 w ) (3.71)
To obtain charging free energies in kcal/mol, the Wigner constant must be multiplied by a factor to convert e 2 /Å (in CHARMM this is the CCELEC constant = 332.0716). More information can be found in [START_REF] Simonson | Concepts and protocols for electrostatic free energies[END_REF]. Here we report some important observations. First, the charging free energy extracted from a simulation ∆G(L) is smaller than the one obtained in the large box limit. This is intuitive considering that the positive ion will interact with the neutralizing gellium, which is a uniform negative charge density spread over the simulation box volume. Second, this analytical formula is valid for spherical ions in cubic boxes. For more complex solutes and different simulation boxes, the correction must be computed using continuum theory. Numerical schemes can be found in [START_REF] Rocklin | Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: An accurate correction scheme for electrostatic finite-size effects[END_REF]. The Chapter 4

PDZ-peptide binding specificity with polarizable free energy simulations 1 Introduction

Despite the abundant experimental data, our understanding and ability to rationally engineer PDZ/peptide interactions is limited. For the best-studied domains, experiments have revealed the binding thermodynamics of 20-40 protein or peptide variants, but only a few x-ray structures are available. Even with x-ray structures, many details are unclear. Protein:peptide binding affinity and specificity depend on many factors, difficult to quantify using experiments. These are short-range interactions in the binding pocket, longer-range electrostatic interactions, dielectric shielding by protein and solvent, ordered waters in the binding site, the structure and flexibility of the unbound peptide, and the conformational dynamics of the receptor.

Computer simulations are a complementary tool and several approaches have been proposed to predict the binding energetics for protein/peptide complexes. Some of them are semi-empirical and use a Linear Interaction Energy or "PB/LIE" free energy function based on molecular mechanics plus an implicit solvent. This function has several contributions (usually van der Waals, continuum electrostatics and a solvent accessible area term) weighted by different factors, which are optimized to reproduce binding affinities for a training dataset.

Of particular interest are more expensive, nonempirical methods that predict binding free energy differences from molecular dynamics (MD) simulations without any adjustable parameters. These methods are formally exact and are not optimized for a particular test system. Nonempirical methods are based on thermodynamic perturbation theory, and extract free Chapter 4. PDZ-peptide binding specificity with polarizable free energy simulations energies from series of simulations where a mutating side chain of interest is alchemically mutated from one type to another. This class of methods is often referred to as free energy perturbation or FEP [START_REF] Zwanzig | High-temperature equation of state by a perturbation method. i. nonpolar gases[END_REF] ; Straatsma & McCammon [1992] ; [START_REF] Kollman | Free energy calculations: Applications to chemical and biochemical phenomena[END_REF]).

There are mainly two sources of errors that can influence the prediction of binding free energies with the FEP. The first is related to sampling. FEP is expensive and usually needs long simulations which may last several hours or days, even when performed on a supercomputer. Sampling requirements may be especially high for ionic mutations, since they involve long range interactions and long timescales that characterize dielectric relaxation in proteins [START_REF] Simonson | Electrostatics and dynamics of proteins[END_REF]). A second source of errors is the molecular mechanics force field used to run MD simulations and to extract the free energy differences. As anticipated in the last chapter, a strong limitation is the simple treatment of electronic polarization modeled implicitly through a fixed set of atomic partial charge, as in most of the widely-used, "additive" force fields. The use of a polarizable force field could be especially important when comparing ligands that differ by one or more ionic side chain mutations. However, when using more sophisticated models the computational time needed to accomplish sufficient conformational sampling is increased.

Experimental data is known for 50 Tiam1:peptide complexes, 25 of them are peptide variants [START_REF] Shepherd | The tiam1 PDZ domain couples to syndecan1 and promotes cell-matrix adhesion[END_REF][START_REF] Shepherd | Distinct ligand specificity of the tiam1 and tiam2 PDZ domains[END_REF] ; Liu et al. [2013]). Binding free energies are within a 2.2 kcal/mol range, with mean experimental uncertainty of 0.2 kcal/mol. Some peptides have high or non measurable dissociation constant. Previous FEP calculations were run using Amber ff99SB and produced good results in the prediction of relative binding affinities for six nonionic mutations. Results are listed in our recently published paper [START_REF] Panel | Accurate PDZ/peptide binding specificity with additive and polarizable free energy simulations[END_REF]). Errors were low, with mean unsigned error (MUE) of 0.32 kcal/mol, root mean square deviation (RMSD) of 0.37 kcal/mol and correlation coefficient of 0.84. For other six ionic mutations, mean errors with the additive force field were larger, and there was one very large error of 3 kcal/mol [tab. 4.1]. Among these six mutants, four of them displayed errors of more than 1 kcal/mol: three are mutations of the peptide position P 4 E4K, E4L, K4L and one mutation was in the protein K912E. We use italics to indicate mutations, for example E4K will refer to the alchemical Sdc1→ Sdc1-E4K mutation.

In the present chapter, we report alchemical free energy simulations of the Tiam1 PDZ domain using two additive force fields (Amber ff99SB and Charmm C36) and the Drude polarizable force field. Our main goal was to determine the accuracy of free energy perturbation (FEP) in characterizing the binding energetics of PDZ/peptide complexes. We were also interested in determining if the explicit threatment of electronic polarization can improve results for ionic mutations, for which additive force fields are known to be imprecise. In this work, we studied four ionic mutations using the Charmm additive force field C36 and the Drude pola-

Methods

rizable force field, in order to investigate if the explicit threatment of electronic polarizability can improve predictions. All calculations (with polarizable or additive force field) relied on the Dual Topology Approach. To our knowledge, FEP has never been used to study protein:ligand binding free energy changes with the Drude polarizable force field.

Simulations were done using periodic boundary conditions and particle mesh Ewald (PME) summation for long-range electrostatic interactions [START_REF] Darden | Computational biochemistry and biophysics, treatment of long-range forces and potential[END_REF]). Special care was taken to extrapolate the computed free energies correctly to the thermodynamic limit of an infinitely large simulation box. Applying the Drude polarizable force field to the four ionic mutations with the largest errors, two of the errors were significantly reduced and the largest one decreased from 3.0 to 1.3 kcal/mol. Thus, the polarizable force field led to a good accuracy for the ionic mutations, although not as good as that previously obtained for the nonionic mutations with ff99SB. Overall, we have established that FEP can be used successfully to understand and potentially engineer PDZ/peptide binding specificity, and that a polarizable force field is necessary to obtain good accuracy for some of the ionic mutations. The FEP methodology can thus be used in a predictive way, in synergy with experiment and other computational approaches, to design PDZ/peptide specificity.

Methods

Structural models and simulation setup

Complexes involving mutants of the Tiam1 PDZ domain or the Sdc1 peptide were built from the wild-type/Sdc1 x-ray structure (pdb access code 4GVD) using the SCWRL4 program [START_REF] Krivov | Improved prediction of protein sidechain conformations with scwrl4[END_REF]). For each complex, the mutated side chain and nearby side chains were positioned by SCWRL4, with the rest of the structure held fixed. There were six missing residues in the wild-type/Sdc1 x-ray structure, in the exposed β1-β2 and β2-β3 loops, which were rebuilt using the program MODELLER (Fiser et al. [2000]). The unbound peptides were all modeled by extracting the bound conformation from the corresponding complex.

Simulation setup was carried out using the CHARMM program. The system was first built using the C36 force field. Each complex or unbound peptide was immersed in a cubic box of size L =80 Å of pre-equilibrated TIP3P water, deleting solvent molecules within the volume occupied by the solute, using a distance threshold of 3.0 Å from the complex heavy atoms. Then the box was trimmed, from cubic to truncated octahedron, deleting solvent molecules at the cube edges. The wild type/Sdc1 system was neutralized using sodium ions. Solvated systems were equilibrated using CHARMM and the additive force field C36, in a step-wise procedure.

Chapter 4. PDZ-peptide binding specificity with polarizable free energy simulations 500 ps of dynamics were run increasing timestep and system temperature, decreasing initial restraints applied to the solute heavy atoms. Then 40 ns of MD were run. After this first, "additive" equilibration, Drude particles and Lone pairs were added using a CHARMM script based on the drude-prepper tool of the CHARMM-GUI (Jo et al. [2008[START_REF] Jo | CHARMM-GUI 10 years for biomolecular modeling and simulation[END_REF]). The script reads coordinates of the additive system, generates the Drude structure file (PSF) and centers each Drude particle at the position of its reference atom. At this point, it was necessary to relax Drude particles, since the generated Drude positions were unphysical. A first energy minimization was run (200 steps) to relax the system with restrained centers and unrestrained Drude particles. A second minimization (100 steps) was run relaxing the whole system without restraints. Starting from the minimized conformation, 20 ns of NPT equilibration followed by 100 ns of production dynamics were run. Minimization and MD were run using NAMD 2.12. Molecular dynamics simulations with the C36 and Drude force fields were done at room temperature and pressure, using Langevin dynamics with a Langevin Piston Nosé-Hoover barostat (Feller et al. [1995]). Long-range electrostatic interactions were treated with a PME approach with tin-foil boundary conditions [START_REF] Darden | Computational biochemistry and biophysics, treatment of long-range forces and potential[END_REF]). With Drude, the MD timestep was 1 fs, instead of 2 fs used for the simulations with C36. The Drude particles had the default mass of 0.4 atomic units, and their temperature was controlled with a Langevin thermostat with T D = 1 K.

Alchemical MD simulations

The FEP approach [START_REF] Simonson | Free energy simulations come of age: Proteinligand recognition[END_REF] ; Tuckerman [2007] ; [START_REF] Jorgensen | Perspective on free energy perturbation calculations for chemical equilibria[END_REF]) was used to calculate the binding free energy differences between variants A and B using equation 3.39. Two separate alchemical simulations were performed. In the first simulation, the mutating side chain was alchemically transformed in the bound state, represented by the protein:peptide complex in water solution. In the second simulation, the transformation was performed in the unbound state, represented by the unbound peptide (or unbound protein) in solution. From these two simulations, we extracted free energy differences ∆G bound (A → B) and ∆G unbound (A → B). The relative binding free energy for a generic mutation A → B was then obtained by difference

∆∆G bind (A → B) = ∆G bound (A → B) -∆G unbound (A → B)
Let us assume we are comparing the wild type Sdc1 peptide to one of its point mutants. For this pair, we consider a hybrid peptide, which has two side chains at the mutated position, one of each type. The energy function U depends on a coupling coordinate λ that scales selected electrostatic and van der Waals energy terms. When λ = 0 (respectively, 1), the 2. Methods mutant (respectively, wild type) side chain is decoupled from its surroundings, retaining only its covalent interactions with its own backbone. For intermediate λ values, both side chains are present, with intermediate weights. U has the form

U (λ) = U LJ (λ) + U elec (λ) (4.1)
Considering two atoms i and j at distance r ij , their Lennard Jones interactions are modified using a coupling parameters λ LJ

U LJ (r ij ; λ LJ ) = λ LJ ǫ ij      R min,2 ij r 2 ij + δ(1 -λ LJ )   6 -   R min,2 ij r 2 ij + δ(1 -λ LJ )   3    (4.2)
The electrostatic potential was modified using another coupling paramer λ elec . For simulations with C36:

U elec (λ) = λ elec q i q j r ij (4.3)
while for the Drude polarizable force field were scaled charge-charge, charge-dipole and dipoledipole interactions:

U elec (λ) = λ elec q A,i • q A,j |r A,i -rA,j | + q D,i • q D,j |r D,i -rD,j | + q D,i • q A,j |r D,i -rA,j | + q A,i • q D,j |r A,i -rD,j | (4.4)
The van der Waals interactions are thus modified using a soft-core functional form, ensuring a gradual transformation. The shifting parameter δ was set to its default value in NAMD δ = 5. Proceeding in this way, clashes at endpoints are avoided when λ is close to 0 or 1 [START_REF] Beutler | Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations[END_REF] ; [START_REF] Zacharias | Separation-shifted scaling, a new scaling method for lennard-jones interactions in thermodynamic integration[END_REF]).

Calculations were done using the alch facility of NAMD. When performing simulations with the Drude force field, particular attention is also needed to balance appearing and disappearing electrostatic interactions. This is related to the stability of MD simulations, since one has to avoid the polarization catasfrophe occurring when a Drude oscillator deviates too much from its reference atom. In NAMD, simulations often become unstable when one atom-Drude distance exceeds the default threshold of 0.2 Å. We observed that when the LJ repulsion between a group of atoms is annihilated too early compared to electrostatics, positive charges belonging to nearby groups can get too close to the negative Drude charges, attracting them away. As a consequence, these oscillators exceed the threshold distance and NAMD prints an error message interrupting the MD simulation. This situation shows up often when water We set the NAMD parameter alchElecLambdaStart=0.4 or 0.5, depending whether the simulations were run using the additive or the Drude force field. This means the van der Waals terms of the new side chain are introduced starting at λ = 0, but the electrostatic terms are introduced starting at λ = 0.4 or λ = 0.5. With a value of 0.4, electrostatic terms of the disappearing side chain are completely decoupled at λ = 0.6, while its van der Waals terms disappear at λ = 1. In figure 4.1, we show the protocol adopted for Drude simulations.

Methods

Figure 4.2 -Water-acceptor interaction with partially-annihilated LJ interactions. Water hydrogens approach the hydrogen bond acceptors, increasing displacements of their Drude particle.

Drude dual topology

The structure files (PSF) for alchemical MD were built using special topology files called dual topologies applied on the mutating residue. To simulate a mutation between two different side chains A and B, A and B were both attached to the same C α atom. Bonded and nonbonded interactions between the two side chains were excluded. In figure 4.3 we give an example of dual residue called HEK, composed of glutamate and lysine side chains attached to the same backbone. In the first part of the file, we list atoms of the HEK backbone, shared by the two side chains. Backbone atoms form a neutral, zero-charge group. In the second part of the file we list the side chains, attached to the same C α through two bonds CBE-CA and CBK-CA. At each ATOM line we indicate the atom name, type, charge, polarizability and Thole factor, as usual in the Drude force field. Then a list of atom names belonging to the other side chain is given. These define nonbonded interactions between GLU and LYS that will be excluded. These exclusions are applied for all values of λ. In the given example, we used "[...]" to abbreviate the lines present in the real topology file, where all the atoms of the other side chain are listed.

In HEK, the C β atoms (CBE and CBK) are depolarized (they do not carry a Drude particle) even if according to the foorce field they should be polarizable. This modification was necessary to exclude spurious dipole-dipole interactions between these two atoms. 2. Methods free energy to depolarize them is computed separetly, as explained below. Other undesired, spurious interactions between the two side chains must be removed when atom types of CA, CBE and CBK form an angle for which the force field expects an harmonic term. The same holds when bridging atoms are involved in a dihedral term of the energy function. This does not occour too often in C36 or Drude and is more frequent in Amber, since Amber has fewer atom types. In this case, undesired terms must be manually deleted within CHARMM before producing the system PSF. The so defined non-bonded lists are exploited to compute non-bonded interactions as follows (Harder et al. [2006]). Interactions between core charges q A are excluded between 1-2 and 1-3 pairs. Interactions between Drude particles q D and core charges q A are excluded between 1-2 and 1-3 pairs. Interactions between all core charges q A and all Drude particles q D are computed for 1-4 pairs and beyond. Interactions between Drude charges q D of 1-2 and 1-3 pairs are computed and screened according to their Thole parameters, using equation 3.21. With NAMD, interactions between Drude particles are also screened for nonbonded 1-4 pairs and beyond, if their distance is lower than a threshold value set by default to 5.0 Å. However, the screening between 1-4 pairs is disabled when running the FEP using the alch keyword.

Spurious interactions with Drude and dual topology

Chapter 4. PDZ-peptide binding specificity with polarizable free energy simulations Summing up, as mentioned in the official CHARMM documentation "If two atoms do not see each other in the non-polarizable force field, their dipoles (and only their dipoles, not their net partial charges) will see each other in the polarizable force field". The particular bond hierarchy is used to compute the screened dipole-dipole interactions between 1-2 and 1-3 pairs, as proposed by Thole [1981]. However, this introduces a spurious interaction in the dual topologies. In the HEK example, C β atoms CBE and CBK are in a 1-3 relationship as they are both bonded to the backbone atom C α . As a consequence, their Drude particles DCBE and DCBK have a 1-3 relationship: their interaction is computed and screened according to their Thole parameters, using equation 3.21. This is a spurious interaction, since in the FEP the two side chains of the mutating residue should not interact. It occours because the Drude force field is only partially supported by the alch facility of Namd, which does not correctly modify the nonbonded-lists in order to exclude the dipole-dipole interaction between the two

C β atoms.
As a consequence, the current Namd implementation of alch does not allow to run FEP simulations with Drude using dual topologies where the bridging atom is linked to two atoms carrying Drude particles. However, the problem can be solved following an ad-hoc protocol. Our dual topologies were built removing the Drude particle attached to both CB atoms of the two side chains. These atoms were depolarized. In practice, we condensed the total charge of the CB+DCB couple to the CB atom Q = q A and we have set the CB polarizability and Thole factor to zero [fig. 4.5]. Alchemical simulations were then run using these modified topology files, where the spurious interaction were removed. Removing C β polarizability has a free energy cost that must be considered. We introduced additional transformations at the beginning and the end of the mutation process, where the

Methods

C β Drude charges were gradually removed (for A) or introduced (for B). These simulations correct the alchemical free energy difference for the transformation with dual topology λ = 0 → 1, as shown in figure 4.6. We defined coupling parameters µ(A) and µ(B) that scaled the polarizabilities on the two C β atoms. When µ(A)=1, the C β of A had its full polarizability; when µ(A)=0, its polarizability had been removed. To change the polarizability α of an atom, we adjusted the atom+Drude charges q A and q D using the relationship

α = q 2 D k D q A + q D = Q (4.5)
We varied µ(A) from 1 to 0 in four MD windows. During these windows, side chain A was fully coupled (the coupling coordinate λ was zero), whereas side-chain B was decoupled and the polarizability of its C β was turned off. We then did the alchemical A→B transformation with both C β polarizabilities turned off. Finally, we turned on the C β polarizability of B,µ=0→1.

For the polarizability transformations, we used MD windows where µ(A) (respectively, µ(B)) was set to 0, 0.3, 0.7 and 1 (20 ns each). The value of the polarizability was controlled by editing charges in the NAMD protein structure file. Free energies were computed with TI and added to the λ=0→1 results.

PME correction

Using PME under periodic boundary conditions, the electrostatic potential of the infinite, periodic system is infinite unless the simulation boxes are grounded. In alchemical simulations

Chapter 4. PDZ-peptide binding specificity with polarizable free energy simulations of ionic mutations, the transformation of the mutating side chain along the path λ = 0 → 1 is a charging process and the endpoints cannot both be neutral. There are two possible solutions, previously discussed in the last chapter. The first is to introduce one or more counterions in the simulation box, which are alchemically transformed during the FEP in such a way to neutralize the system for each value of λ. The free energy can be easily corrected keeping the counterion far from the solute, since it will contribute independently to the charging free energy. A more common strategy is to introduce the gellium, that has an important effect: it leads to an electrostatic potential, averaged over the box volume, that is zero at all times.

In other words, it grounds all simulation boxes to zero. However, the periodic system with the gellium is not a faithful representation of the true physical endpoints. These correspond to a more realistic situation, where the solute (the protein:peptide complex or the unbound peptide) is solvated in a macroscopic but finite solvent volume, such the large spherical cluster of water molecules, like shown in figure 4.7. In the macroscopic clusters, the average electrostatic potential per box φ cluster is not zero, it is close to the vapor/liquid boundary potential Φ lv . For TIP3P water, Φ lv = -520 mV [START_REF] Lin | An overview of electrostatic free energy computations for solutions and proteins[END_REF] ; Harder & Roux [2008]) while for SWM4-NDP its value is -545 mV [START_REF] Lamoureux | A polarizable model of water for molecular dynamics simulations of biomolecules[END_REF]). The grounded PME potential of a simulation box φ box is so shifted upwards, by a quantity close to this boundary potential. The exact value of this shift depends on the fraction of box occupied by the water and on the nature of the solute. Calculations [START_REF] Lin | An overview of electrostatic free energy computations for solutions and proteins[END_REF]) showed that for a small spherical ion solvated in a cubic box, the shift scales as L -3 where L is the size of the box. For more heterogeneous and different solutes, the shift will be 2. Methods different even if the same box size is used.

The potential shift for the protein/peptide complex and of the unbound peptide simulations will be different, since these simulations are performed in different boxes, and the nature of the solute is different. As a consequence, inserting a charge on the peptide in the simulation box of the complex (bound state), this charge will be subjected to a shifted potential which is different respect to the one in the simulation box of the unbound peptide. A mutation that inserts a net charge onto the peptide in both bound and unbound states requires a free energy correction, which is computed following the protocol described in Lin et al. [2014] ; [START_REF] Simonson | Concepts and protocols for electrostatic free energies[END_REF].

For a given solute and mutation endpoint we compute, from the MD trajectory, the mean potential in the solvent region of the box, defined to be > 12 Å away from the solute or its periodic images. For a small solute like a peptide, it is close to zero, the overall box average with PME. For a protein solute, the overall box average is an average over a solvent region (far from the solute) and a protein region (the rest of the box):

V p Φ p + V w Φ w = 0 (4.6)
Here, V w and V p are the volumes of the solvent and protein regions and the brackets indicate a spatial average of the potential φ over either region. Thus,

Φ w = V p V w Φ p (4.7)
The potential in the solvent region of a protein solute is offset, instead of matching that in the peptide simulation (zero). With both TIP3P and Drude SWM4 water, the offset is positive, making it artificially easier to insert a positive charge into the protein than the peptide.

For a given solute, we denote the solvent potential φ w (solute). If the mutation introduces a net charge q on the peptide, this charge will then experience a potential that is shifted by φ w (complex) in the bound state and φ w (peptide) in the unbound state. Multiplying q by the difference gives the free energy correction. In practice, we average the free energy correction over the simulations of the initial and final peptides (before and after the mutation): peptide, complex, peptide ′ and complex ′ , respectively. The potential shift correction is thus:

∆∆G Φ P M E = q 2 Φ w (complex) + Φ w (complex ′ ) -Φ w (peptide) + Φ w (peptide ′ ) (4.8)
Chapter 4. PDZ-peptide binding specificity with polarizable free energy simulations

In effect, we average over a calculation where q is inserted into the initial peptide and one where q is removed from the final peptide. In practice, the two values are very similar. We verified that the potential is almost uniform throughout the solvent region, and very similar for complex and complex ′ , supporting the idea that φ w (solute) accurately represents the PME potential shift.

The potential throughout the simulation box was obtained for each frame of each endpoint trajectory by the VMD plugin PMEpot (Humphrey et al. [1996]), which prints out the potential values on a cubic grid spanning the box. Time and spatial averaging were done with a shell script. MD simulations with a cubic box (of the same volume) were specifically done for each system to allow use of the plugin (which does not support a truncated icosahedral box). With Drude, for a +1 charge insertion an artificial contribution of ∆∆G φ = -0.7 kcal/mol was obtained and must be corrected for. For the same process C36 has a smaller shift of -0.39 kcal/mol, close to the one obtained earlier for ff99SB which equals -0.45 kcal/mol.

Results

We computed binding free energies differences ∆∆G for four mutations using the additive force field C36 and the Drude polarizable force field. The same variants were previously studied using the Amber ff99SB additive force field [START_REF] Panel | Accurate PDZ/peptide binding specificity with additive and polarizable free energy simulations[END_REF]). Results are shown in table 4.1, including the ∆∆Gs of two other ionic mutations previously obtained using ff99SB. Peptide positions are numbered backward from the C-terminus (position 0 or P0). In the table we list all contributions to the binding free energy differences. Care was taken here to extrapolate the MD results to the thermodynamic limit of a very large simulation box. The Drude PME correction Φ P M E was 0.70 kcal/mol for E4L, -0.70 for K4L and 1.40 kcal/mol for E4K. For K912E/Sdc1, the correction was zero, since the V p /V w ratio in bound and unbound states was essentially the same. PME corrections for additive force fields were smaller, because the mean electrostatic potential in TIP3P water is smaller than that in Drude SWM4-NDP water. With Drude, there were additional transformations µ A , µ B at the beginning and end of the mutation that removed/restored the C β polarizabilities, computed with TI. These contributions were between -0.3 and 0.1 kcal/mol. The alchemical steps computed using BAR λ = 0 → 1 contributed between -1.0 and 2.1 kcal/mol with Drude, while for C36 and ff99SB they were always positive and between 0.8 and 3.2 kcal/mol. Some of the Drude contributions were analyzed also with TI, obtaining the same result. Indeed, even if BAR is not formally correct when used with Drude (as described before in the text), the moderate errors for the bound and unbound states are mostly canceled when computing the binding free energy differences ∆∆G = ∆G bound -

Results

∆G unf ound . We note that all three contributions to the Drude binding free energy change λ = 0 → 1, Φ P M E and µ A , µ B were of roughly the same order of magnitude.

Mutation Contribution Exp. Drude ff99SB C36 E4K µ A ,µ B 0.1 - - λ = 0 → 1 -1.0 0.8 1.0 Φ P M E
1.4 0.9 0.8 ∆∆G bind 0.8 0.5 1.9 1.8 err.

-

-0.3 1.1 1.0 E4L µ A ,µ B -0.2 - - λ = 0 → 1 1.4(1.3) 3.1 3.2 Φ P M E 0.7 0.5 0.4 ∆∆G bind 0.6 1.9(1.8) 3.6 3.5 err. - 1.3 3.0 2.9 K4L µ A ,µ B -0.3 - - λ = 0 → 1 2.1(2.1) 1.4 2.2 Φ P M E -0.7 -0.5 -0.4 ∆∆G bind -0.2 1.1(1.1) 0.9 1.8 err. - 1.3 1.1 2.0 K912E µ A ,µ B 0.7 - - λ = 0 → 1 2.1 2.0 1.6 Φ P M E 0.0 0.0 0.0 ∆∆G bind 1.0 2.
= 0 → 1 - 1.5 - Φ P M E - 0.9 - ∆∆G bind 1.3 - 2.3 - err. - - 1.0 - K912E/Caspr4 λ = 0 → 1 - 0.5 - Φ P M E - 0.0 - ∆∆G bind 0.7 - 0.5 - err. - - -0.2 - rmsd - - 1.50 - mue - - 1.
23 -Table 4.1 -Relative binding free energies for ionic mutations obtained with polarizable and additive force fields. Specific contributions are shown on separate lines. Alchemical free energy differences λ = 0 → 1 were obtained using BAR. Two values were also computed with TI, values are shown in parenthesis. The table is separated into two parts: Rmsd and Mue are first computed considering the first 4 mutations for which Drude and C36 results are known (E4K, E4L, K4L and K912E). Then, rmsd and mue are calculated for the six ionic mutations previously studied using ff99SB (Panel et al. [2018]).

The overall Drude binding free energy changes had errors of -0.3, 1.3, 1.3, and 1.6 kcal/mol, respectively, for an RMS error of 1.2 kcal/mol. The same variants with ff99SB and C36 gave RMS error of 1.8 and 1.9 kcal/mol. Thus, the RMS error was reduced with Drude. The improvement was mainly due to E4L, which had the largest error with both the additive force fields. The C36 force field gave ∆∆G values for E4K and E4L within 0.1 kcal/mol of AMBER ff99SB and a similar RMS error. The closure errors for the thermodynamic cycle (Sdc1/E4K/E4L/Sdc1 ) with Drude and C36 were very small, 0.3 and 0.0 kcal/mol, respectively, suggesting convergence was good. The ff99SB closure error was slightly larger, 0.8 kcal/mol. For the E4L and E4K ionic mutations, previous very long simulations were run with ff99SB (495 and 770 ns, respectively) to establish that the errors were due to force-field limitations, not sampling. The agreement between C36 and ff99SB supports the idea that the large E4K errors were due to the lack of explicit polarizability.

For the other five ionic mutations, agreement with experiment was still rather good using AMBER ff99SB, with an RMSD of 0.90 kcal/mol. Errors of -1 kcal/mol were also seen recently with this force field for ionic mutations in the tyrosyl-tRNA synthetase enzyme binding its amino acid substrate (Simonson et al. [2016]).

The main interactions that stabilize the different complexes were identified by structure analysis. With the AMBER additive force field, the Sdc1-E4 side chain makes a salt bridge with either Arg871 or Lys912 half of the time [fig. 4.9], whereas with the Drude polarizable force field, a Glu4-Lys921 salt bridge is present only 1/4 of the time as shown in figure 4. [START_REF] Goldstein | Evolutionary perspectives on protein thermodynamics[END_REF]. No salt bridge is present in the x-ray complex. Thus, the additive force field overpopulates the salt-bridge conformations, leading to an overstabilization of wild-type/Sdc1 relative to the E4K and E4L mutants. In the mutant E4K and E4L complexes, there are no salt bridges involving position P4, either with AMBER ff99SB or with Drude. K4 mostly makes polar interactions with solvent and L4 is hydrophobic. Although the Drude model is less affected by salt-bridge overstabilization, it still predicts an E4-Lys912 interaction 1/4 of the time, and this appears to account for the positive ∆∆G errors for the E4L and K912E mutations obtained with Drude.

WT:Sdc1

WT:E4K WT:E4L

Conclusions

To our knowledge, this study was the first application of the Drude force field to compute relative binding affinities for a set of protein:ligand complexes. The described protocol can be applied to any alchemical FEP when a dual topology is preferable. Technical issues were highlighted, mainly related to the actual NAMD implementation of FEP. Simulations were needed to compute the additional contributions µ A , µ B to restore the deleted CB polarizabilities. This increases the computational cost of the Drude FEP, which is already bigger (a factor 4) compared to Amber ff99SB or C36. However, future versions of the alch facility of Namd could solve this problem, deleting the spurious interactions. Overall Drude results are

Conclusions

encouraging: the E4L error was reduced from 3.0 to 1.3 kcal/mol, the E4K error was reduced from 1.1 to 0.3 kcal/mol; the K4L error was reduced by 0.7 kcal/mol compared to C36, and was similar to the ff99SB error. Only the K912E/Sdc1 error increased somewhat with Drude, from 1.0 kcal/mol with ff99SB to 1.6 kcal/mol. K912 is highly exposed to solvent and does not form a salt bridge. Established additive force fields are well parameterized for this situation, whereas we can speculate that the more recent Drude parameters may still require some tuning. Further studies could focus on this aspect. We also recognize that the TIP3P water model employed in the additive simulations can play an important role, and may contribute to the observed errors. SWM4-NDP is a four-point model, which can also adapt its dipole moment in reponse to perturbations of nearby charges. Moreover, the possibility that a Drude FEP using SWM6-NDP could produce better results is an interesting open question, since SWM6 is recognized to reproduce water properties even better than SWM4 (Yu et al. Other methods like replica exchange [START_REF] Kokubo | Two-dimensional replica-exchange method for predicting protein-ligand binding structures[END_REF]) increase sampling by thermally activating some or all degrees of freedom. Here, we used standard MD with long simulations, and we verified convergence by doing redundant calculations. As anticipated before in the text, the closure errors suggest that the C36 and Drude FEP were well converged, even if this result does not give information about the convergence of the µ A , µ B simulations. Indeed, these last contributions cancel out in the thermodynamic cycle since they were added and then substracted. Overall, the accuracy obtained indicates FEP can be used in a design strategy, where many variants obtained with empirical energy functions are filtered further with more rigorous FEP simulations before being tested experimentally. For ionic mutations in buried regions, electronic polarization is likely to contribute and the use of the Drude force field could be preferred. ). Syndecan1 is a member of the syndecans family, high conserved transmembrane and cytoplasmic domains containing four conserved tyrosine residues that may be phosphorylated [START_REF] Reiland | Pervanadate activation of intracellular kinases leads to tyrosine phosphorylation and shedding of syndecan-1[END_REF]).

In the case of Syndecan1, the phosphorylation of a tyrosine at C-terminal position P -1 has an important role in cell adhesion (Sulka et al. [2009]).

The Tiam1 PDZ domain can bind both the wildtype Syndecan1 peptide (Sdc1, sequence TKQEEFYA) and its phosphorylated form (pSdc1, sequence TKQEEF(pTyr)A). Their affinities to Tiam1 were studied experimentally, revealing that the phosphorylation of Tyr -1 in Sdc1 does not change considerably its binding affinity to Tiam1 ∆∆G bind (Sdc1 → pSdc1) = -0.21 kcal/mol (Liu et al. [2013]). Crystal structures of the wildtype Tiam1:Sdc1 complex (pdb access code 4GVD) and of the Tiam1:pSdc1 complex (pdb access code 4GVC) show that there is no Tiam1 rearrangement after binding, with an rmsd between the two backbone conformations of just 0.26 Å. However, the interactions between amino acid side chains in the binding pocket are different [fig. 5.1]. The tyrosine at position P -1 rotates of almost 90 • after phosphorylation, where the phosphate group interacs strongly with the Tiam1 residue K879 and also with T857 [fig. 5.1]. On the other hand, in the Tiam1:Sdc1 complex these interactions are different with Tyr -1 interacting with N876, E -4 and K -6 thorugh an hydrogen bonding network. Previous molecular dynamics simulations of the Tiam1:Sdc1 and Tiam1:pSdc1 complex (Nicolas Panel, results are not yet published) suggested that in the bound state the dianionic phosphotyrosine pTyr 2-is dominant respect to its monoanionic form pTyr -. During MD, monoanionic phosphotyrosine did not conserved the crystal structure configuration, moving towards the tyrosine orientation in the WT complex. On the other hand, dianionic phosphotyrosine conserved the crystal structure orientation, interacting strongly with K879.

Phosphates in biology

Phosphate groups are essential in biochemistry. They are components of nucleic acids and also present in certain protein side chains, whose interactions with solvent, metal ions and ionic side chains help control protein folding and binding (Westheimer [1987]). Phosphates are in ATP, they are added/removed enzymatically as a function of cellular condition, and they facilitate important PPIs in a number of signaling pathways (Pawson [1995]). Often, the presence or absence of the phosphate group makes the protein competent or incompetent to bind another molecule [START_REF] Shepherd | Distinct ligand specificity of the tiam1 and tiam2 PDZ domains[END_REF]). Indeed, the phosphorylated species act as chemical entities different from natural amino acids, making intramolecular and intermolecular hydrogen bonds with PPI partners. During evolution, nature chose phosphates because they were abundant on Earth, because they are solubile in water, and because they have versatile chemical properties (Hunter [2012]). Quoting the author, "Life as we know it could not have Computer simulations can give detailed insights into the structure and interactions of phosphate groups. To exploit this capability it is essential to have accurate force field models of these moieties. However, developing accurate simulation models is challenging for ionic interactions such as those involving phosphate groups.

Earlier additive results, need for polarizability

Phosphate-containing molecules will usually be more polarized when exposed to solvent or close to an ion, and less polarized when removed from solvent or ions. For this reason, additive force fields can give large errors, since they are usually parameterized to describe ionic interactions between groups that are highly-exposed to aqueous solvent. Previous calculations of the standard binding free energy between Mg 2+ and hydrogen phosphate P 2- i = HPO 4 were done using the Charmm additive force field C27, overstimating of an order of magnitude the experimental value of -3.7 kcal.mol [START_REF] Satpati | Free energy simulations of a gtpase: Gtp and gdp binding to archaeal initiation factor 2[END_REF]). Other previous calculations of the relative binding free energy for the phosphorylation of Tyr -1 in the Tiam1:Sdc1 complex (not published) revealed that the additive force field Amber ff99SB strongly overstimates the stability of pTyr in the bound state, giving a ∆∆G(Sdc1 → pSdc1) = -12.1 kcal.mol, where the experimental value is -0.2 kcal/mol.

If additive force fields give large deviations from experiments, remarkable results were obtained using polarizable models. Mg 2+ :phosphate binding free energies were studied with AMOEBA, giving excellent agreement with experiment [START_REF] Kumar | Structure and thermodynamics of mg:phosphate interactions in water: A simulation study[END_REF][START_REF] Kumar | Corrigendum: Structure and thermodynamics of mg:phosphate interactions in water: A simulation study[END_REF]). The AMOEBA force field was also used to reproduce the phosphate binding mode of a phosphatebinding protein, giving binding affinities in excellent agreement with experimental measurements [START_REF] Qi | Elucidating the phosphate binding mode of phosphate-binding protein: The critical effect of buffer solution[END_REF]). Dimethyl phosphate (DMP) and its interactions with Mg 2+ were studied with the Drude polarizable force field, giving good agreement with experiment for structure and dynamics (Lemkul & MacKerell [2016a]). Another recent study revealed that when compared to additive force fields, the Drude polarizable force field better estimates relative binding free energies for a series of ionic mutations of the Tiam1:Sdc1 complex [START_REF] Panel | Accurate PDZ/peptide binding specificity with additive and polarizable free energy simulations[END_REF]). All these studies confirm that to describe phosphate interactions, the explicit treatment of electronic polarization becomes necessary.

Methyl phosphate as a model

Here, we extended the Drude polarizable force field to three other phosphate compounds and their interactions with Mg 2+ , then we developed a model for the dianionic phosphotyrosine residue (pTyr 2-) [fig. 5.2]. We fist developed Drude parameters for methyl phosphate (MP), monoanionic and dianionc MP -≡ CH 4 PO - 4 , MP 2-≡ CH 3 PO 2- 4 . If DMP is analogous to the phosphate groups within DNA and RNA polymers, methyl phosphate (MP) is analogous to the phosphate groups at DNA and RNA termini. We validated our MP model computing MP:Mg 2+ standard binding free energies, for which the additive force field gave large errors. We also developed a model for mono hydrogen phosphate P 2- i = HPO 4 -, for which experimental Mg 2+

binding data is available. After validation, MP was used to derive parameters for dianionic phosphotyrosine, for which the additive force field overstimated its interactions in the binding pocket of the Tiam1 PDZ domain. Our parametrization strategy was inspired by the method used to develop the pTyr 2-model for the existing Charmm C36 force field (Feng et al. [1996]), where dianionic methyl phosphate was used to optimize a residue patch that transforms the standard tyrosine to pTyr 2-. The same approach could be used in the future to develop a Drude polarizable model for monoanionic phosphotyrosine (pTyr -) but also for other phosphatecontaining residues, like phosphoserine and phosphothreonine. Several programs and tools were used to obtain and to fit the QM data. We finally validated our models computing Mg 2+ :phosphate standard binding free energies in aqueous solution.

QM optimized geometry

Electrostatic parameters optimization QM electrostatic potential maps QM polarizability and dipole moment QM solute:water/ions interaction energy scans

Internal parameters optimization

QM dihedral energy scans QM conformational energies Phosphate:Mg 2+ binding free energies ). GAUSSIAN is the program traditionally used to compute target QM quantities for Charmm C36 and Drude force field parametrization. Several scripts to produce the QM data were already available, suggesting us to use the GAUSSIAN to compute most of the quantities. On the other hand, PSI4 is more recent and is becoming more and more used.

It is free and has a flexible python interface. Some of the GAUSSIAN scripts were so rewritten for Psi4, and will be part of our parametrization toolbox.

The GAAMP and DGENFF tools

Parameter optimization was done using several tools. Some of them are programs borrowed from the GAAMP distribution (Huang & Roux [2013]). GAAMP "General Automated Atomic Model Parameterization" is a webserver for automatic force-field parametrization of small molecules. It has a simple interface, allowing the user to load an input structure and then automatically obtain the desired parameters for various additive force fields (for example Charmm C36 or Amber) and for the Drude polarizable force field. However, we preferred to run a more careful manual parametrization, where several utilities written in C++ were extracted from the GAAMP server and then run locally on our computers. Other steps of the parametrization procedure were carried out computing molecular-mechanics quantities with the CHARMM program (version 42b2), together with bash and python scripts. CHARMM scripts were extracted from an ongoing first release of the Drude "DGENFF", analogous to the well-known toolbox for the Charmm force field parametrization CGENFF (Vanommeslaeghe

et al. [2009]).

Parametrization strategy

Starting from a PDB file containing structure information we generated the initial Drude topology and parameter files using GAAMP. Input geometries were optimized with GAUSSIAN, and then used to compute the target QM data (see Methods: QM calculations). .4, left), proceeding in two main phases. In the first phase, electrostatic potential maps and polarizability from QM were fitted with the MM energy function by adjusting the atomic charges, polarizabilities and Thole factors. In a second phase, these optimized parameters were further adjusted in order to reproduce the QM solute:water interaction energies, together with the QM molecular dipole moment. For details about the parameter fitting procedure, see the section "Methods: fitting the QM quantities". Once electrostatic parameters of MP -were obtained, we optimized those of MP 2-(2D structure shown in figure 5.4, center). We fitted the MP 2-atomic charges to all the QM data at once (electrostatic potential maps, water scans, molecular dipole moment and polarizability) by means of a grid search in MM parameter space. Once electrostatic parameters of MP -and MP 2-were obtained, we transferred them to P 2- i = HPO 2- 4 (2D structure shown in figure 5.4, right) and then reoptimized with a grid 96 2. Methods: overview search.

In this study we did not optimized LJ parameters. Indeed, most atom types were taken from the existing dimethyl phosphate (DMP) in the nucleic acid Drude force field; the main exception was an oxygen (O1) typed according to the lipid Drude force field. We also studied DMP with QM and with the existing Drude force field, to better evaluate the quality of our fit for methyl phosphate. For both MP -and MP 2-, the overall QM/MM agreement was similar to that obtained for DMP.

For both MP -and MP 2-, we then went on and combined the optimized electrostatic parameters with existing Drude bonded parameters. All steps stayed as close as possible to existing Drude types/parameters, with input also from C36 where relevant (eg, a few bonded parameters). Few dihedrals were scanned with QM. The total energy for several dihedral configurations was computed with QM and then compared to MM. In this phase, few parameters were optimized to improve agreement. At the end of the parametrization, we tested the overall MM model further by comparing MM conformational geometries and energies to QM.

Final MP and P 2- i models were validated computing their standard binding free energies with Mg 2+ . The magnesium binding model is described below. Affinities were computed with alchemical free energy simulations (see Methods: phosphate:Mg 2+ binding free energies). Dianionic phosphotyrosine was parametrized starting from existing Drude force field information and the MP 2-final set of parameters. The pTyr 2-side chain was modeled using a model compound, the phosphorylated form of the p-Cresol (see Results: Drude model for dianionic phosphotyrosine). QM ESP maps (unperturbed+perturbed) and solute-water interaction energies were fitted. Phosphotyrosine parametrization is still underway, but the model was partially validated running polarizable molecular dynamics of the Tiam1:pSdc1 complex. The MP -/2-: Mg 2+ binding is governed by the equation:

Phosphate:Mg 2+ binding model

(Mg 2+ ) w + (MP -/2-) w ∆G --→ bind (Mg 2+ MP -/2-) w (5.1)
where the subscript w is used to indicate that the binding reaction takes place in pure water. For P 2- i , the equation is similar. There are two possible bound states: Outer Shell (from now called OS) and Internal Shell (from now called IS) shown in figure 5.5. In IS, the phosphate interacts directly with magnesium ion, while in OS this interaction is mediated by a shell of water molecules. Experimental affinities for anionic and dianionic phosphates are known for H 2 PO - 4 (-1.7 kcal/mol) and P 2- i = HPO 2- 4 (-3.7 kcal/mol) (Verbeeck et al. [1984], Alberty & Goldberg [1992]). To obtain the binding free energies, we have run a series of alchemical MD simulations. The detailed protocol is explained below in the text (see Methods: phosphate:Mg 2+ binding free energies).

This study is motivated by two facts. First, methyl phosphate parameters were built using the established parametrization procedure, where most of the calculations involved constrained geometries and interactions in vacuum. To better validate our model it is important to check its behaviour during molecular dynamics in SWM4-NDP water. Binding affinities can give important insights about reiability of the developed parameters, and how this simple polarizable model can reproduce biochemical properties. Second, despite several phosphatecontaining moieties were already parametrized in the official Drude force field, the correct threatment of their interactions with charged ions are still challenging. In previous studies, a strong effort has been made to develop consistent Drude force field parameters to balance 3. Methods: QM calculations interactions between Mg 2+ and nucleic acids (Lemkul & MacKerell [2016a]). Several ad-hoc parameters were used to optimize interaction energies and geometries of complexes involving nucleic acid moieties, water and magnesium. However, they were not extensively tested. Magnesium binding free energy simulations are thus a stringent test not just for our new models, but also for the existing Drude force field parameters.

3 Methods: QM calculations Quantum calculations were performed using as reference the QM optimized geometry. Optimization was done with GAUSSIAN, using MP2/6-31G*. We extracted the optimized geometry from the GAUSSIAN output and then prepared several input files, used to calculate the following QM data.

QM electrostatic potential maps

Ten concentric grids based on Connolly surfaces were positioned around each solute using the program CGRID from GAAMP, totalling around 2000 grid points and ranging up to about 5 Å from the solute atoms. An example of such a grid is given in figure 5.6, left. GAUSSIAN allows to automatically generate an uniformly distribute cubic grid with the molecule centered in it, however the present strategy allow to reduce the number of points and to better adapt to the solvent accessible surface of different molecules. Points which are far from the center are inserted with a lower density factor respect to points distributed over closer Connelly surfaces. 67 perturbing charges were distributed around MP -and 69 around MP 2-[fig. 5.6, right]. Grid layers were alterned with perturbing charges layers. The electrostatic potential was then evaluated for each grid point i and for each perturbing charge j, totalling 67(69) perturbed ESP maps i φ QM ji for i = 1, ..., n grid plus one unperturbed map i φ QM i . During the QM calculation, the model compound was kept fixed at its optimized geometry, which is not relaxed for different perturbations. As grid points, perturbing charges were distributed uniformly on connelly surfaces, except some of them which were placed at specific positions following the protocol described by Anisimov et al. [2005]. A first group of charges sits along molecule's bond lines, at the intersection with the surfaces. Then a second group of charges was positioned around oxygen or nitrogen atoms when they form covalent bonds with two other atoms. Finally, the charge layers were covered respecting a minimum distance threshold between couples of charges of 1.5 Å. The electrostatic potential was evaluated with GAUSSIAN, using MP2/6-31G*.

QM polarizability and dipole moment

We computed molecular polarizability of MP -, MP 2-and P 2- i using MP2/6-31G*. The "exact polarizability" (au 3 ) was extracted from the GAUSSIAN output and converted to Å 3 (1 au = 0.529177249 Å). The dipole moment was computed with PSI4 using RIMP2/6-31G*. Only diagonal elements of the 3x3 polarizability matrix were considered because they are considerably greater than off-diagonal elements, which were considered negligible.

QM solute-water and solute-ion interaction energy scans

We start describing the protocol for water scans. Water poses were generated with GAUSSIAN (scan keyword) or with a CHARMM script, varying only the solute-water distance in the [1.4:3.0] Å range, in 0.1 Å steps. During each scan, the solute and water internal geometries were held fixed. Internal coordinates were specified in a Z-matrix and the solute-water distance was controlled by an r OH variable. For each pose, the interaction energy was computed with PSI4 using RIMP2/cc-pVQZ with BSSE corrections. To define water poses relative to a specific acceptor/donor, a DUMMY particle was added to the solute-water system [fig. 5.7, left]. This auxiliary atom was then used to specify angles and dihedrals that define the internal geometry of the compound:water complex. In the figure, the water hydrogen H1, the phosphate oxygen O2 and the dummy atom DUM form together an angle of 90 • . Solute-ion scans for magnesium (Mg 2+ ) and sodium (Na + ) were performed following an equivalent protocol, except that QM calculations with PSI4 were done using MP2/aug-cc-pvtz with BSSE corrections.

MP -has four hydrogen bond acceptors and one donor: O1, 02, O3, O4, and H2. Although O3 and O4 have equivalent types, they are not symmetrical in the optimized structure so we did a separate scan for each atom. For MP 2-, we scanned the three acceptors: O1, O2, O3 (equivalent to O4), while for P 2- i we scanned acceptors O1, O2 and O3. Interactions were also scanned for the bridging oxygen of dimethyl phosphate (DMP).

QM dihedral scans

Starting from the QM optimized geometry, selected dihedrals were scanned using angles ψ ∈ [-π : π] with a step of 15 • . We scanned O1-P1-O2-H2 and C1-O1-P1-O2 dihedrals for MP -and the C1-O1-P1-O2 dihedral for MP 2-. Scans were performed with GAUSSIAN (scan keyword) or with PSI4, using MP2/6-31+G(d). In the GAUSSIAN input file, geometries were defined by a z-matrix where the target dihedral angle a variable. For each configuration visited during the scan, the system energy was computed after geometry relaxation, performed keeping the scanned dihedral fixed at its value φ. The QM energy and geometries were saved for eachi value of ψ, to be compared afterwards with MM.

QM vibration frequencies

We performed vibrational analysis using Gaussian with the freq keyword. We used MP2/augcc-pVTZ for geometry optimization and frequencies calculation. From the QM output we extracted frequences and vectors for all the normal modes and also force constants in cartesian coordinates, represented in matrix form (the Hessian). This 3N x3N matrix has entries H ijkl where i,j = 1, ..., N is are atom IDs and k,l = x, y ,z.

H ijkl = ∂U QM ∂k∂l (5.2)
The QM force constants matrices of MP -and MP 2-were used to compare QM and MM energies for random displacements generated by combination of the extracted QM modes. Indeed, for small coordinate displacements i d i (where i = 1,...,Natoms) around the molecule's optimized geometry, the QM energy can be estimated from the force constants matrix using the harmonic approximation

Ũ QM = 1 2 ijkl d ik H ijkl d jl (5.3)
where i,j are atoms IDs and k,l = x, y ,z.

4 Methods: fitting the QM quantities

Fitting QM potential maps

The QM electrostatic potential maps were fitted using the GAAMP "drude-fitcharge" module. In addition to the above constraints on the molecule geometry, restraints were applied to charges, polarizabilites and Thole parameters to avoid large deviations from physically reasonable reference values. The target function χ 2 is composed of several terms:

χ 2 = χ 2 ESP + χ 2 CG + χ 2 ALP HA + χ 2 T HOLE
(5.4)

The first term tries to match the QM and MM potential maps. In the presence of n grid points and n pert perturbation charges,

χ 2 ESP = 10 4 2n pert n grid   n pert ngrid i=1 (φ M M i -φ QM i ) 2 + npert j=1 n grid i=1 (φ M M ji -φ QM ji ) 2   (5.5)
4. Methods: fitting the QM quantities

The other terms are the restraints, weighted by the constant factors w CG , w α and w τ (set to default values). Reference values for charges are those specified in the initial topology file {q 0 i }:

χ 2 CG = w CG n n i=1 f (q i ,q 0 i ) (5.6)
where n is the number of atoms and

f (q i ,q 0 i ) =      0 |q i -q 0 i | ≤ 0.03 (|q i -q 0 i | -0.03) 2 otherwise (5.7)
Polarizabilities α i and Thole parameters τ i were restrained using

χ 2 ALP HA = w α n np i α i -α 0 i 2 χ 2 T HOLE = w τ n np i τ i -τ 0 i 2 (5.8)
where n is the number of polarizable atoms. Reference polarizabilities α 0 i were taken from [START_REF] Miller | Additivity methods in molecular polarizability[END_REF]) then scaled by 0.7. Reference Thole factors τ 0 i were set to 1.3 for all atoms; this corresponds to average screening parameter τ ij = τ i + τ j = 2.6, as originally proposed by Thole.

Fitting the molecular polarizability

Using GAUSSIAN, one obtains α in the {xyz} frame. Its diagonal components are denoted α xx , α yy , α zz . Atomic polarizabilities and Thole factors were manually adjusted to reproduce the QM values of α xx , α yy , α zz (and thus their sum, the polarizability trace). The atomic charges were held fixed (since α is independent of the charges). During adjustment, we crosschecked that the polarizability and Thole modifications did not affect too much the quality of the potential maps (the value of χ 2 ESP ). The Drude molecular polarizability was obtained using the CHARMM script polar_efield (from DGENFF). Diagonal elements of αmol in the xyz frame were computed by measuring the variation of the molecular dipole moment under the influence of an external electric field of intensity E = 10 10 V/m. When measuring the electronic polarizability, the atoms were kept fixed in the molecule's optimized geometry. Drude particles were relaxed before computing the dipole moment in the presence or absence of the applied field. The field was applied using the CHARMM command pull efield 1E10 xdir <float> ydir <float> zdir <float>

The molecular mechanics calculation of α xx , α yy , α zz was done in 4 steps. First, the molecular dipole moment μ0 was computed in the absence of an external field, with the CHARMM command "coor dipole oxyz". Next, values of the molecular dipole moment μx , μy , μz were computed in the presence of three different external fields, oriented along x, y, and z, respectively. The nine components of the molecular polarizability tensor were then obtained using

α ij = µ i j -µ 0 j E (5.9)
4.3 Fitting solute-water interaction energies and the molecular dipole Two strategies were used. For MP -, we adjusted the atomic charges manually, without changing the polarizabilites or Thole values. For one atom (the bridging oxygen O1) we departed from the Drude DMPN vdW types (see Results). After adjusting the charges, we checked the quality of the potential maps, the dipole and the polarizability. For MP 2-, we optimized the atomic charges to fit all the QM data at once, through a simple grid search in parameter space. Starting from MP -parameters, we transferred atomic polarizabilities and thole factors to MP 2-. Atomic charges were readjusted to match the total charge. Then we fitted molecular polarizability (not dependent on atomic charges), adjusting polarizability and Thole factors manually. As for MP -, targeted molecular polarizabilities were scaled by a factor 0.85. Once α and τ factors were obtained, we went on performing a grid search for many values of atomic charges to fit water interaction energy scans for acceptors O1, O2, O3. An equivalent protocol was used to derive electrostatic parameters of P 2- i .

5 Methods: phosphate:Mg 2+ binding free energies

Free energy perturbation protocol

To obtain the MP -/2-:Mg 2+ and P 2- i :Mg 2+ binding free energies, we simulated the alchemical annihilation or "decoupling" of magnesium alone in solution and bound to the phosphate solute, over a series of MD simulations. For MP -and MP 2-, we considered both IS and OS binding. For P 2- i , we considered OS binding. From each simulation, the binding free energy was obtained using the Double Decoupling Method (DDM). Here we describe the procedure for methyl phosphate, given that for P 2- i we used exactly the same protocol. The thermodynamic cycle is displayed in figure [5.8]. The starting point was a standard-state ideal-dilute solution of the unbound MP + Mg 2+ at concentration ρ 0 = 1M , the endpoint was a standard-state ideal-dilute solution of the complex (MP -/2-Mg 2+ ) (bound state). The reaction is decomposed in several steps labelled with integers from 1 to 6 and described in the text.

The binding reaction was decomposed in six steps, connecting starting and ending states: [1] We changed the concentration to the very low value ρ ∞ . The corresponding free energy change was:

∆G(ρ 0 → ρ ∞ ) = -2KT lnρ ∞ /ρ 0 (5.10) [2]
We removed the Mg 2+ charge and LJ interactions alchemically, in two steps. First, magnesium charge and polarizability were turned off. The corresponding free energy difference was ∆G elec . Then, Lennard-Jones interactions of the uncharged Mg were removed giving a second contribution ∆G LJ .

[3] A flat-bottomed harmonic restraint was introduced to keep Mg 2+ close to MP. The re-simulations (point 4) the solute was represented by the MP-Mg 2+ complex. OS and IS states were prepared separetly, putting the Mg 2+ ion at average distance d=(a+b)/2 from the P atom, where a and b vary between OS and IS. Each solute was then solvated in a cubic box of side L=31 Å, containing around 1000 pre-equilibrated SWMD4-NDP water molecules. Solvent within the solute volume was deleted using a distance threshold of 3.0 Å from the solute heavy atoms. The system was then minimized with NAMD, first with restrained heavy atoms and unrestrained Drude particles, then without restraints. The minimized system was finally equilibrated, running 200 ps of MD.

All MD similations (equilibration and alchemical) were carried out using NAMD (version 2.12) at room temperature and pressure, using Langevin dynamics with a Langevin Piston Nosé-Hoover barostat (Feller et al. [1995]). The simulations used periodic boundary conditions and a Particle Mesh Ewald electrostatic treatment [START_REF] Darden | Computational biochemistry and biophysics, treatment of long-range forces and potential[END_REF]). Bonds to hydrogen atoms were constrained with the SHAKE method [START_REF] Ryckaert | Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes[END_REF]). LJ interactions were shifted to zero for separations between 10 and 12 Å. The MD timestep was 1 fs and the oscillators temperature was T D =1 K. Throughout the simulations, the center of mass of MP (or P 2- i ) was restrained to the box center using an harmonic potential.

Alchemical free energy simulations

Alchemical free energy changes for steps 2 and 4 were extracted with TI [eqn. 3.51] from series of MD simulations. Decoupling of electrostatic interactions of Mg 2+ in the bound or unbound state was performed by scaling its charges (the atom charge q A and its Drude charge q D ), its polarizability α and Thole factor τ . MD simulations were run using the potential function of equation 4.2, for different values of the coupling parameter λ elec =1, 0.9, 0.8, ..., 0.1, 0. Decoupling of LJ interactions was done using the alch facility of NAMD, with a series of simulations using λ LJ =1, 0.9, 0.8, 0.6, 0.4, 0.2, 0. The electrostatic uncharging(respectively, charging) free energies at concentration ρ ∞ were extrapolated analitically for L → ∞, using equation 3.71

∆G(L → ∞) = ∆G(L) + q 2 |ζ| 2ǫ W L + 2πq 2 R 2 3L 3 ǫ W -1 ǫ W + O(R 2 L -3 ǫ -3 W ) (5.14)
where ǫ W is the solvent dielectric constant, L is the box side length in the NPT simulation (on average 31.2Å) and ζ = -2.837297. The correction was computed considering two distinct values of Mg 2+ radius for bound and unbound states, respectively R=1.2 Å and R=1.5 Å.

6 Results

Initial MP model and atom types

Starting topology and parameter files for MP -and MP 2-borrowed as much information as possible from the Drude dimethyl phosphate (DMPN) of the official Drude nucleic acids topology files. We give existing C36 and Drude information first.

-----------------c36 methyl phosphate information------------------------------------------------------------------------------- ---------Drude dimethyl phosphate, from nucleic acids----------------------------------RESI DMPN - The initial MP -is given below, where atom O2 (the hydroxyl oxygen) has a Lone Pair. Atom types of the initial Drude model were transferred from "nucleic acid" DMPN.

---------Initial Drude methyl phosphate MP- -------------------------------------------

Results

We start giving results for the MP -parametrization, describing in detail all the intermediate sets of parameters. Results for MP 2-and P 2- i will be given later, as they were obtained all at once with a grid search. MP -parametrization was carried out following three steps executed in sequence, resumed in table 5.1.

Set Construction

Target data 1 adjusted q i , α i , τ i to minimize χ 2 with GAAMP QM electrostatic potential maps 2 charges from set 1, manually adjusted α i and τ i QM molecular polarizability 3 α i and τ i from set 2, manually adjusted charges QM dipole moment and water scans 5.2 -Electrostatics Parameters for MP -. Some of the parameters were fixed, other are redundant for atoms belonging to equivalent groups, as described in the text.

In table 5.2 we indicate the list of MP -electrostatic parameters to optimize. We imposed two groups of equivalent atoms: (H11, H12, H13) and (O3, O4). Each group shares the same charge, polarizability and Thole parameters. For the restrained electrostatic potential (RESP) fit, we further reduced the number of parameters by fixing some of the charges at their DMPN values: q(H11)=q(H12)=q(H13)=0.026 and q(C1)=0.0725. The charge of atom O1 is a function of all the other charges. Initial Thole factors were set to 1.3 and weights for restraints were w CG = 3.0, w α = 0.4, w τ = 0.2 (the defaults). The performance and parameters after fitting (set 1) are given below in

Fitting the MP -polarizability

We first computed the diagonal elements of the polarizability tensor using parameter set 1 (post-ESP fit). All the components were understimated compared to the target value, which is the QM result scaled by 0.85. We manually adjusted polarizabilities and Thole parameters. We found that α was well reproduced by transferring parameters from DMPN and modifying only the O1 Thole value. Fit results are in table 5

Fitting the MP --water radial interaction energy scans

Starting from the set 2 parameters, we manually modified charges to better reproduce energy scans and also the molecular dipole moment. During adjustment, we periodically checked the quality of the RESP fit. Below, we report QM (target) and MM interaction energies e min at the minimum distance r min for the MP -acceptor and donor atoms [tab. 5.7]. We also report the optimized charges, "set 3" [tab. 5.9]. The solute-water interaction for acceptor O1 was poorly fitted initially. The problem was related to the O1 type/vdW parameters we had transferred from the existing Drude atom type OD30BN. Indeed, the results suggest that the vdW parameters of atom type OD30BN (Drude FF, nucleic acids) are not ideal for the MP --water interaction when scanning the acceptor O1. Therefore, we studied the same interaction in DMP -, which carries the same atom type on its bridging oxygen.

Results

Atom

We optimized the DMP geometry with Gaussian using HF/6-31G* and then we scanned the water interaction energy, focussing on the the same acceptor O11 (equivalent to O1 in MP -). Interaction energies were obtained with PSI4 using RIMP2/cc-pVQZ with BSSE corrections. Computing interaction energies with the MM energy function and the existing Drude DMPN parameters, we obtained the same error: the interaction e min was too low and r min was strongly underestimated. Minimum energy at minimum distance were more like those of MP -r min = 1.5 Å and e min = -11.544 kcal/mol, where QM values for DMPN were r min = 1.9 Å and e min = -10.996 kcal/mol. Results are shown in figure 5.9. At this point, we repeated the MM scan using a different set of parameters: we assigned the vdW parameters of a different atom type OD30B (Drude FF, lipids) and we transferred to DMP the electrostatic parameters obtained at the end of the MP -optimization (set 3). This transfer required a small charge adjustment to respect the total charge of DMP -. With this substitution, the error was completely removed and the MM radial scan fit well the QM target data. Values were r min = 1.8 Å and e min = -10.784 kcal/mol, more like those obtained with QM. At this point, we replaced the atom type OD30BN with the atom type OD30B also for MP -. Then repeated the charge adjustment. With this new set of charges, the radial interaction energy scans were better fitted for O1, O3, O4 and the donor H, while for O2, the result was slightly worse but still acceptable. To better evaluate the quality of the final set of parameters, the resulting fit for MP -will be given below, together with results for MP 2-and P 2- i .

6.5 Electrostatic parameters optimization for MP 2-and P 2- i Parameters for dianionic methyl phosphate (MP 2-) and hydrogen phosphate P 2- i were obtained performing a grid search in parameter space. Initial models were created transfering MP -atom types and electrostatic parameters, rapidly adjusted to match the total charge. Polarizabilities and Thole factors were manually adjusted to reproduce the QM polarizability. Then the charges were optimized performing a grid search, that fits all QM data at once. The full QM/Drude comparisions are shown below MP 2-and P 2- i , including MP -data for comparision. In table 5.10, we give the fit of molecular polarizabilities and dipole moments. In table 5.11 we give the result of the fitted water scans. The configurations used for the scans are given in figure 5.10 for MP -(MP 2-and P 2- i water poses were similar). In figures 5.11 and 5.12 we plot all the scans for MP -and MP Figure 5.12 -MP 2 --water interaction Grey: QM; Black: MM

Fitting the dihedral energy scans

Electrostatic parameters to compute the MM energy are those optimized above. Internal parameters are those of DMPN except for the hydroxyl torsion (no present in the existing Drude force field) which has been borrowed from C36 and then reoptimized. In all scans, the molecule geometry energy was relaxed at each value of ψ in the selected interval, both in QM and MM calculations. More precisely, the scan has been performed first using QM. For each dihedral value, the QM relaxed conformation has been saved and then re-minimized using MM. We obtained perfect agreement for MP 2-and small deviations for MP -dihedrals. Scans (QM vs MM) are given in figure 5.13. Local minima differs for ∼ 1 kcal/mol and energy barriers are well reproduced using the MM energy function.

Conformational energies of MP -and MP 2-

For both MP -and MP 2-we used optimized QM geometries and the force constant matrix in cartesian coordinates [eqn. 5.2] to estimate the total QM energy Ũ QM [eqn. 5.3] for 200 random conformations. The obtained energies were then compared to the ones calculated with the MM energy function applied on the same 200 random conformations. For each conformation, displacements d i were obtained as random combination of the QM normal modes, extracted from the Gaussian output. We used all modes except the trivial ones (rotations, translations of the molecule) and the ones relative to hydrogen(s) stretching.

Results

Figure 5.15 -Ion scan with magnesium performed on dimethyl phosphate (DMPN). Two positions were considered: direct interaction (Yellow) and inverse direction (Blue)

Results show that the NBFIX/NBTHOLE parameters of table 5.12 were optimized to reproduce the inverse scan for the bridging oxygen, while they produce poor agreement with QM for the "direct" scan of the same interaction. This was a reasonable choice, since the inverse interaction is more favorable compared to the direct one. Moreover, minimum energies for both directions were not well reproduced. On the other hand, using the lipid atom type the energy is better fitted, even if the minimum distance could be improved. O1 scans shows improvement also for MP. This suggested us to use the lipid atom type instead of the NBFIX.

To better check model transferability we also computed scans for sodium, using the parameters developed by [START_REF] Savelyev | All-atom polarizable force field for DNA based on the classical drude oscillator model[END_REF]

Mg 2+ :phosphate binding free energies

Experimental Mg 2+ binding free energies are available for HPO - 2 and H2PO -but not for methyl phosphate. Nevertheless, we expect that the methyl substitution has a modest effect on the free energy, since the magnesium ion binds on the opposite side of the molecule, as shown in figure 5.5 where IS and OS binding modes were extracted from snapshots of the MD simulations. This is further supported by the similar OS binding free energies computed below for MP 2-and P 2- i . The MP -IS complex mostly coordinated Mg 2+ through a single oxygen (mean Mg 2+ -O distance of 2.07±0.05 Å). The MP 2-IS complex mostly had bifurcated magnesium coordination, with two of the phosphate oxygens interacting closely with Mg 2+ (mean distances of 2.09±0.07 Å).

The IS complexes were stable during 1 ns of unrestrained MD, indicating that they are separated from OS by a free energy barrier. The MP 2-OS complex was also stable during 1 ns of MD, with no tendency to unbind or shift to IS. The MP -OS complex had a shorter lifetime, estimated from a 10 ns unrestrained MD simulation [fig. 5.16]. At the simulation concentration (55 mM), the MP -OS binding free energy was 0.4 kcal/mol (see below), giving an expected long-time limit of 33% for the OS population. The observed population during the 10 ns of MD was 27%, indicating that rough convergence was achieved for the binding/unbinding rates (which is sufficient for our purpose). The complex remained bound for intervals of about 275 ps on average, rebinding a few hundred picoseconds later. Evidently, the barrier separating the OS bound and unbound states is small, probably on the order of 2 kcal/mol. Transitions to IS were not observed. The OS lifetime is roughly comparable to the rate of OS/unbound transitions obtained for P - i with the AMOEBA force field (which was not measured precisely, [START_REF] Kumar | Structure and thermodynamics of mg:phosphate interactions in water: A simulation study[END_REF][START_REF] Kumar | Corrigendum: Structure and thermodynamics of mg:phosphate interactions in water: A simulation study[END_REF]).

For P 2- i , we only considered OS binding, for several reasons. OS binding was predominant for MP 2-in the present simulations, see below; we expect P 2- i to behave similarly. It also proved difficult in this work to obtain a good force field model for short-range P 2- i -Mg 2+ interactions using the same van der Waals parameters as for MP. More extensive parameter optimization with one or two new atom types would be necessary. Since modeling inorganic phosphate was not our goal in this work, we left a full parameter optimization and IS binding study for future work. The computed free energies and their components are given in table 5.15. The unbound Mg 2+ decoupling free energy, 412.3 kcal/mol, is 1.8 kcal/mol larger than the result (corrected for box size effects) of earlier simulations (Lemkul & MacKerell [2016b]) that used a smaller model and shorter runs. Repeating the entire decoupling calculation (51 ns) for the Mg 2+ -MP -complex gave the same result within 0.1 kcal/mol, indicating good convergence. For MP -, the computed binding free energies were -2.1 kcal/mol for the OS state, 0.0 kcal/mol for the IS state, and -2.1 kcal/mol overall. All three values are within 0.8 kcal/mol of the AMOEBA predictions for H2PO -. The overall binding free energy is 0.4 kcal/mol away from the experimental H2PO -value of -1.7 kcal/mol (Verbeeck et al. [1984], Alberty & Goldberg [1992]). Binding is predicted to be dominated by the OS state, which has a 97% occupancy, close to the AMOEBA H 2 PO -prediction (95% OS). For MP 2-, computed binding free energies were -2.8 kcal/mol for the OS state, -0.5 kcal/mol for the IS state, and -2.8 kcal/mol overall. The overall value is 0.9 kcal/mol away from the experimental HPO 2- value of -3.7 kcal/mol. The bound state is predicted to be 87% OS and 13% IS. Earlier ab initio calculations with a dielectric continuum solvent also predicted that OS was favored over IS, but with a smaller difference. For P 2- i , the OS binding free energy was -4.0 kcal/mol, close to the experimental value. The OS values for MP 2-and P 2- i show that the methyl substitution has a modest effect, as expected, weakening binding by 1.2 kcal/mol. We note that for MP 2-, P 2- i , and unbound Mg 2+ , the correction for the MD box size was small but distinctly nonnegligible (over 1 kcal/mol). It was zero for MP -, since introducing Mg 2+ changes the total charge from -1 to +1 (same absolute charge magnitude). Overall, for all three complexes considered, the computed binding free energies are in good agreement with experiment. This contrasts dramatically with the additive, Charmm C36 force field, which overestimated the Mg 2+ binding free energy by an order of magnitude [START_REF] Satpati | Free energy simulations of a gtpase: Gtp and gdp binding to archaeal initiation factor 2[END_REF]). The results are consistent with mostly OS binding for all the complexes. 7. Classical Drude model for dianionic phosphotyrosine

Results

We note that a couple of atom types differ respect to those described during optimization. Their names were changed to avoid overlapping parameters/types in the official force field, including those for interactions with Mg 2+ and the existing force field parameters for the nucleic acids bridging (OD30BN) and hydroxyl (OD31A) oxygens. The first was replaced by the new type OD30D, the second by the new atom type OD31D. They have the same LJ parameters of the existing OD30B (lipids) and OD31A (nucleic acids) types.

Classical Drude model for dianionic phosphotyrosine

As anticipated in the Methods section, after validation of the methyl phosphate and P i model we parametrized the dianionic phosphotyrosine pTyr 2-. The tyrosine residue in the existing Drude force field was derived from the model compound p-Cresol C 7 H 8 O (CRES, 3D structure shown in figure 5.18). The official Drude topology files for CRES and TYR are given in figure 5.17. Atom types and electrostatic parameters of CRES were completely transferred to the TYR side chain, then the charges of C β hydrogens HB1 and HB2 were readjusted. To parametrize dianionic phosphotyrosine, we adopted a similar strategy. We considered the dianionic, phosphorylated form of pCRESOL (residue CREP), 3D structure shown in

QM quantities

We computed target QM quantities, based on the QM optimized geometry. All calculations were run with the same QM softwares and methods used for MP. During optimization, the CZ-OH-P-O2 dihedral was constrained to 180 • . Then we computed unperturbed and perturbed QM ESP maps using a grid of 2178 points and 64 perturbing charges. We scanned the solute:water interaction energy for acceptors OH and O3, then for donors P and HE1. Geometries used for the scans are shown in figure 5.20. Using PSI4, we scanned the dihedral CE1-CZ-OH-P. During the scan we kept the CZ-OH-P-O2 dihedral constrained to 180 • , as for geometry optimization.

Electrostatic parameters optimization

Parameter optimization was done following an existing procedure, used in a previous study to develop the phosphotyrosine model of the additive force field C36 (Feng et al. [1996]). We considered pTyr like a patch for the tyrosine residue, which shares with TYR the atoms belonging to the methyl group and the first part of the TYR ring, as shown in figure 5. [START_REF] Findless | Design of artificial riboswitches as biosensors[END_REF]. Charges of atoms CE1(2), HE1(2), CZ, OH, P, O2(3,4) were readjusted performing a grid search in order to obtain the better agreement between Drude and the QM ESP maps and QM solute-water interaction energy scans. After optimization, we obtained a χ 2 ESP =0.25 in good agreement with the QM data. Results for the water scans are listed in table 5.16. We give topology information for the final, optimized CREP and pTyr 2-models in the extended format. Atom names are the ones used for the pTYr (TP2) residue of the C36 force field (Feng et al. [1996]). Starting from the crystal stucture, the Tiam1:pSdc1 complex was solvated in a large box of size L=80 Å containing pre-equilibrated TIP3P water. Then we deleted solvent molecules within the volume occupied by the solute, using a distance threshold of 3.0 Å from the complex heavy atoms. The box was trimmed, from cubic to truncated octahedron, deleting solvent molecules at the cube edges. Drude and lone particles were added with CHARMM, then the Drude particles were minimized with restrained atomic centers. The system was finally relaxed without restraints.

Equilibration MD was done using the Drude force field, in a step-wise procedure: first restraints were applied on the solute backbone and side chain atoms, then just on the solute backbone. After equilibration, we have run 20 ns of producion MD. Throughout the simulations, the center of mass of Tiam1 was restrained to the box center using an harmonic potential. MD was done using NAMD (version 2.12) at room temperature and pressure, using Langevin dynamics with a Langevin Piston Nosé-Hoover barostat (Feller et al. [1995]). We used periodic boundary conditions and a PME electrostatics [START_REF] Darden | Computational biochemistry and biophysics, treatment of long-range forces and potential[END_REF]). Bonds to hydrogen atoms were constrained with SHAKE [START_REF] Ryckaert | Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes[END_REF]). LJ interactions were shifted to zero for separations between 10 and 12 Å. MD was run using a timestep of 1 fs and Drude oscillators temperature T D =1 K.

We extracted configurations from the production simulation, where pTyr 2-strongly interacts with K879. Average and crystal structures are showed in figure 5.23. The average phosphorus-nitrogen distance during the 20 ns of MD was 3.42 Å, close to the 3.61 Å of the crystal structure. The average value of the χ 1 dihedral (N-CA-CB-CG) was 168.4 degrees, slightly bigger than in the crystal structure where it is 158.5 degrees. In the crystal structure, phosphotyrosine interacts also with the Tiam1 residue T857. In the simulation, this interaction is weaker, with an average phosphate-threonine side chain distance which is more than 3 Å larger respect to the crystal structure.

Conclusion

A classical Drude model for methyl phospate and its interactions with Mg 2+ was developed. The parameters obtained were similar to those derived earlier for DMP, Lemkul & MacKerell [2016a] with some significant differences. Agreement between the QM and MM data was similar to that obtained earlier for DMP.

The MP model was used to compute Mg 2+ standard binding free energies. Results were in good agreement with experiments, indicating that OS state was predominant for both anionic forms of MP. This is the first time the predominant bound state has been clearly identified for MP 2-and MP -. For this and similar binding processes, atomistic simulations are a powerful tool, and the polarizable force fields perform definetively better then additive additive force fields. This has now been confirmed for phosphate-Mg 2+ binding free energies, and was confirmed also for several other ionic interactions [START_REF] Panel | Accurate PDZ/peptide binding specificity with additive and polarizable free energy simulations[END_REF]).

The parametrization and simulation of MP with the Drude polarizable force field is a step towards modeling phosphorylated side chains in proteins. The Drude force field is applicable to

Conclusion

large biomolecules and provides the speed and accuracy for long and wellconverged simulations. The accuracy obtained here for phosphate-Mg2+ binding indicates that the Drude force field can be applied to elucidate molecular recognition processes involving phosphorylated proteins, which are essential for biological signalling.

Résumé

Les interactions protéine-protéine (IPPs) médient la signalisation cellulaire. Leur ingénierie peut fournir des informations et conduire au développement de molécules thérapeutiques. Les IPPs sont souvent médiés par domaines spécialisés, avec différents modes de reconnaissance (grandes surfaces de contact, peptides). Leur caractérisation avec des modèles physiques est complexe à cause des nombreux événements qui ont lieu après le binding: changements de conformation, réorganisation des molécules de solvant, redistribution des charges ou protonation-déprotonation des chaînes latérales.

Le travail de thèse est focalisé sur les domaines PDZ, importants médiateurs de IPPs. Elles lient les 4-10 résidus C-terminaux de protéines cibles. Elles lient aussi les peptides correspondants, qui peuvent servir de systèmes modèle ou d'inhibiteurs. Nous avons développé deux approches computationnelles basé sur des modèles physiques et les avons appliquées au domaine PDZ de la protéine Tiam1, un facteur d'échange pour la protéine Rac, impliqué dans la protrusion neuronale. Sa cible est la protéine Syndecan1. Des affinités expérimentales sont connues pour le peptide C-terminal, noté Sdc1, et plusieurs mutants; elles ont servi pour tester les calculs.

Dans une première phase, nous avons développé des modèles performants "haut débit" basé sur l'échantillonnage Monte Carlo adaptatif. Contrairement aux différentes méthodes computationnelles existantes, les nouveaux modèles permettent l'étude d'un grand nombre de variants de protéines (grands ensembles des séquences et des structures) et de réduire significativement les temps de calcul. Nous avons appliqué notre méthodes pour de dessin computationnel haut débit de Sdc1 dans le complexe Tiam1-Sdc1. Une simulation Monte Carlo est faite où les chaines latérales de la protéine et du peptide peuvent changer de conformères et certaines positions peuvent muter. Le solvant est implicite. Le paysage énergétique est applati par la méthode adaptative de Wang-Landau, de sorte qu'un vaste ensemble de variants est échantillonné. Effectuant des simulations distinctes du complexe et du peptide seul nous avons obtenu les énergies libres relatives d'association de 75,000 variants en heure CPU sur une machine de bureau. Les valeurs sont compatibles avec les quelques données expérimentales disponibles. Notre modèle haut débit est implémenté dans la nouvelle version de Proteus, le Protein-protein interactions (PPIs) mediate complex signaling networks in cells. Engineering them can provide understanding and is a recognized strategy for drug design. PDZ domains are among the most widespread domains mediating PPIs. They recognize the 4-10 Cterminal amino acids of their target proteins. They can also bind the corresponding peptides, which can thus serve as inhibitors or model systems. We have developed and tested two computational approaches to characterize and engineer PDZ:peptide recognition. We applied them to the PDZ domain of the Tiam1 protein, a Rac GTP exchange factor involved in neuronal protrusion and axon guidance. Its natural target protein is the Syndecan1 protein. Experimental affinities are available for the Cterminal Syndecan1 peptide, denoted Sdc1, and several mutants; this data was used for testing and benchmarking. We first developed a novel high throughput strategy for protein and peptide design. A Monte Carlo simulation was done where protein and peptide side chains could change conformations and selected positions could mutate. Solvent was modeled implicitly. The energy landscape was adaptively flattened, following the Wang-Landau method, so that a very diverse set of sequences was sampled. By performing separate simulations for the PDZ:peptide complex and the unbound peptide, we could recover the relative binding free energies of around 100,000 peptide variants using just one hour of CPU time on a desktop machine. The computed affinities were consistent with (sparse) available experimental data. Next, we developed a much more detailed and accurate simulation model. The solute and solvent were described in full atomic detail, using molecular dynamics simulations and a state-of-the-art force field that explicitly accounts for electronic polarization: the Charmm Drude force field. Polarizability is expected to be important for PDZ:peptide binding, since residues that form the binding interface are transferred from a solvent-rich to a solventpoor environment upon binding. We performed alchemical free energy simulations to compare the binding free energies of four peptide variants, which all differ in one or two ionic side chains. The calculations gave good agreement with experiment. In contrast, the Charmm and Amber additive force fields, which treat electronic polarization implicitly, gave poorer agreement. These calculations represent the first example of relative protein:ligand binding free energies computed with alchemical free energy simulations and a polarizable force field. Finally, to allow future modeling of phosphorylated peptide variants, we extended the Drude force field to include both methyl phosphate and phosphotyrosine. The force field gave excellent agreement with experimental phosphate--magnesium binding free energies.

  [2011], Steiner & Caflisch [2012], Blöchliger et al. [2015] and Sensoy & Weinstein [2015].

1. 1 Figure 1 . 1 -

 111 Figure 1.1 -Bonded interactions Bond (left), Angle (center) and Dihedral (right).

1 .Figure 1 . 2 -

 112 Figure 1.2 -Bonded interactions Urey-Bradley (left) and Improper (right).

Figure 1 . 3 -

 13 Figure 1.3 -Solvent accessible surface obtained rolling a sphere of radius r which represents a rigid water molecule in contact with the solute. For surface calculations in Proteus we use a probe radius r = 1.5Å.

Figure 1 . 5 -

 15 Figure 1.5 -GB interaction energy as function of the distance between a couple of charges r. Blue line: interaction for a couple of exposed charges, low b i b j . Green line: interaction for a couple of buried charges, high b i b j . Red line: interaction for intermediate value of b i b j

Figure 1 . 6 -

 16 Figure 1.6 -Surface overlap between three buried residues. Individual surfaces are colored in green, red, blue.

Toulbar2

  is a C++ solver for cost function networks (CFN) developed at INRA in Toulouse (Allouche et al. [2014]). The program uses a cost function to find the GMEC from a matrix containing interactions between couples of residues (Traoré et al. [2013]). These interactions, however, must be computed using an external program like Proteus or Rosetta-fixbb (fixed backbone design).

Figure 1 . 7 -

 17 Figure 1.7 -General Proteus CPD protocol of 3 steps, as described in the text

Figure 1 . 8 -

 18 Figure 1.8 -Unfolded and Folded states as described in the text. The unfolded state is modeled as an extended polymer where side chains belonging to different residues do not interact with pairwise terms.

Figure 1 . 9 -

 19 Figure 1.9 -Selected rotamers for side chains of ASP and GLU (left, right).

Figure 1 . 10 -

 110 Figure 1.10 -Proteus energy matrix calculation as described in the text. Matrix elements contain the interaction energy between couples of positions i and j when occupy two rotamers Rot i and Rot j of two types A and B

Figure 1 . 12 -

 112 Figure 1.12 -Bound and unbound states for two ligands L and L ′

44 ) 5 .Figure 1 . 14 -

 445114 Figure 1.14 -Titration of two ligands by variation of their relative concentrations. The binding free energy is calculated extrapolating the midpoint of the sigmoidal function (dashed line).

  THE JOURNAL OF CHEMICAL PHYSICS 149, 072302 (2018) Adaptive landscape flattening in amino acid sequence space for the computational design of protein:peptide binding Francesco Villa, Nicolas Panel, Xingyu Chen, and Thomas Simonson a) Laboratoire de Biochimie (CNRS UMR7654), Ecole Polytechnique, Palaiseau, France (Received 12 January 2018; accepted 16 March 2018; published online 4 May 2018)

FIG. 2 .

 2 FIG.2. Relative binding free energies ∆∆G from adaptive Monte Carlo (dashed line), experiment (black dots), and MD free energy simulations (gray dots). The MC values correspond to 50 000 sampled sequences, numbered by increasing ∆∆G values. For two sequences, only an experimental lower bound is known (interval indicated by vertical dashed lines). Labels indicate the mutations that define the experimental sequences (relative to the wildtype Sdc1).

  FIG.3. Left: Sequence logos from Monte Carlo and from an experimental library of peptides that bind Tiam1.53 Each column corresponds to a peptide position; position P0 is the C-terminus; the last five positions are shown. Amino acid types are shown with heights proportional to their abundancy and colored by physico-chemical categories. Right: The 30 strongest-binding MC sequences (ClustalW colors). The consensus sequence is shown at the bottom.

FFigure 2 . 1 -

 21 Figure 2.1 -Free energy surface (pictorial representation) as function of sequence s. Left: unbiased landscape, where the GMEC sits on the bottom and sequences A and B are separated by a free energy difference of several kcal/mol. Right: Biased simulation, with an almost flatten landscape. All sequences have similar probabilities and small free energy difference. Supposing perfect bias convergence, all sequences have the same probability (in the flattened landscape) and all their free energy differences are equal to zero.

Figure 3 . 2 -

 32 Figure 3.2 -Left: atom-drude + lone pair system. Right: SWM4-NDP Drude polarizable water model

3 .

 3 Point charge models other atoms in the residue. In the following example, the polarizability of the first oxygen depends on the positions of atoms C, CL and N. ANISOTROPY O C CL N A11 0.697 A22 1.219

Figure 3 . 3 -

 33 Figure 3.3 -Atom-Drude distance during MD simulation is shown in the two plots. Left: simulation where the Drude oscillator follows a normal regime, where the distance between

Figure 3 . 4 -

 34 Figure 3.4 -Interpolation of derivatives in thermodynamic integration. Calculated derivatives are represented with circles, the interpolated curve using a dashed line.

  [2002] ; Maragliano & Vanden-Eijnden [2006] ; Abrams & Tuckerman [2008] ; Cuendet & Tuckerman [2014]

65 )Figure 3 . 5 -

 6535 Figure 3.5 -Absolute binding free energy (∆G bind ) thermodynamic cycle. Top, left: bound state in solution (RL) s . Top, right: unbound state in solution, R s + L s . Bottom, left: restrained-bound state R s : L g with ligand in the gas phase. Bottom, right: unrestrainedunbound state R s + L g with the ligand in the gas phase.

+∞ 0 4πr 2 Figure 3 . 6 -

 0236 Figure 3.6 -Relative binding free energy ∆∆G(L → L ′ ) thermodynamic cycle. The transformations endpoints represent bound and unbound states with two different ligands L and L ′ .

Chapter 4 .

 4 PDZ-peptide binding specificity with polarizable free energy simulations molecules interact with hydrogen bond acceptors [fig. 4.2].

Figure 4 . 1 -

 41 Figure 4.1 -Alchemical transformation of two side chains A and B with Drude. In the figure, se use the red dashed line to indicate electrostatic interactions which are decoupled/coupled at λ elec =0.5. Scaling of LJ interactions is represented with blue dashed lines: side chain A is canceled between λ LJ = 0 → 1, while side chain B is inserted in the same interval.

Figure 4 . 3 -

 43 Figure 4.3 -Glutamate/Lysine Dual Topology (residue HEK)

Figure 4 . 4 -

 44 Figure 4.4 -Bond hierarchy in the Drude force field for a diatomic o triatomic molecule. Atoms are represented in yellow, while their Drude particles are represented in blue. Interactions are indicated with arrows.

Figure 4 . 5 -

 45 Figure 4.5 -Model residue used for the Drude dual topology definition. C β atoms are depolarized, in the sense that its Drude particle is completely removed from the system. Polarization and Thole factors are set to zero, and the oscillator charge is condensed in the reference atom.

Figure 4. 6 -

 6 Figure 4.6 -Thermodynamic cycle for the alchemical relative binding free energy of two side chains A and B. The transformations ∆G λ are carried out with the dual topology, where CB atoms are depolarized. The other legs in the thermodynamic cycle ∆G µ are the transformations used to compute the additional contributions µ(0 → 1), as described in the text.

Figure 4 . 7 -

 47 Figure 4.7 -Grounded simulation boxes and the more realistic macroscopic clusters, which are not grounded as described in the text. The two boxes used to simulate bound and unbound states are neutrilized by a different grounding potential. For the unbound state, the potential shift is close to Φ lv .

Figure 4 . 8 -

 48 Figure 4.8 -Closure errors for relative binding free energies calculated with Polarizable or additive force fields

Figure 4 . 9 -

 49 Figure 4.9 -Structural views based on the Amber ff99SB MD simulations for the wild-type/Sdc1 complex (top) and its E4K (middle) and E4L (bottom) mutants. The 3D structures are the mean structures from the MD.

Figure 4 . 10 -

 410 Figure 4.10 -Structural views based on the Drude polarizable MD simulations for the wild-type/Sdc1 complex (top) and its E4K (middle) and E4L (bottom) mutants. The 3D structures are the mean structures from the MD.

  [2013]). Sophisticated methods have been proposed to accelerate conformational sampling. They include adaptive biasing of the energy surface with metadynamics(Laio & Parrinello [2002]), λ-dynamics[START_REF] Hayes | Adaptive landscape flattening accelerates sampling of alchemical space in multisite lambda dynamics[END_REF]), or the orthogonal space random walk[START_REF] Zheng | Random walk in orthogonal space to achieve efficient free-energy simulation of complex systems[END_REF]).

1

 1 Phosphotyrosine in Tiam1 bindingPhosphorylation is among the most diffused post-translational modifications that regulate PPIs in cell signalling networks[START_REF] Khoury | Proteome-wide post-translational modification statistics: frequency analysis and curation of the swiss-prot database[END_REF] ;[START_REF] Pawson | Protein phosphorylation in signaling -50 years and counting[END_REF]

Figure 5 . 1 -

 51 Figure 5.1 -Tiam1-Sdc1 and Tiam1-pSdc1 structures superposition. The wildtype complex is shown in green, the phosphorylated complex is shown in cyan. From crystal structures 4GVD and 4GVD.

Figure 5 . 2 -

 52 Figure 5.2 -Methyl phosphate, hydrogen phosphate and phosphotyrosine structures.

2

  Methods: overview Parametrization was done following the established Drude parametrization strategy described in Anisimov et al. [2005] ; Lemkul et al. [2016]. The empirical, molecular-mechanics (MM) Drude energy function was used to fit target quantum-mechanical (QM) quantities [fig. 5.3].

Figure 5 . 3 -

 53 Figure 5.3 -Parametrization procedure overview.

Figure 5 . 4 -

 54 Figure 5.4 -Methyl phosphate and hydrogen phosphate 2D structures. Left: MP -. Center: MP 2-. Right: P 2- i

Figure 5 . 5 -

 55 Figure 5.5 -Internal and Outer shell binding for methyl phosphate.

Figure 5 . 6 -

 56 Figure 5.6 -QM electrostatic potential maps for MP 2-. Left: Grid for evaluation of electrostatic potential. Right: perturbing charges, in green.

Figure 5 . 7 -

 57 Figure 5.7 -Water-solute scans. Left: water pose built using a dummy particle (shown in blue) on the MP 2-acceptor O2. Right: QM interaction energy scan for the MP 2-acceptor O2.

Figure 5 . 8 -

 58 Figure 5.8 -Thermodynamic cycle used to compute absolute binding free energies for Mg 2+ and methyl phosphate complexes (vertical arrow on the left, indicated using the label Bind).The reaction is decomposed in several steps labelled with integers from 1 to 6 and described in the text.

Figure 5 . 9 -

 59 Figure5.9 -Water scan for DMP. Left: interaction energy as function of the water-acceptor O11 distance. Several lines represent the scan performed using different parameters set. In blue, we show the scan performed using the official Drude parameters for DMPN (nucleic acids topology). In red, we show the scan performed after substitution of the atom type OD30BN with the lipid atom type OD30B. Charges were reoptimized, giving perfect agreement showed in the black line. QM is shown in grey. Right: geometry of the complex used for the scan.

Figure 5 . 14 -

 514 Figure 5.14 -Conformational energies. Left: QM vs MM comparision. Right: superposition of the 200 random configurations generated by random linear combination of the QM normal modes. For both MP -and MP 2-QM/MM energies have correlation coefficient of 0.9. Root mean square deviations are 0.3 kcal/mol for MP -and 0.6 kcal/mol for MP 2-.

Figure 5 . 16 -

 516 Figure5.16 -Magnesium-phosphate binding-unbinding Black line : distance between Mg 2+ and P atoms, extracted from a 10ns simulation of the MP -:Mg 2+ complex without restrain potential. Red ,dashed line: limit for OS bound state where the magnesiumphosphorus distance is between 4.1 and 5.9Å

Figure 5 . 17 -

 517 Figure 5.17 -Drude topology files for p-Cresol (left) and tyrosine (right). Equivalent atoms are shown on the same line.

Figure 5 .

 5 Figure 5.18p-Cresol 3D structure

figure 5 .

 5 figure 5.19. We created the initial model combining atom types and electrostatic parameters of the existing p-Cresol and of our MP 2-model.

Figure 5 . 19 -

 519 Figure 5.19 -Phosphorylated p-CRESOL (CREP) from non-phosphorylated p-CRESOL (CRES). The phosphorylated form can be seen as a patch applied to the non-phosphorylated compound. Atoms on the right of the dashed line are modified, while atoms on the left are kept as in CRES. Groups of equivalent atoms are (HB1-HB2-HB3), (CD1,CD2), (CE1-CE2), (HD1-HD2), (HE1-HE2) and (O2-O3-O4).

Figure 5 . 20 -

 520 Figure 5.20 -Phosphorylated p-Cresol Water scans

Figure 5 . 21 - 2 Figure 5 . 22 -

 5212522 Figure 5.21 -CREP 2--water interaction Grey: QM; Black: MM

Figure 5 . 23 -

 523 Figure 5.23 -pTyr 2--K879 interactions during 20 ns of MD. Left: crystal structure (light blue) and average structure from the simulation (green). Right: phosphorus-nitrogen distance during MD. Frames were extracted each 40 ps. The straight line indicates the distance in the crystal structure (3.61 Å).

  Les interactions protéine-protéine (IPPs) médient la signalisation cellulaire. Leur ingénierie peut fournir des informations et conduire au développement de molécules thérapeutiques. Les domaines PDZ sont des médiateurs importants de IPPs. Elles lient les 4--10 résidus Cterminaux de protéines cibles. Elles lient aussi les peptides correspondants, qui ppeuvent servir de systèmes modèle ou d'inhibiteurs. Nous avons développé deux approches computationnelles et les avons appliquées au domaine PDZ de la protéine Tiam1, un facteur d'échange pour la protéine Rac, impliqué dans la protrusion neuronale. Sa cible est la protéine Syndecan1. Des affinités expérimentales sont connues pour le peptide Cterminal, noté Sdc1, et plusieurs mutants; elles ont servi pour tester les calculs. Nous avons d'abord développé une méthode de dessin computationnel haut débit. Une simulation Monte Carlo est faite où les chaines latérales de la protéine et du peptide peuvent changer de conformères et certaines positions peuvent muter. Le solvant est implicite. Le paysage énergétique est applati par la méthode adaptative de Wang-Landau, de sorte qu'un vaste ensemble de variants est échantillonné. Effectuant des simulations distinctes du complexe et du peptide seul nous avons obtenu les énergies libres relatives d'association de 75,000 variants en heure CPU sur une machine de bureau. Les valeurs sont compatibles avec les quelques données expérimentales disponibles. Ensuite, nous avons développé une approche beaucoup plus détaillée et réaliste. Soluté et solvant sont décrits par un champ de force atomique, qui représente explicitement la polarisation électronique: le champ de force Drude de Charmm. La polarisabiité peut être importante car les résidus de l'interface PDZ:peptide passent, lors de l'association, d'un environement riche en solvant à un autre pauvre en solvant. Nous avons fait des simulations alchimiques d'énergie libre pour comparer quatre variants du peptide qui diffèrent par une ou deux chaines latérales ioniques. Les résultats sont en bon accord avec l'expérience. Les chamops de force additifs Charmm et Amber, qui représentent la polarisabilité implicitement, donnent un moins bon accord. Ces calculs sont le premier exemple de simulations alchimiques d'énergies libre d'association relatives protéine:ligand avec un champ de force polarisable. Enfin, pour une modélisation future de peptides phophorylés, nous avons étendu le champ de force Drude pour inclure le méthyl phosphate et la phospho tyrosine. Il en résulte un excellent accord vaec les affinités expérimentales phosphate:magnésium.

TABLE I .

 I Mutation energies for the unbound peptide (kcal/mol). Values for each amino acid type and peptide position, with Ala as the reference, averaged over the unbound peptide simulation. "Mean" indicates an average over positions ☞4, . . ., ☞1.

				∆E
	Type	Mean	P ☞4	P ☞3	P ☞2	P ☞1
	ALA	0.0	0.0	0.0	0.0	0.0
	ARG	☞46.9	☞47.3	☞48.1	☞48.7	☞48.6
	ASN	☞15.8	☞16.1	☞16.9	☞16.2	☞16.3
	ASP	☞18.2	☞18.5	☞19.1	☞17.5	☞17.8
	CYS	☞0.4	0.1	☞0.7	☞0.4	☞0.5
	GLN	☞12.6	☞12.5	☞14.1	☞12.7	☞12.8
	GLU	☞15.7	☞15.1	☞16.3	☞14.3	☞14.7
	HIS	14.8	15.2	15.4	14.5	14.1
	ILE	2.9	1.8	1.3	2.2	2.2
	LEU	☞3.9	☞4.2	☞4.8	☞4.7	☞4.7
	LYS	☞5.5	☞6.0	☞6.7	☞7.3	☞7.1
	MET	☞1.9	☞2.2	☞2.9	☞2.4	☞2.4
	PHE	☞1.1	☞1.3	☞2.2	☞0.9	☞1.2
	SER	☞4.4	☞4.4	☞4.7	☞4.3	☞4.5
	THR	☞7.1	☞6.8	☞7.9	☞7.1	☞7.4
	TRP	☞0.7	☞2.2	☞1.9	☞1.2	☞1.6
	TYR	☞6.8	☞7.1	☞8.4	☞6.7	☞7.2
	VAL	☞3.1	☞3.8	☞4.5	☞3.5	☞3.9

TABLE II .

 II Experimental and computed binding free energies for PDZ:peptide complexes (kcal/mol), with wildtype Sdc1 (WT) as the reference.

	Peptide	Mutations	Expt.	MC	|Error|
	TKQEEFYA	WT	0.0	0.0	. . .
	TKQEEFTA	Y ☞1 T	0.7 a	☞0.3	1.0
	TKQEDFYA	E ☞3 D	1.8 a	☞0.3	2.1
	TKQEDFTA	E ☞3 DY ☞1 T	0.9	0.1	0.8
	TKQETFKA	E ☞3 TY ☞1 K	1.3	0.5	0.8
	TKQLEFYA	E ☞4 L	0.6	0.6	0.0
	TKQKEFYA	E ☞4 K	0.8	1.8	1.0
	TKQEENYA	F ☞2 N	>1.3 b	2.6	0.0
	TKQEEEYA	F ☞2 E	>1.3 b	0.8	0.5
	TKQEEIYA	F ☞2 I	0.8	7.7	6.9

a From MD free energy simulations.

66

b Only an experimental lower bound is available.

54

listed in Table

II

. Comparing the MC to the reference values, there are just two large errors. For the E ☞3 D peptide single mutant, the MC value is ☞0.3 kcal/mol, compared to an MD

Table 2 .

 2

	PB/LIE MC

1 -PB/LIE binding free energies for 14 selected peptide sequences. Each sequence specifies the amino acids at positions P -4,-3,-2,-1,0 . Specific contributions to the empirical ∆G are given on different columns. All terms are in kcal/mol, except values in the SA column (the surface area) which are in Å 2 . PB, vdW, and SA contributions were weighted by optimized parameters to compute ∆G [eqn. 2.1]. Optimized weights were α=0.02 (vdw), β=0.25 (elec) and γ=-4 cal/mol/Å 2 (Panel et al.

[2017]

). Relative binding free energies of each sequence were obtained subtracting the ∆G of the WT sequence EEFYA.

Table 3 . 1 -

 31 Drude-C36 comparision for MD run using NAMD 2.12. System details are given in the text.

		C36	Drude
	Number of atoms	31 994	53 393
	Timestep	2 fs	1 fs
	Number of cores 224 (8x28) 224 (8x28)
	Performance	75 ns/day 18 ns/day

The

  Chapter 4. PDZ-peptide binding specificity with polarizable free energy simulations

	RESI HEK 0.00
	! Backbone
	GROUP
	ATOM N	ND2A2	-0.382 ALPHA -1.942 THOLE 0.250
	ATOM HN	HDP1A	0.272
	ATOM CA	CD31C	0.169 ALPHA -0.960 THOLE 1.078
	ATOM HA	HDA1A	-0.017
	ATOM C	CD2O1A	0.497 ALPHA -0.675 THOLE 0.295
	ATOM O	OD2C1A	0.000 ALPHA -0.651 THOLE 0.310
	ATOM LPOA LPDO1	-0.312
	ATOM LPOB LPDO1	-0.227
	! Glutamate sidechain
				!excluded interactions:
	ATOM CBE CD32B	-0.066	cbk hb1k hb2k cgk [...]
	ATOM HB1E HDA2A	0.033	cbk hb1k hb2k cgk [...]
	ATOM HB2E HDA2A	0.033	cbk hb1k hb2k cgk [...]
	ATOM CGE CD32A	-0.190 ALPHA -2.528 THOLE 1.414	cbk hb1k hb2k cgk [...]
	ATOM HG1E HDA2A	0.004	cbk hb1k hb2k cgk [...]
	ATOM HG2E HDA2A	0.004	cbk hb1k hb2k cgk [...]
	ATOM CDE CD2O2A	0.708 ALPHA -1.016 THOLE 0.899	cbk hb1k hb2k cgk [...]
	ATOM OE1 OD2C2A	0.003 ALPHA -0.699 THOLE 2.399	cbk hb1k hb2k cgk [...]
	ATOM OE2 OD2C2A	0.003 ALPHA -0.699 THOLE 2.399	cbk hb1k hb2k cgk [...]
	ATOM LP1A LPD	-0.383	cbk hb1k hb2k cgk [...]
	ATOM LP1B LPD	-0.383	cbk hb1k hb2k cgk [...]
	ATOM LP2A LPD	-0.383	cbk hb1k hb2k cgk [...]
	ATOM LP2B LPD	-0.383	cbk hb1k hb2k cgk [...]
	! Lysine sidechain
				!excluded interactions:
	ATOM CBK CD32A	-0.066	cbe hb1e hb2e cge [...]
	ATOM HB1K HDA2A	0.033	cbe hb1e hb2e cge [...]
	ATOM HB2K HDA2A	0.033	cbe hb1e hb2e cge [...]
	ATOM CGK CD32A	-0.156 ALPHA -1.660 THOLE 1.300	cbe hb1e hb2e cge [...]
	ATOM HG1K HDA2A	0.078	cbe hb1e hb2e cge [...]
	ATOM HG2K HDA2A	0.078	cbe hb1e hb2e cge [...]
	ATOM CDK CD32A	-0.156 ALPHA -1.660 THOLE 1.300	cbe hb1e hb2e cge [...]
	ATOM HD1K HDA2A	0.078	cbe hb1e hb2e cge [...]
	ATOM HD2K HDA2A	0.078	cbe hb1e hb2e cge [...]
	ATOM CE	CD32A	0.043 ALPHA -1.656 THOLE 0.895	cbe hb1e hb2e cge [...]
	ATOM HE1 HDA2C	0.143	cbe hb1e hb2e cge [...]
	ATOM HE2 HDA2C	0.143	cbe hb1e hb2e cge [...]
	ATOM NZ	ND3P3A -0.349 ALPHA -1.298 THOLE 0.895	cbe hb1e hb2e cge [...]
	ATOM HZ1 HDP1B	0.340	cbe hb1e hb2e cge [...]
	ATOM HZ2 HDP1B	0.340	cbe hb1e hb2e cge [...]
	ATOM HZ3 HDP1B	0.340	cbe hb1e hb2e cge [...]

Table 5 . 1 -

 51 Electrostatic parameters optimization workflow for MP -, with intermediate and final parameter sets. Each set was optimized fitting different QM target data, as described in the text.6.2 Fitting the MP -potential maps

	Atom Parameters	Atom Parameters
	O1	polarizability, Thole	O3	charge, polarizability, Thole
	P1	charge, polarizability, Thole	O4	-redundant -
	O2	charge, polarizability, Thole	C1	-fixed -
	LPAO2 charge	H11 -fixed -
	LPBO2 charge	H11 -redundant -
	H2	charge	H11 -redundant -
	Table

table 5

 5

	.3 and table 5.4.

Table 5 . 3 -

 53 Restrained fit quality for MP -including all contributions to the target function from equations 5.5, 5.6 and 5.8. The target function is close to the ESP fit quality, meaning that optimized parameters don't deviate too much respect to applied restraints.

	Atom charge	α	τ	Atom Charge	α	τ
	O1	-0.479 -0.360 1.305	O3	-0.738 -0.497 1.391
	P1	1.122 -1.254 1.395	O4	-0.738 -0.497 1.391
	O2	0.097 -0.457 1.194	C1	0.072 -1.642 0.767
	LPAO2 -0.290			H11	0.026
	LPBO2 -0.290			H12	0.026
	H2	0.166			H13	0.026

Table 5 .

 5 4 -Intermediate parameters set 1 for MP -, "post-ESP" fit.

Table 5 . 5 -

 55 Molecular polarizabilities in (Å 3 ) fitted using initial values (intermediate parameters set 1, "post-ESP"fit) and then optimized with the intermediate parameters set 2.

	.5 and optmized parameters (set2) in table

Table 5 . 6 -

 56 Intermediate parameters set 2 for MP -, after fitting the molecular polarizability.

		Charge	α	τ	Atom Charge	α	τ
	O1	-0.479 -0.901 0.831	O3	-0.738 -0.931 1.083
	P1	1.122 -0.974 2.098	O4	-0.738 -0.931 1.083
	O2	0.097 -0.927 1.300	C1	0.072 -1.642 0.862
	LPAO2 -0.290			H11	0.026
	LPBO2 -0.290			H12	0.026
	H2	0.166			H13	0.026

Table 5 . 7 -

 57 Solute-water interaction energies for MP -, fitted using the intermediate set (after fitting the molecular polarizability) and then optimized with the intermediate parameters set 3. We highlight the strong deviation of r min for acceptor O1, the methyl phosphate bridging oxygen.

			QM	MM, set 2	MM, set 3
	atom r min	e min	r min	e min	r min	e min
	O1	1.9 -11.663 1.5 -12.675 1.5 -12.663
	O2	1.9 -10.583 1.9	-9.456	1.8 -10.892
	O3	1.8 -13.433 1.7 -13.030 1.7 -13.759
	O4	1.8 -12.934 1.7 -13.206 1.7 -14.076
	H	2.0	0.073	2.3	3.036	2.0 0.5518
			µ x	µ y	µ z	µ
	QM		-4.2395 2.5428 -1.0626 5.0565
	MM, set 2 -2.7443 1.6167 -0.5710 3.2358
	MM, set 3 -3.1737 2.5779 -0.5337 4.1234

Table 5 . 8 -

 58 Molecular dipole moment for MP -(Debyes), obtained using the intermediate set 2 parameters and then optimized with the parameters set 3.

Table 5 . 9 -

 59 Parameters set 3 for MP -, after fitting the water interaction energy scans. The intermediate set 2 is shown for comparison.

	6. Results

Table 5 .

 5 2-, QM vs MM. 10 -MP -, MP 2-and P2-

	molecule method	α xx α yy α zz	α	method	µ x	µ y	µ z	µ
		QM	8.28 7.26 6.92 7.49	QM	-4.24 2.54 -1.06 5.06
	MP -	QM×0.85 7.04 6.17 5.88 6.36	MM	-2.74 1.62 -0.57 3.23
		MM	7.09 5.96 5.88 6.31
		QM	9.83 7.80 7.92 7.49	QM	-4.57 -0.45 0.01 4.59
	MP 2-	QM×0.85 8.36 6.63 6.73 7.24	MM	-4.23 0.47 0.00 4.26
		MM	7.59 6.50 6.49 6.86
		QM	7.05 6.72 6.69 6.82	QM	0.12 -1.10 -2.36 2.61
	P 2-i	QM×0.85 5.99 5.71 5.68 5.80 MM 4.75 4.80 4.47 4.67	MM	0.08 -0.67 -1.10 1.29
			i	molecular polarizabilities (Å 3 ) and dipole moments (D)
							QM		MM
			molecule atom r min e min r min e min
				O1	1.9 -11.66 1.8 -11.40
				O2	1.9 -10.58 1.8 -10.93
			MP -	O3	1.8 -13.43 1.7 -13.66
				O4	1.8 -12.93 1.7 -13.93
				H	2.0	0.07	2.0	0.55
				O1	1.7 -21.08 1.7 -21.18
			MP 2-

O2 1.6 -22.75 1.6 -22.70 O3 1.6 -22.88 1.6 -22.80 O1 1.7 -20.77 1.6 -21.00 P 2- i O2 1.6 -22.41 1.6 -23.15 O3 1.6 -22.61 1.6 -22.62

Table 5 .

 5

	6. Results

11 -MP -, MP 2-and P 2- i water interaction distances (Å) and energies (kcal/mol) Figure 5.10 -Solute-water energy scan geometry for MP -.

Table 5 .

 5 2-QM/MM energies have correlation coefficient of 0.9. Root mean square deviations are 0.3 kcal/mol for MP -and 0.6 kcal/mol for MP 2-. different set of parameters. First, we used all the special NBFIX and NBTHOLE parameters from the existing force field, listed in table 5.12. Then we have run a second round. Instead of using the official NBFIX/NBTHOLE for the bridging oxygen, we applied LJ parameters of the Drude lipid atom type (OD30B) as we did to correct the water scans. Results are given in table 5.13. 12 -Pair-specific parameters for Mg 2+ and Na + in Drude, to balance their interactions with phosphates and water. Values are extracted fromLemkul & MacKerell [2016a] ;

					type		NBFIX	NBTHOLE
							ǫ(kcal/mol) r(Å)
		Drude charge of water DOH2		-0.07500 2.65500	-
	Mg 2+ bridging oxygen	OD30B(N) -0.40592 4.40157	1.151
		nonbridging oxygen	OD2C2C	-0.31211 2.59990	1.449
	Na +	nonbridging oxygen	OD2C2C	-0.04696 3.49668	1.449
		hydroxyl oxygen	OD31A		-0.06000 3.07000	1.820
	Savelyev & MacKerell [2014]
					QM	NBFIX	OD30B
	molecule atom type	r min	e min	r min	e min	r min	e min
	DMPN	O1	OD30BN 1.8 -278.27 3.35 -152.20 1.9 -265.45
			O1-inv OD30BN 3.0 -313.31 3.2 -254.16 2.45 -303.32
			O3	OD2C2C 1.8 -299.47 1.9 -261.07 1.9 -260.87
	MP -	O1	OD30B	1.9 -282.93 3.25 -168.71 1.9 -297.55
			O2	OD31A	1.9 -258.80 2.2 -258.40 1.95 -262.60
			O3	OD2C2C 1.8 -290.17 2.15 -272.85 1.95 -272.75
			O4	OD2C2C 1.8 -298.76 2.16 -295.15 1.95 -295.10
	MP 2-	O1	OD30B	1.8 -480.05 3.1 -315.90 1.9 -519.62
			O2	OD2C2C 1.7 -483.41 1.9 -508.07 1.9 -469.15
			O3	OD2C2C 2.0 -485.08 1.9 -502.42 1.9 -463.35

Table 5 .

 5 13 -Ions scans with Mg 2+ Distances are in Å and interaction energies in kcal/mol.

Table 5 .

 5 . Results are shown in table 5.14. We obtained good QM/Drude agreement, for both distances and interaction energies. 14 -Ions scans with Na + Distances are in Å and interaction energies in kcal/mol.

	QM	NBFIX

Table 5 .

 5 15 -Mg 2+ -phosphate binding free energies from simulations and experiment In kcal/mol. a At the standard state concentration. b Obtained by subtracting the Mg 2+ value from the value in "total" column. c Values are for H 2 PO -

		recharge Mg 2+
	solute or complex	size (ρ M D ) correction ∆G elec	restore LJ a apply ∆G LJ restraint total b ∆G bind Expt. c
	Mg 2+ :MP -(OS) -415.6 Mg 2+ :MP -(IS) -414.4	0.0 0.0	+0.8 +1.0	+0.4 +1.1	-414.4 -412.3	-2.1 0.0	-1.7
	Mg 2+ :MP 2-(OS) -417.5 Mg 2+ :MP 2-(IS) -416.5	+1.0 +1.0	+1.0 +1.6	+0.4 +1.1	-415.1 -412.8	-2.8 -0.5	-3.7
	Mg 2+ :P 2-i (OS) Mg 2+	-418.5 -414.1	+1.0 +0.9	+0.8 +0.9	+0.4 -	-416.3 -412.3	-4.0 -	-3.7 -
					4 ≡ P -i and HPO 2-4 ≡ P 2-i	Verbeeck
	et al. [1984] ; Alberty & Goldberg [1992].

  The final models for methyl phosphate were included in the official Drude polarizable force field. We give now topologies and parameters, as they are in the new force field release

						6. Results
	RESI MP_1 BONDS	-1.00	! CH4O4P Methylphosphate, anionic
	GROUP OD30D PD1AN 240.00 1.600	! Villa et al., 2018
	ATOM O1 CD33C OD30D 335.00 1.440 OD30D -0.546 ALPHA -0.901 THOLE 0.811	!	H11
	ATOM P1 OD31A PD1AN 237.00 1.610 PD1AN 1.122 ALPHA -0.974 THOLE 2.098	!	|
	ATOM O2 OD2C2C PD1AN 525.00 1.500 OD31D 0.000 ALPHA -0.927 THOLE 1.100	! H13--C1--H12
	ATOM LPA HDP1A OD31A 536.50 0.970 LPDNA1 -0.300	!	|
	ATOM LPB CD33C HDA3A 322.00 1.111 LPDNA1 -0.300	!	O1
	ATOM H2	HDP1A	0.360	!	|
	ATOM O3 ATOM O4 ANGLES ATOM C1 CD33C OD30D PD1AN OD2C2C OD2C2C CD33C ATOM H11 HDA3A OD30D PD1AN OD31A ATOM H12 HDA3A OD30D PD1AN OD2C2C 98.90 113.00 -0.778 ALPHA -0.921 THOLE 1.083 -0.778 ALPHA -0.921 THOLE 1.083 -0.008 ALPHA -1.642 THOLE 0.862 40.00 120.50 0.076 48.10 108.00 0.076 ATOM H13 HDA3A HDA3A CD33C OD30 60.00 109.50 0.076 BOND O1 HDP1A OD31A PD1AN 30.00 115.00 40.00 2.3500 ! Urey-Bradley term ! O4==P1==O3 (-) ! | ! O2 ! \ ! H2 P1 BOND O1 OD31A PD1AN OD2C2C 98.90 111.00 C1 BOND O2 OD2C2C PD1AN OD2C2C 120.00 125.00 P1 BOND O3 P1 HDA3A CD33C HDA3A 35.50 108.40
	BOND O4 BOND O2 DIHEDRALS BOND C1 CD33C OD30D PD1AN P1 H2 H11 BOND C1 H12 CD33C OD30D PD1AN	OD2C2C 0.100 3 OD31A 0.500 3	0 .0 ! from nucleic acids 0.0
	6.10 Final MP and P 2-i BOND C1 H13 CD33C OD30D PD1AN OD31A BOND O2 LPA PD1AN OD30D CD33C HDA3A	0.950 2 topology and parameters 0.0 0.000 3 0.0
	We give first the final model for hydrogen phosphate P 2-i =HPO 2-BOND O2 LPB OD30D PD1AN OD31A HDP1A 0.2000 1 180.0 4 , for which vdW optimiza-LONEPAIR relative LPA O2 P1 H2 distance 0.35 angle 110.00 dihe 90.00 OD30D PD1AN OD31A HDP1A 1.0000 2 0.0 tion should be improved LONEPAIR relative LPB O2 P1 H2 distance 0.35 angle 110.00 dihe 270.00 HDP1A OD31A PD1AN OD2C2C 0.2000 3 0.0
	ANISOTROPY O2 P1 LPA LPB A11 0.76473 A22 1.16239 CD33C OD30D PD1AN OD2C2C 0.100 3 0.0
	RESI HP_2 -2.000 ! HPO4 Hydrogenphosphate, dianionic ! PATCH FIRST NONE LAST NONE	H
	GROUP NONBONDED E14FAC 1.000000		!	|
	ATOM P1 RESI MP_2 PD1AN OD30D 0.0 -0.1700 1.020 ALPHA -1.244 -2.00 ! CH3O4P Methylphosphate, dianionic THOLE 1.478 ! 1.7700 ! from DMP @ lipids: new type for MP O1
	ATOM O2 GROUP PD1AN 0.0000 -0.3200 OD2C2C -0.860 ALPHA -0.921 ! Villa et al., 2018 THOLE 0.803 ! 1.9000 ATOM O3 OD2C2C -0.860 ALPHA -0.921 THOLE 0.803 ! (-)O4--P1--O3(-) | ATOM O4 OD2C2C -0.860 ALPHA -0.921 THOLE 0.803 ! | ATOM O1 OD31D 0.000 ALPHA -0.927 THOLE 0.900 ! O2 ATOM LPA LPDNA1 -0.320 ATOM LPB LPDNA1 -0.320 ATOM H HDP1A 0.200 ! H11 ATOM O1 OD30D -0.685 ALPHA -1.291 THOLE 1.091 ! CD33C 0.0000 -0.0780 1.9400 | ATOM P1 PD1AN 1.050 ALPHA -1.244 THOLE 1.678 ! OD31A 0.0000 -0.1500 1.7650 H13--C1--H12 ATOM O2 OD2C2C -0.860 ALPHA -0.951 THOLE 1.083 ! OD2C2C 0.0000 -0.0700 1.8650 | ATOM O3 OD2C2C -0.860 ALPHA -0.951 THOLE 1.083 ! HDP1A 0.0000 -0.0100 0.4000 O1 ATOM O4 OD2C2C -0.860 ALPHA -0.951 THOLE 1.083 ! HDA3A 0.0000 -0.0240 1.3400 | ATOM C1 CD33C 0.020 ALPHA -1.642 THOLE 0.862 MAGD 0.0000 -0.0500 1.1264156 ! (-)O4==P1==O3(-) LONEPAIR relative LPA O1 P1 H distance 0.35 angle 110.00 dihe 90.00 LONEPAIR relative LPB O1 P1 H distance 0.35 angle 110.00 dihe 270.00 BOND O1 P1 BOND O2 P1 BOND O3 P1 BOND O4 P1 BOND O1 H ATOM H11 HDA3A 0.065 ! || ATOM H12 HDA3A 0.065 ! NBFIX O2 ATOM H13 HDA3A MAGD DOH2 -0.07500 2.65500 ! Lemkul & MacKerell, J.Phys.Chem.B 2016 0.065 MAGD OD2C2C -0.31211 2.59990 ! Lemkul & MacKerell, J.Phys.Chem.B 2016 BOND O1 P1 MAGD LPD -0.16504 2.20163 ! Lemkul & MacKerell, J.Phys.Chem.B 2016 BOND O1 C1 SODD OD2C2C -0.04696 3.49668 ! nucleic acids, Savelyev-2014 BOND O2 P1 SODD OD30BN -0.02510 3.38168 ! nucleic acids, Savelyev-2014 BOND O3 P1 BOND O4 P1 SODD LPD -0.09000 2.92000 ! nucleic acid bases BOND O1 LPA BOND C1 H11 SODD OD31A -0.06000 3.07000 ! ethanol, ser, thr, Hui-2015 BOND O1 LPB BOND C1 H12 MAGD OD30D -0.40592 2.25157 ! new, Villa et al, 2018 PATCH FIRST NONE LAST NONE BOND C1 H13 SODD OD30D -0.02510 3.38168 ! transferred from DMP OD30BN
	PATCH FIRST NONE LAST NONE
	NBTHOLE
	MAGD	OD2C2C 1.44900 ! Lemkul & MacKerell, J.Phys.Chem.B 2016
	MAGD	OD30D	1.15100 ! Lemkul & MacKerell, J.Phys.Chem.B 2016 (OD30BN) )
	MAGD	ODW	1.51567 ! Lemkul & MacKerell
	SODD	OD31A	1.82000 ! ethanol, ser, thr, Hui-2015
	MAGD	OD31A	1.95000 ! new, Villa et al, 2018

Table 5 .

 5 16 -CREP-water interaction distances (Å) and energies (kcal/mol) 129 Chapter 5. Classical Drude Model for methyl phosphate and phosphotyrosine

			18.02 1.8 -17.93
	O3	1.7 -19.50 1.8 -19.48
	P	3.3 -24.18 3.2 -24.15
	HE1 2.8	3.74	2.6	3.50

′i → t i . The total energy change for the mutation

Remerciements

Je tiens à remercier les membres du jury Amedeo Caflisch, Jean-Philip Piquemal, Nathalie ACKNOWLEDGMENTS Discussions with David Mignon (Ecole Polytechnique) and Ernesto J. Fuentes (U. of Iowa) are gratefully acknowledged.

correction is obtained by computing a charging free energy for a periodic and a non-periodic system and then taking the difference: ∆G(L → ∞) = ∆G(L) + ∆G P E (L → ∞) -∆G P E (L) (3.72) straint potential was, for both OS and IS states: where r is the distance between Mg 2+ and the MP phosphorus atom P. For IS, [a; b] was set to [0.0 : 3.5Å] and c=4 kcal/mol/Å 2 . For OS we used [4.1 : 5.9Å] and c=8 kcal/mol/Å 2 . The restraint free energy was obtained analitically, using

where [4] We alchemically reintroduced interactions between Mg and MP [START_REF] Pokala | Energy functions for protein design I: Efficient and accurate continuum electrostatics and solvation[END_REF] We restored the standard-state concentration, with the corresponding free-energy contribution KT ln(ρ ∞ /ρ 0 ).

[6] We removed the restraint and the corresponding free energy cost was negligible.

At the end, all terms involving ρ ∞ (steps 1,3,4) cancel out, leaving a logarithmic dependency on the standard-state concentration kT lnρ 0 . Just free energies for steps 2 and 4 were computed with alchemical free energy differences. Details about these simulations are given below.

Simulations setup

System preparation was done with the CHARMM program. For unbound state simulations (point 2) the solute was one Mg 2+ ion placed at the center of the box. For the bound state

Results

Fitting QM interactions with Mg and Na ions

Pair-specific NBFIX and NBTHOLE terms were optimized by the Drude force field developers to tune Lennard-Jones and electrostatic interactions between magnesium and nucleic acids compounds in water (Lemkul & MacKerell [2016a]). Others specific terms for sodium were developed earlier by [START_REF] Savelyev | All-atom polarizable force field for DNA based on the classical drude oscillator model[END_REF]. Details about parameter optimization and QM target data can be found in the publications. We recall that NBFIX parameters define pair-specific LJ interactions for a couple of atom types i and j. Two variables indicating minimum energy ǫ ij min and minimum radius r ij min are used to compute the LJ interaction between atoms of type i and j at distance r. The use of NBFIX does not alter LJ interactions between atoms of type i or j and all other types k = i,j in the force field. NBTHOLE is analogous to NBFIX. A pair-specific screening parameter τ ij is used to define the smeared charge distribution [eqn. 3.22] for the calculation of electrostatic interaction between a couple of atoms of type i and j.

We give the list of these optimized parameters in table 5.12. A first NBFIX was applied between magnesium (atom type MAGD) and the Drude charge attacched to the oxygen atom of SWM4-NDP water (atom type DOH2). The authors modified this interaction after observing undesired behaviour of Mg 2+ -water complexes during MD simulations. More precisely, they noticed that after few nanoseconds of dynamics the standard Drude force field parameters produced frequent exchange of water molecules in complex with Mg 2+ , where one ion interacts with six SWM4-NDP molecules. This should not happen, because these exchanges should take place at the microsecond timescale [START_REF] Neely | Rate of water exchange from hydrated magnesium ion[END_REF]). The introduced NBFIX was then used to damp the too-strong dipole-dipole repulsion between water molecules, which is induced by the too-strong dipole reponse of water molecules when interact with Mg 2+ . The authors also observed that the new parameters improved the agreement between QM and MM Mg 2+ -water interaction energies. We conclude that this NBFIX should always be applied in our simulations. Two couples of NBFIX and NBTHOLE parameters were introduced to better reproduce QM geometries for several nucleic acid moieties (Ade, Gua, Cyt, Ura, DMP) in complex with Mg 2+ . The authors considered Inner-shell and Outer-shell coordinated complexes, parameters were optimized by a grid search. We note that the value of r min for the phosphate bridging oxygen (O1) is considerably big when compared with the one for the nonbridging oxygen.

We performed QM/Drude scans with magnesium and sodium on methyl phosphate and for the bridging oxygen of the existing DMP from the official nucleic acids Drude force field (DMPN). For DMPN, we scanned the bridging oxygen in two directions. First, imposing "direct interaction" with Mg 2+ , then in the opposite direction when magnesium is closer to the phosphate nonbridging oxygens [fig. Dans une seconde phase, nous avons développé une approche beaucoup plus détaillée et réaliste, avec un coût de calcul plus élevé. Cette approche est basée sur des simulations de dynamique moléculaire et implique l'utilisation de techniques de calcul d'énergie libre. Soluté et solvant sont décrits par un champ de force atomique, qui représente explicitement la polarisation électronique: le champ de force Drude de Charmm. La polarisabilité peut être importante car les résidus de l'interface PDZ:peptide passent, lors de l'association, d'un environement riche en solvant à un autre pauvre en solvant. Nous avons fait des simulations alchimiques d'énergie libre pour comparer quatre variants du peptide qui diffèrent par une ou deux chaines latérales ioniques. Pendant chaque transformation alchimique, les chaînes latérales d'intérêt sont mutées en utilisant un chemin d'états non physiques. Les transformations sont réalisés avec plusieurs simulations indépendantes, nécessitent l'utilisation de plusieurs processeurs en parallèle. Nos calculs sont principalement effectués sur le superordinateur OC-CIGEN du CINES. Les résultats sont en bon accord avec l'expérience. Les champs de force "additifs" Charmm et Amber, qui représentent la polarisabilité implicitement, donnent un moins bon accord. Ces calculs sont le premier exemple de simulations alchimiques d'énergies libre d'association relatives protéine:ligand avec un champ de force polarisable. Enfin, pour une modélisation future de peptides phophorylés, nous avons étendu le champ de force Drude pour inclure le méthyl phosphate et la phospho tyrosine. Il en résulte un excellent accord vaec les affinités expérimentales phosphate:magnésium.