Paramétrisation du champ de force CHARMM22 pour les ligands 8 et 308

Télécharger le fichier pdf d’un mémoire de fin d’études

Méthodes de simulation de dynamique moléculaire

L’ensemble des simulations de dynamique moléculaire présentée dans cette étude ont été réalisées en utilisant le programme CHARMM [41] avec l’algorithme d’intégration des équations newtoniennes du mouvement leapfrog [49] et un pas de temps d’une femtoseconde (fs). Afin d’utiliser un pas de temps d’une femtoseconde, les distances des liaisons covalentes formées avec des atomes d’hydrogène ont été fixées à leurs valeurs de référence (r0) en utilisant l’algorithme SHAKE [50]. Les contraintes appliquées par SHAKE permettent d’éliminer les fluctuations de ces liaisons covalentes qui interviennent à des fréquences élevées, i.e. dans des temps très courts (0,5 fs). En utilisant SHAKE, l’énergie correspondant aux liaisons covalentes formées avec des atomes d’hydrogène n’intervient pas dans l’énergie finale.
Parfois, il est nécessaire d’appliquer d’autres types de contraintes afin de fixer ou restreindre les positions des atomes ou les distances séparant deux atomes non-liés. Dans le programme CHARMM, la commande FIX permet de fixer la position des atomes en utilisant des constantes de force élevées. Pour restreindre les fluctuations autour d’une valeur de référence pour les distances séparant deux atomes non-liés, les contraintes de type NOE sont utilisées par le programme CHARMM. Ces contraintes permettent de définir les valeurs minimale et maximale acceptées pour la distance de séparation. Pour les contraintes NOE appliquées dans cette étude, les valeurs des paramètres KMIN, KMAX et FMAX sont respectivement égales à 100, 100 et 500 kcal/mol/Å2.
Lors des simulations de dynamique moléculaire, deux représentations du solvant peuvent être utilisées : représentation implicite et représentation explicite. Contrairement à la représentation explicite, dans une représentation implicite du solvant, les molécules d’eau ne sont pas incluses dans le système. Cette représentation implicite permet de tenir compte de certains effets relatifs au solvant (e.g., effet d’écran des charges atomiques et/ou friction entre soluté et solvant) tout en diminuant, en général, les temps de calcul nécessaires. Pour la représentation explicite du solvant, il est possible de diminuer les temps de calcul en limitant le nombre de molécules d’eau incluses dans le système.
Certains processus biologiques interviennent dans des intervalles de temps non accessibles par les méthodes classiques de simulation. Il est parfois nécessaire d’introduire une force dans le système afin de simuler ces processus biologiques dans des temps de calcul raisonnables. Dans les sous sections suivantes, les méthodes de simulation utilisées dans cette étude employant une solvatation implicite ou explicite et celles permettant d’introduire une force dans le système sont décrites.

Simulations de dynamique moléculaire aux limites stochastiques ou Stochastic Boundary Molecular Dynamics, SBD

La méthode Stochastic Boundary molecular Dynamics (SBD) [53, 35] a été imaginée dans le but de limiter le nombre de molécules du solvant contenues dans le système pour des simulations d’un processus localisé dans l’espace, ne faisant intervenir qu’une partie du système [54].
Cette méthode est particulièrement intéressante pour l’étude précise d’une région spécifique d’un système comme par exemple le site d’interaction protéine-ligand. Elle permet, lors des calculs, de ne pas prendre en compte la plupart des atomes n’intervenant pas directement dans le phénomène étudié. Elle est également très utile pour vérifier l’importance de certaines molécules d’eau dans la structure d’une protéine ou d’un complexe sans avoir à inclure de nombreuses molécules n’ayant pas la même importance. La méthode SBD a été maintes fois utilisée pour l’étude de systèmes biologiques et il a été montré qu’elle permet une description correcte des mouvements et interactions intervenant au niveau du site actif [55].
La méthode SBD repose sur un découpage du système en trois régions : réaction, tampon et réservoir (voir figure 2.1). Les atomes de la région de réaction sont soumis à une dynamique classique sans contrainte tandis que les mouvements de ceux de la région tampon sont régis par l’équation de Langevin (équation 2.10) e. La région réservoir entoure la zone d’interaction (régions tampon et de réaction) et fournit un champ de force statique où les atomes sont contraints, ce qui permet de maintenir les propriétés d’équilibre des atomes de cette zone. Des molécules d’eau f sont incluses seulement dans la zone d’interaction et peuvent diffuser librement entre les régions de réaction et tampon [35]. Afin d’éviter l’apparition de problèmes liés aux effets de bord, il est nécessaire de définir avec précaution les parties du système contenues dans les différentes régions.

Simulations de dynamique moléculaire modifiées par ajout d’une force ou Biased Molecular Dynamics, BMD

Les processus biologiques d’intérêt, comme la formation/dissociation de complexe dans les conditions physiologiques, sont associés à des barrières énergétiques relativement élevées. Par conséquent, ils surviennent sur des intervalles de temps de la milliseconde à la seconde, voire la minute [56]. Comme ces intervalles de temps ne sont pas accessibles par des simulations de dynamique moléculaire classique, la force quadratique supplémentaire dépendante du temps implémentée dans la méthode Biased Molecular Dynamics (BMD) [36, 37] est introduite afin d’expulser le ligand du site de fixation de la protéine. Dans cette méthode la force dépend du changement d’une coordonnée de réaction (RC) choisie et est introduite dans le système seulement lorsque cette coordonnée de réaction diminue. Contrairement aux autres méthodes, le temps de décomplexation est aléatoire puisque la force n’est pas introduite continuellement, comme attendu lors d’une décomplexation naturelle [37]. De plus, l’introduction de la force implémentée dans la méthode BMD permet d’obtenir une décomplexation qui intervient à travers de mécanismes similaires à ceux observés lors de la décomplexation naturelle, tel que démontré dans l’étude de Curcio et al. [57].

Calculs de chimie quantique

Afin de déterminer des paramètres pour les ligands 8 et 308 dans le champ de force CHARMM22, des calculs ab initio ont été nécessaires pour obtenir les paramètres de référence. Pour le ligand 8, tous les calculs ab initio, i.e. HF/3-21G*, HF/6-31G* et MP2/6-31G*//HF/6-31G*, ont été effectués en utilisant le programme Gaussian03 [58], tandis que pour le ligand 308, les programmes Gaussian03 et GAMESS-US [59] ont été utilisés afin de comparer les résultats obtenus avec ces deux programmes. Afin de construire les fichiers d’entrée nécessaires aux calculs ab initio, le programme MOLDEN [60] a été utilisé. Afin de déterminer les contributions des interactions électroniques à courte distance dans l’état intermédiaire, des calculs d’orbitales moléculaires ont été effectués avec l’hamiltonien PM5 du programme MOZYME disponible dans la version commerciale de MOPAC (MOPAC2006) [61]. Afin de diminuer les temps de calcul, les molécules d’eau ont été remplacées par l’approche continue COSMO [62] en utilisant une constante diélectrique de 78,4.

Procédure générale de paramétrisation du champ de force CHARMM22

La procédure de paramétrisation d’une nouvelle molécule dans un champ de force se déroule en quatre étapes principales (voir figure 3.1).
La première étape de la paramétrisation correspond à la recherche des paramètres manquants pour la molécule dans le champ de force CHARMM22. Par exemple, dans la molécule M de la figure 3.1 qui n’a pas été paramétrée dans CHARMM22, il est possible de trouver des groupes qui possèdent des correspondances avec une partie d’une molécule déjà paramétrée. Lors du développement de CHARMM22, les développeurs ont favorisé la possibilité d’utiliser pour le groupe d’atomes CH2SH des paramètres déjà existants tels que ceux de la chaîne latérale de la cystéine, ceci au détriment de la précision de ces paramètres pour la molécule M entière. Ainsi, il est possible de diminuer le nombre de paramètres à déterminer en utilisant de telles correspondances entre groupes d’atomes. En revanche, la partie HCOCO de la molécule M ne présente aucune correspondance avec une partie d’une molécule déjà paramétrée dans CHARMM22. Il est donc nécessaire de déterminer de nouveaux paramètres pour cette partie. Cette première étape de la paramétrisation conduit à définir des molécules modèles telles que D de la figure 3.1 qui vont représenter les parties pour lesquelles des correspondances ne peuvent être trouvées. Ces molécules modèles correspondent à des fragments de la molécule de départ dont les parties manquantes sont remplacées par des atomes d’hydrogène ou des groupements méthyle (voir figure 3.1). Cette fragmentation s’effectue de sorte à minimiser la taille des molécules modèles. En fait, plus la taille de la molécule modèle est importante, plus il est difficile d’obtenir des données expérimentales pour l’ensemble de la molécule ou d’effectuer des calculs ab initio dans des temps raisonnables (ceux-ci dépendant également des ressources disponibles). Toutefois, selon le niveau de précision souhaité, cette fragmentation ne doit pas se faire au détriment du respect de l’environnement chimique autour des parties pour lesquelles les paramètres doivent être déterminés. Par exemple pour la molécule M, en utilisant la fragmentation présentée sur la figure 3.1, il n’est pas nécessaire de déterminer de nouveaux paramètres pour le terme Einteractions liées pour la partie CH2-SH, du fait de la correspondance de ce groupe d’atomes avec la chaîne latérale de la cystéine. Toutefois, l’environnement chimique du grou-
1. Recherche des paramètres à déterminer
– Recherche de correspondance
– Définition de molécules modèles
2. Initialisation des nouveaux paramètres
– Assignation des types atomiques
– Attribution des valeurs initiales par analogie
3. Détermination des paramètres de référence
– Calculs au niveau ab initio ou mesures sur la structure cristallographique
4. Optimisation des nouveaux paramètres
– Optimisation des paramètres des termes non-liés
– Optimisation des paramètres des termes liés

Procédure d’optimisation des nouveaux paramètres pour CHARMM22

La procédure d’optimisation, décrite par MacKerell [67], est un processus itératif qui nécessite la vérification de tous les paramètres après modification de l’un d’entre eux. Lors de la procédure d’optimisation, les paramètres du champ de force sont ajustés de manière à obtenir une convergence des données obtenues par mécanique moléculaire (MM) vers celles de référence. Différents critères de convergence sont utilisés et dépendent du paramètre optimisé (critères présentés dans les sous-section suivantes).
Les procédures d’optimisation de chaque paramètre ainsi que les critères de convergence sont ensuite détaillés.

Optimisation des paramètres de Einteractions non-liées

Après l’obtention des géométries initale et de référence pour les molécules modèles, les premiers paramètres à être optimisés sont ceux des termes Evan der Waals et Eélectrostatiques. Ces paramètres sont toujours les premiers à être optimisés. Les méthodes utilisées pour optimiser les charges partielles atomiques et les paramètres de van der Waals sont présentées ci-dessous accompagnées des critères de convergence usuels.
Paramétrisation des charges partielles atomiques Afin d’optimiser les valeurs des charges partielles atomiques, deux méthodes sont principalement utilisées [67, 65]. La méthode utilisée dans cette étude repose sur la reproduction en mécanique moléculaire des énergies et distances d’interaction calculées en HF/6-31G* entre une molécule modèle et une molécule d’eau [70, 71]. Pour générer ces énergies d’interaction en HF/6-31G* et en mécanique moléculaire (MM), il est nécessaire de placer une molécule d’eau en interaction avec l’atome dont la charge partielle doit être paramétrée respectivement dans la géométrie de référence et dans la géométrie minimisée b. L’optimisation de chaque charge partielle atomique étant associée à l’utilisation d’une molécule d’eau placée en acceptrice ou en donneuse de liaison hydrogène ainsi que l’illustre la figure 3.2 qui supperpose deux exemples d’un tel placement. Après le placement initial de la molécule d’eau (dinitial), l’énergie d’interaction HF/6-31G* entre la molécule d’eau et la molécule modèle est calculée pour chaque valeur de distance de séparation comprise dans un intervalle choisi. Les géométries de la molécule modèle et de la molécule d’eau ne sont pas optimisées durant ces calculs, un scan dit rigide de l’énergie d’interaction est réalisé. À partir du modèle moléculaire équivalent construit en MM, il s’agit ensuite de reproduire la courbe d’énergie d’interaction HF/6-31G* c, plus particulièrement les valeurs d’énergie et de distance associées au minimum de cette courbe [67]. Afin de pouvoir être utilisé comme profils de référence, les valeurs obtenues par des calculs HF/6-31G* doivent être corrigées. Premièrement, les valeurs d’énergie obtenues pour les espèces non chargées doivent être multipliées par un facteur de 1,16 d. Ce facteur est utilisé pour compenser l’absence de traitement spécifique de la polarisation et conserver un équilibre entre les interactions soluté-solvant et solvant-solvant lors de l’utilisation du modèle d’eau TIP3P [67] e. Deuxièmement, pour l’interaction entre la molécule d’eau et un site polaire ne portant pas de charge nette, la distance de séparation associée à une énergie minimale (doptimal) est diminuée de 0,20 Å. Cette diminution permet de reproduire correctement les propriétés du solvant pur [73] et de tenir compte du manque de traitement explicite de la polarisation ainsi que de l’absence du traitement de la force attractive de dispersion dans la méthode HF [74].

Paramétrisation du champ de force CHARMM22 pour les ligands 8 et 308

Les ligands 8 et 308 (figure 3.3) comportent des parties pour lesquelles les paramètres nécessaires aux calculs de l’énergie potentielle sont manquants dans le champ de force CHARMM22g. Afin de pouvoir étudier les complexes FKBP12-8 et FKBP12-308, il a donc été nécessaire de paramétrer, i.e. définir et optimiser, ces deux ligands en utilisant la procédure décrite par Foloppe et MacKerell [79, 67] et détaillée dans le chapitre précédent.
Dans ce chapitre, le choix des molécules modèles utilisées pour définir les nouveaux paramètres est commenté, puis la définition des paramètres est décrite. Ensuite, les résultats de la procédure d’optimisation de ces paramètres sont présentés pour chaque type de contribution à l’énergie potentielle (voir sous-section 1.1. p. 13). Enfin, les résultats des tests de minimisation énergétique qui ont permis de contrôler les géométries des ligands paramétrés sont analysés afin de valider les nouveaux paramètres.

Recherche des paramètres à déterminer pour les ligands 8 et 308

Dans un premier temps, des équivalences ont été recherchées pour les atomes, les liaisons covalentes, les angles de valence, les angles dièdre et chaque partie du ligand afin de limiter le nombre de paramètres à optimiser.
L’observation des structures des ligands 8 et 308 permet de faire ressortir un découpage en quatre parties de chaque ligand (figure 3.3). Ces deux ligands possèdent une partie centrale, appelée coeur et trois parties périphériques. Pour le ligand 8, les parties périphériques corresg pondent à Ph1 par une chaîne aliphatique, Ph2 et iPe ou pour le ligand 308, à Ethe, iBu et Tol h.
À partir de ces découpages, des correspondances apparaissent entre les parties périphériques et des molécules déjà paramétrées dans CHARMM22. Par exemple, la partie iBu de 308 correspond à la chaîne latérale de la leucine (figure 3.4 B). Les correspondances trouvées pour les deux ligands sont présentées sur la figure 3.4. Pour la partie de coeur, il est possible de voir que les correspondances ne sont pas aussi directes que pour les parties périphériques. De même, pour la partie iPe de 8 la correspondance n’est pas exacte (figure 3.4 A). Ainsi, il est nécessaire d’expliquer les raisons qui nous ont conduits à utiliser malgré tout ces correspondances ou à nous en inspirer pour certains paramètres, notamment les paramètres des termes Evan der Waals, Eliaisons covalentes, Eangle de valence et EUrey-Bradley.

Paramètres pour Eangles de valence et EUrey-Bradley

Pour les paramètres permettant de définir les angles de valence, comprenant les paramètres des termes Eangles de valence et EUrey-Bradley, la même démarche que pour les paramètres des liaisons covalentes a été utilisée. Ainsi, la plupart des paramètres ont été transférés à partir de paramètres existants dans CHARMM22 (voir tableau A.5 p. 201 en Annexes). Toutefois, excepté pour les angles dièdres des parties Ph1-C17 et Ph2 de 8 et iBu de 308, les valeurs du paramètre 0 ont été tirées des valeurs d’angle de valence mesurées dans la géométrie du ligand optimisée en HF/6-31G*. Lorsque les correspondances ne sont pas exactes (e.g., pour les parties de coeur), les originaux de choix effectués sont commentés dans la suite. De plus, l’absence de correspondance satisfaisante pour les angles de valence C3 −\S4 − C5 et N7 −\S8 − C17 du ligand 308 avec les molécules incluses dans CHARMM22 a conduit à la détermination de nouveaux paramètres.
Ligand 8 Aucune correspondance exacte n’a été trouvée dans CHARMM22 pour les angles de valence C8 −\C9 − O4, C8 −\C9 − C10, O3 −\C8 − C9 et N7 −\C8 − C9 (figure 3.4). En revanche, en comparant les paramètres issus de la LPDB avec des paramètres contenus dans CHARMM22 pour des molécules présentant des similitudes avec la fonction α-dicétone, des valeurs équivalentes ont été trouvées. Ainsi, les paramètres K de la LPDB ont pu être utilisés pour les quatre angles précédents.
De la même manière, bien qu’il y ait des similitudes entre un résidu proline et le cycle pipécolinyle C2N7, les paramètres K des angles de valence C6 −\N7 − C2 et C1 −\C2 − H ont été pris dans la LPDB afin de tenir compte de la différence de taille de ces cycles. Le paramètre K de l’angle de valence N7 −\C8 − O3 a également été pris dans la LPDB. En outre, l’absence d’équivalence exacte pour l’angle de valence H20 −\C15 − O1 a conduit à utiliser le paramètre K trouvé dans la LPDB pour cet angle.
Ligand 308 Les angles de valence du cycle C2SN7 ont été traités de trois manières différentes. Contrairement à ce qui a été fait pour 8, les paramètres K des angles de valence formés avec l’atome 308-N7 ont été copiés sur ceux trouvés dans la proline, l’utilisation de cette correspondance ayant été jugée nécessaire pour ce ligand. Afin de tenir compte de la différence entre des cycles non-aromatiques formés de cinq et de six atomes, les paramètres K pour les angles de valence C2 −\C3 − S4, H2 −\C3 − S4, H3 −\C3 − S4, S4 −\C5 − H4, S4 −\C5 − H5 et S4 −\C5 − C6 ont été copiés sur ceux du cyclohexane. Enfin, la présence de l’atome 308-S4 dans le cycle C2SN7, a conduit à la détermination des paramètres des termes Eangles de valence et EUrey-Bradley pour l’angle de valence C3 −\S4 − C5.

Table des matières

1 Introduction 
2 Méthodologie 
1. Simulation de dynamique moléculaire
1.1. Les champs de force
1.2. Méthodes de simulation de dynamique moléculaire
1.2.1. Simulations de dynamique moléculaire Langevin, LD
1.2.2. Simulations de dynamique moléculaire aux limites stochastiques ou Stochastic Boundary Molecular Dynamics, SBD
1.2.3. Simulations de dynamique moléculaire modifiées par ajout d’une force ou Biased Molecular Dynamics, BMD
2. Calculs de chimie quantique
3 Paramétrisation du champ de force CHARMM22 24
1. Procédure générale de paramétrisation du champ de force CHARMM22
1.1. Procédure d’optimisation des nouveaux paramètres pour CHARMM22
1.1.1. Optimisation des paramètres de Einteractions non-liées
1.1.2. Optimisation des paramètres de Einteractions liées
2. Paramétrisation du champ de force CHARMM22 pour les ligands 8 et 308
2.1. Recherche des paramètres à déterminer pour les ligands 8 et 308
2.1.1. Paramètres pour Evan derWaals
2.1.2. Paramètres pour Eélectrostatiques
2.1.3. Paramètres pour Eliaisons covalentes
2.1.4. Paramètres pour Eangles de valence et EUrey-Bradley
2.1.5. Paramètres pour Eangles dièdres et Eangles dièdres impropres
2.1.6. Résumé des nouveaux paramètres déterminés pour les deux ligands
2.2. Découpage des ligands en molécules modèles
2.3. Initialisation des paramètres à optimiser
2.4. Optimisation des nouveaux paramètres
2.4.1. Charges partielles atomiques optimisées
2.4.2. Paramètres optimisés pour les angles de valence
2.4.3. Paramètres optimisés pour les angles dièdres
3. Validation de la paramétrisation
4 Les états complexés FKBP12-8 et FKBP12-308 73
1. Procédure de simulation des états complexés EC1 et EC2
1.1. Préparations des structures initiales SI1 et SI2
1.2. Simulations des états complexés EC1 et EC2
1.2.1. Simulations LD
1.2.2. Simulations SBD
2. Analyse des états complexés
2.1. Analyse des structures de FKBP12 et du ligand dans EC1 et EC2
2.1.1. RMSD
2.1.2. Analyse des structures secondaires de FKBP12
2.1.3. Analyse des fluctuations des atomes de FKBP12 et du ligand dans EC1 et EC2
2.2. Analyse des interactions formées entre FKBP12 et le ligand dans EC1 et EC2
2.2.1. Interactions permanentes
2.2.2. Interactions transitoires entre deux groupes d’atomes
2.3. Conclusions des analyses de FKBP12 et du ligand dans EC1 et EC2
5 Protocole de décomplexation et étude des états intermédiaires 
1. Simulations de décomplexation et stabilisation des ensembles de structures métastables
2. Les simulations SBD comme dernière étape de la recherche des états intermédiaires
3. Caractéristiques des états intermédiaires
3.1. Structure des EI
3.2. Position du ligand dans le site de fixation de FKBP12
3.3. Mobilité de FKBP12 et du ligand
3.4. Analyse des interactions formées dans les EI
3.4.1. Analyse des interactions permanentes
3.4.2. Analyse des interactions transitoires
3.4.3. Analyse des orbitales moléculaires
3.4.4. Résumé des interactions formées dans les EI
6 Discussion 
1. Discussion des caractéristiques principales des EI
2. Transposition du modèle de l’état intermédiaire à d’autres ligands de FKBP12
3. Les EI permettent d’expliquer les différences d’affinité entre certains ligands de FKBP12
4. Importance des interactions de type π-π dans
la reconnaissance moléculaire
5. Propriétés structurales importantes du coeur des ligands de FKBP12
7 Conclusion 
Annexes 
Bibliographie 

Télécharger le rapport complet