Algorithme d'optimisation par colonie de fourmis

Faculté des Sciences Département d’Informatique

Algorithme d’optimisation par colonie de fourmis : développement et application à la prédiction ab initio de la structure native des protéines Fabian TEHEUX

Mémoire présenté sous la direction du Prof. Marianne ROOMAN et du Prof. Esteban ZIMANYI en vue de l’obtention du grade de Licencié en Informatique Année académique 2005–2006

Table des matières Préface 1

2

7

Introduction 1.1 Les protéines . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Structure primaire . . . . . . . . . . . . . . . . . . . 1.1.2 Structure secondaire . . . . . . . . . . . . . . . . . . 1.1.3 Structure tertiaire . . . . . . . . . . . . . . . . . . . . 1.2 Le repliement des protéines . . . . . . . . . . . . . . . . . . 1.3 La prédiction de la structure native des protéines . . . . . . 1.3.1 Objectif . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Aperçu général des méthodes de prédiction . . . . 1.3.3 Les impératifs de la prédiction de structures natives Présentation des outils et méthodes existants 2.1 Les modèles utilisés . . . . . . . . . . . . . . . . . . 2.1.1 Modélisation de la protéine . . . . . . . . . 2.1.2 Discrétisation de l’espace conformationnel 2.2 L’exploration de l’espace conformationnel . . . . . 2.2.1 Monte Carlo . . . . . . . . . . . . . . . . . . 2.2.2 Le Recuit Simulé . . . . . . . . . . . . . . . 2.2.3 Recherche Tabou . . . . . . . . . . . . . . . 2.2.4 Algorithmes évolutionnaires . . . . . . . . 2.2.5 Colonies de fourmis . . . . . . . . . . . . . 2.2.6 Méthodes hybrides . . . . . . . . . . . . . . 2.3 Les fonctions énergétiques . . . . . . . . . . . . . . 2.4 L’évaluation des résultats d’une prédiction . . . . 2.5 Récapitulatif . . . . . . . . . . . . . . . . . . . . . .

2

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . .

9 9 10 11 12 13 17 17 18 19

. . . . . . . . . . . . .

20 20 20 23 25 25 26 27 28 32 39 39 40 41

TABLE DES MATIÈRES 3

4

5

Nos outils et méthodes 3.1 Modélisation de nos protéines . . . . . . . . 3.2 Discrétisation de notre espace de recherche 3.3 Exploration de l’espace conformationnel . . 3.4 Fonction énergétique utilisée . . . . . . . . 3.4.1 Les potentiels de torsion . . . . . . . 3.4.2 Les potentiels de distance . . . . . . 3.4.3 Utilisation combinée des potentiels 3.5 Evaluation de nos solutions . . . . . . . . . 3.6 Récapitulatif . . . . . . . . . . . . . . . . . .

3

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

Notre algorithme : AC-ProPre 4.1 L’espace conformationnel . . . . . . . . . . . . . . . 4.2 Les potentiels énergétiques . . . . . . . . . . . . . . 4.3 La structure et les calculs effectués sur les protéines 4.3.1 La structure protéique tridimensionnelle . . 4.3.2 Les méthodes de calcul . . . . . . . . . . . . 4.4 La colonie de fourmis . . . . . . . . . . . . . . . . . . 4.5 Récapitulatif . . . . . . . . . . . . . . . . . . . . . . . Tests et résultats 5.1 Les configurations d’exécution de notre algorithme 5.2 Résultats obtenus . . . . . . . . . . . . . . . . . . . . 5.2.1 Evaluation énergétique de notre algorithme 5.2.2 Calcul du RMSD . . . . . . . . . . . . . . . . 5.2.3 RMSD Vs énergie libre . . . . . . . . . . . . . 5.3 Notre meilleure prédiction : une RMSD de 3,5222Å

. . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . .

42 43 45 46 48 48 49 50 51 51

. . . . . . .

52 52 53 54 55 56 59 60

. . . . . .

62 62 63 63 64 69 72

Conclusion

76

Bibliographie

78

Table des figures 1.1 1.2 1.3 1.4 1.5 1.6

2.1

2.2 2.3

2.4 2.5

Représentation atomique d’un acide aminé. . . . . . . . . . . . . . Représentation simplifiée d’un acide aminé. . . . . . . . . . . . . . Illustration des angles de torsion φ et ψ dûs aux liaisons chimiques. Exemples de structures secondaires : (a)hélice-α, (b)feuillet-β. . . Exemples de structure tertiaire : l’ubiquitine. . . . . . . . . . . . . Représentation tridimensionnelle du paysage énergétique d’une protéine : (a) un paysage énergétique idéal, (b) un paysage énergétique réel, (c) deux chemins de repliement différents, (d) le paradoxe de Levinthal. [source : Gilquin Bernard : Exploration des méchanismes de repliement des protéines par dynamique moléculaire] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Résultats d’une prédiction bidimensionelle basée sur le modèle HP pour une protéine de séquence ’HHHPPHPHPHPPHPHPHPPH’ (le début de la séquence est symbolisée par ’1’ et les résidus ’H’ sont en noir). . . . . . . . . . . . . . . . . . . . . . . . . . . . . Illustration de l’angle ω. La molécule entourée est un résidu quelconque venu s’attacher au résidu que nous étudions. . . . . . . . Exemple de modification conformationnelle locale : (a) l’état avant la modification conformationnelle, (b) l’état après la modification conformationnelle. . . . . . . . . . . . . . . . . . . . . . . . . . . . Fonctionnement d’un algorithme évolutionnaire ordinaire. . . . . Evolution du comportement des fourmis en fonction d’une modification de l’environnement : (a) situation initiale, (b) un ostacle s’interpose et les fourmis choisissent d’aller d’un côté ou de l’autre avec la même probabilité, (c) l’axe BCD étant plus court, les fourmis arrivent plus vite de l’autre côté de l’obstacle, (d) le nombre croissant de phéromones fait émerger le chemin BCD. . . . . . . . 4

10 11 12 12 13

16

21 22

24 30

33

TABLE DES FIGURES 3.1

5

Représentation simplifiée d’un acide aminé, attention le cas particulier de la Glycine verra ses atomes Cβ et Cµ se trouver au même endroit que l’atome Cα . . . . . . . . . . . . . . . . . . . . . . . . . . Exemple d’angle de valence : en pointillé, l’angle de valence associé à l’atome O. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

4.1

AC-ProPre : déroulement schématique. . . . . . . . . . . . . . . .

61

5.1

Graphiques de l’énergie moyenne des prédictions effectuées avec une configuration de (α = 1, β = 1, ρ = 0.8, ds0 et ds1 ne sont pas pris en compte). Les graphique allignés se passent à un nombre d’itérations équivalent et les graphiques en colonne ont le même facteur d’élitisme. Les quatres courbes représentent chacune un nombre de fourmis d’élite différent. . . . . . . . . . . . . . . . . . Graphiques de l’énergie moyenne des prédictions effectuées avec une configuration de (α = 1, β = 1, ρ = 0.6, ds0 et ds1 ne sont pas pris en compte). Dans ce cas-ci, les diverses courbe montrent un comportement tout à fait disparate. . . . . . . . . . . . . . . . . . . Graphiques de l’énergie moyenne des prédictions effectuées avec une configuration de (α = 1, β = 2, ρ = 0.8, ds0 et ds1 sont pris en compte). Les graphique allignés se passent à un nombre d’itérations équivalent et les graphiques en colonne ont le même facteur d’élitisme. Les quatres courbes représentent chacune un nombre de fourmis d’élite différent. . . . . . . . . . . . . . . . . . Graphiques de l’énergie moyenne des prédictions effectuées avec une configuration de (α = 1, β = 1, ρ = 0.6, ds0 et ds1 sont pris en compte). Dans ce cas-ci, les diverses courbe montrent un comportement tout à fait disparate. . . . . . . . . . . . . . . . . . . Graphiques du RMSD moyen des prédictions effectuées avec une configuration de (α = 1, β = 5, ρ = 0.8, ds0 et ds1 sont pris en compte). L’allure de ces courbes n’évoque rien de particulier. . . . Graphiques de comparaison entre le RMSD et l’énergie libre d’une conformation prédite : (a) prise en compte de toute les prédictions effectuées sans considération des potentiels de distance, (b) limitation aux prédictions présentant un RMSD inférieur à 7Å. . . . . Représentation de la structure tertiaire de la prédiction effectuée par AC-ProPre. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Représentation de la structure tertiaire de l’ubiquitine (résidu 1-35).

3.2

5.2

5.3

5.4

5.5

5.6

5.7 5.8

43

65

66

67

68

70

71 73 74

TABLE DES FIGURES 5.9

Supperposition de la conformation simulée par AC-ProPre avec la conformation native connue. . . . . . . . . . . . . . . . . . . . .

6

75

Préface Depuis ses origines, l’informatique a apporté des solutions dans des domaines aussi variés que complexes. Aujourd’hui, la science y fait une fois de plus appel afin de prédire la structure tridimentionelle des protéines. Dans la vie de tous les jours, les protéines jouent des rôles importants sur notre Terre. Il est donc très intéressant de comprendre comment elles fonctionnent. Malheureusement, il nous faut pour ce faire déterminer la structure tridimentionnelle qu’une protéine adoptera afin de jouer son rôle biologique. Cela n’est évidemment pas aisé ! Les scientifiques se sont donc rués sur des méthodes expérimentales ayant fait leurs preuves, mais ces méthodes sont d’une lenteur que le plus patient des curieux ne saurait supporter. Ainsi, de nombreux chercheurs utilisent aujourd’hui les moyens informatiques afin de prédire ce que l’expérience pourra démontrer. Lorsqu’un jour des méthodes de prédictions seront parfaitement au point, il sera alors possible de créer les médicaments les plus adaptés aux défaillances diverses de l’organisme, d’utiliser au mieux les diverses caractéristiques de ces protéines afin d’embellir notre environnement, etc. Un jour peut-être sera-t-il même possible de créer une protéine en fonction de besoins du moment ! Nous voicis donc devant une frontière scientifique que l’exploration ne fera que repousser... En ce qui nous concerne, nous ferons un premier pas dans l’implémentation d’une méthode de prédiction (basée sur la métaheuristique d’optimisation par colonie de fourmis) appliquée à un espace discret mais proche de celui dans lequel évoluent les protéines réelles. En effet, la littérature actuelle propose quelques articles implémentant ce type de méthode mais elle n’a encore jamais été utilisée avec des modèles de représentations proche de la réalité. Nous découvrirons donc pas à pas toutes les informations nécessaires à la bonne compréhension du fonctionement de cette méthode pour ensuite en découvrir les particularités. 7

« Z  ,      ,              . » (« U   ,     `  ´     ´        ’. ») Gerhard MULDER

Chapitre 1 Introduction 1.1 Les protéines Nous sommes en 1835. Un chimiste hollandais, Gerhard Mulder, effectue une analyse chimique sur le blanc d’œuf et déclare que c’est une substance constituée d’une base de carbone (symbole chimique C), hydrogène (H), oxygène (O) et azote (N), auxquels s’ajoutent un peu de phosphore et de souffre. Il montre également que des substances d’une composition identique sont présentes dans des tissus animaux et végétaux. Les protéines venaient d’être découvertes. Le mot protéine vient du grec prôtos qui signifie premier, essentiel. En effet, les protéines interviennent dans la grande majorité des processus qui régissent le fonctionnement de tout être vivant [1]. Leurs rôles sont aussi variés que complexes : – l’insuline, par exemple, joue un rôle hormonal et transmet des messages au travers de l’organisme ; – l’hémoglobine, protéine cinématographiquement célèbre, s’occupe du transport du dioxygène dans notre corps ; – d’autre protéines de transport existent et agissent dans des processus tels que la photosynthèse ; – certaines protéines jouent un rôle majeur dans la conversion d’énergie chimique en énergie mécanique, dans les muscles par exemple ; – les anticorps dont nous avons tous entendu parlé et qui jouent un rôle prédominant dans le système immunitaire sont encore des protéines. Ils savent faire la différence entre les corps connus et les corps intrus tels que les virus et autres bactéries.

9

CHAPITRE 1. INTRODUCTION

10

Une autre éthymologie, moins probable mais quelque peu plus séduisante, voudrait que le mot protéine fasse référence au dieu grec Protée, dieu qui possédait un pouvoir de polymorphisme infini sur son corps. Un fois de plus, force est de constater que s’il existe des corps naturels présentant une multitude de formes différentes, ce sont bien les protéines !

1.1.1 Structure primaire Toutes les protéines ont un point commun : ce sont des polymères linéaires construits à partir de différentes combinaisons des 20 acides aminés de base. On peut représenter un acide aminé comme à la figure 1.1. L’atome central est appelé carbone α (Cα ). Trois structures chimiques s’y connectent : le groupement amine (−NH2 ), le groupement carboxyle (−COOH) et une chaîne latérale (habituellement symbolisée par R). Chaque acide aminé diffère par cette chaîne latérale. De plus, cette dernière confère à l’acide aminé ses propriétés particulières.

F. 1.1 – Représentation atomique d’un acide aminé. Dans la suite de ce mémoire, nous simplifierons les acides aminés par la représentation suivante (cfr. figue 1.2) où N, Cα , C et O sont les atomes réels de la chaîne principle de l’acide aminé, Cβ est le premier atome de la chaîne latérale et Cµ est un atome fictif utilisé pour simuler la chaine latérale. Le Cµ correspond à la moyenne géométrique de tous les atomes de la chaîne latérale de l’acide aminé (à l’exception des atomes d’hydrogène H). Cette moyenne est calculée à partir d’une base de donnée des structures protéiques connues dont on a extrait les informations de l’acide aminé concerné [2]. Cette représentation simplifiée a (évidemment) une exception : l’acide aminé nommé glycine. En effet, cet acide aminé ne possède qu’un simple atome d’hy-


11

F. 1.2 – Représentation simplifiée d’un acide aminé. drogène en guise de chaîne latérale. Il est donc commun d’établir une correspondance des atomes Cβ et Cµ avec l’atome Cα . Lors de la formation d’une protéine, la liaison entre deux acides aminés se fait au moyen d’une réaction de condensation. Ainsi, il est commun de parler de résidus au lieu d’acides aminés. Le nombre de résidus constituant une protéine varie d’une cinquantaine à plusieurs milliers. On appelle séquence d’une protéine la suite des résidus qui la compose. C’est cette séquence qui constitue la structure primaire de la protéine. Il est remarquable de constater que la séquence contient toutes les informations nécessaires à l’adoption d’une structure spécifique et à l’éxcution de la fonction biologique de la protéine.

1.1.2 Structure secondaire Les atomes sont liés entre eux par un lien chimique et le lien qui lie l’atome C d’un résidu à l’atome N d’un autre résidu est appelé lien peptidique. De par leurs natures, ces liaisons offrent une certaine liberté de rotation aux atomes liés. Toutefois, le lien peptidique entre deux résidus est particulier : il a un caractère partiel de double liaison qui a pour conséquence de contraindre dans un même plan les atomes Cα de deux résidus voisins et les atomes intermédiaires C, O, N. Cela dit, comme nous le montre la figure 1.3, il existe toujours des libertés de rotation autour des liens N − Cα (appelé angle de torsion φ) et Cα − C (angle ψ). De plus, dans le cas d’une représentation de résidus plus détaillée, les chaines latérales jouissent également d’une certaine liberté de rotation. Il en résulte que le nombre de conformations potentiellement accessibles pour une protéine d’une


12

certaine longueur est astronomiquement grand. Et pourtant, certains motifs structuraux réguliers sont observés de manière récurente dans les protéines. Ce sont ces arrangements que l’on appelle structure secondaire. Parmi eux, deux catgories se distinguent : les hélices-α (figure 1.4(a)) et les feuillets-β (figure 1.4(b)). Le squelette (backbone) de ces structures est composé de la suite des atomes Cα de chaque acide aminé constituant la protéine.

F. 1.3 – Illustration des angles de torsion φ et ψ dûs aux liaisons chimiques.

F. 1.4 – Exemples de structures secondaires : (a)hélice-α, (b)feuillet-β.

1.1.3 Structure tertiaire Chaque protéine a la particularité de posséder des caractéristiques physiques et chimiques qui lui sont propres. Toutefois, ces caractéristiques sont spécifiques à la conformation spatiale adoptée par la protéine. La majorité des protéines se replient donc pour adopter une conformation tridimensionnelle unique appelée


13

structure tertiaire dont la particularité est d’utiliser des interactions entre résidus proches dans l’espace mais distants dans la séquence. Cette structure tertiaire, également appellée état natif, est stable grâce aux nombreuses interactions favorables qui s’établissent en son sein. Notons que l’on appelle interacions natives une interaction présente dans la structure native. La figure 1.5 nous montre en exemple la structure tertiaire de l’ubiquitine. Les diverses structures secondaires sont reliées par ce que l’on appelle des coudes.

F. 1.5 – Exemples de structure tertiaire : l’ubiquitine. C’est grâce à l’adoption de cette structure tridimensionnelle particulière que les protéines peuvent réaliser leur fonction biologique. L’hémoglobine, par exemple, transporte le dioxygène dans une cavité enfouie au sein de sa structure tertiaire.

1.2 Le repliement des protéines Une protéine peut adopter différentes conformations, celles de l’état dénaturé étant généralement les moins stables. A l’image d’un ressort trop tendu, une protéine en état instable cherche à atteindre son état le plus stable caractérisé par une énergie libre minimale. Une hypothèse largement admise est que la structure native d’une protéine correspond généralement à sa conformation d’énergie libre minimale. Le repliement d’une protéine est un processus extrêmement spécifique et particulièrement efficace. En effet, la rapidité avec laquelle une protéine trouve sa conformation de plus basse énergie est remarquable. Par exemple, une protéine


14

de 100 résidus peut adopter 2100 (≈ 1030 ) structures tridimensionnelles différentes en supposant que seulement 2 conformations soient accessibles à chaque résidu, ce qui est loin de la réalité. En admettant que le passage d’une conformation à une autre se fait en 10−13 secondes (ce qui correspond au temps necessaire pour une rotation autour d’une liaison peptidique), il faudrait 1017 secondes soit environ trois milliard d’années, pour tester toutes les conformations possibles. Très clairement, la totalité des conformations accessibles à une véritable protéine, sans aucune simplification, dépasse de loin le degré d’abstraction humain. Pourtant, les protéines arrivent à retrouver leur structure native en un temps allant de quelques millisecondes à quelques secondes ! Levinthal en 1969 [3], fut le premier à relever cette incompatibilité que l’on nomme aujourd’hui paradoxe de Levinthal. Il ajoutait que, de toute évidence, les protéines n’explorent pas la totalité de leur espace conformationnel. Actuellement, seuls des modèles théoriques viennent expliquer la rapidité de ce processus. En effet, les mécanismes qui permettent à une séquence d’acides aminés de se replier en une structure tridiensionelle unique sont actuellement mal connus. Le modèle de diffusion-collision propose un mécanisme de repliement hiérarchique qui permet de réduire drastiquement l’espace conformationnel. La théorie soutient qu’il y aurait, dans un premier temps, un repliement des sousdomaines de la protéines, ces sous-domaines étant d’une taille très limitée. Par la suite, des interaction entre sous-domaines se créeraient, formant ainsi des domaines plus important. Cela durerait jusqu’à la formation complète de la structure native de la protéine. Le modèle de l’effondrement hydrophobe quant à lui s’appuie sur la tendance des résidus hydrophobes à se regrouper dans les premiers instants du processus afin d’éviter tout contact avec l’eau. Ce regroupement aurait lieu avant la formation de structures secondaires et aurait une influence prépondérante pour la suite du repliement. En effet, le noyau hydrophobe temporairement formé réduirait considérablement l’espace conformationnel accessible à la protéine. Enfin, le modèle de nucléation-condensation met en avant l’importance des interactions tertiaires (c’est-à-dire entre résidus distants dans la séquence). Afin que ces interactions aient lieu, la formation d’un noyau de repliement est requise. Les fragments de protéine non repliés entrant en contact avec le noyau adopteraient ensuite leur structure native. D’autre modèles théoriques existent et le lecteur intéressé trouvera une brève présentation de la majorité d’entre-eux dans [4]. Durant les années 90, un nouveau modèle de repliement protéique a été


15

développé [5] et fut revisité en 1998 [6]. Sa caractéristique principale est de représenter le paysage énergétique d’une protéine. Il s’agit en fait d’une représentation multi-dimensionnelle de l’énergie libre d’une conformation de la protéine en fonction de sa similarité avec la structure native. Une interprétation tridimensionnelle du paysage énergétique produit un entonnoir communément appelé entonnoir de repliement. L’axe vertical de cet entonnoir représente l’énergie libre d’une conformation et les axes horizontaux représentent les conformations de la protéine accessibles à un certain niveau d’énergie. A la base de cet entonnoir, on trouve la structure native qui est donc d’une énergie libre minimale. Plus on monte dans l’entonnoir, plus le nombre de conformations possibles à un même niveau d’énergie augmente. Le repliement serait donc un processus progressif de modifications conformationnelles entraînant la création d’interactions natives favorables et débouchant sur la structure native (cfr. figure 1.6). Dans une situation idéale, l’entonnoir serait parfait et chaque modification conformationnelle en direction de l’état natif stabiliserait la protéine (figure 1.6(a)). Malheureusement, l’expérience montre que des modifications conformationnelles nécessaires à la création d’interactions natives peuvent être défavorables d’un point de vue énergétique ! Cela introduit des variations dans le paysage énergétique et laisse apparaitre un concept bien connu dans le domaine de l’informatique : le problème des minima locaux (figure 1.6(b)). Un paysage énergétique peut donc s’avérer plus ou moins complexe, ce qui induit qu’il existe différents chemins de repliement entre une conformation d’état dénaturé et une conformation d’état stable (figure 1.6(c)). Il est intéressant de comparer le paradoxe de Levinthal à un paysage énergétique plat [7]. L’idée à l’époque était qu’il existait une conformation native unique et que toute autre conformation était inintéressante. On se rend alors compte que le problème d’origine est mal posé : en effet, il est mathématiquement improbable de tomber dans le trou de l’entonnoir à partir d’un paysage plat en un temps raisonnable (figure 1.6(d)). Ce modèle de paysage énergétique ne remet pas en question tous les modèles précédents, au contraire. Il permet d’expliquer plusieurs phénomènes observés tels que l’existence de chemins de repliement préférés. Cela appuierait certaines théories comme par exemple celles de la création d’un noyau hydrophobe ou d’intermédiaires de repliement. Concrètement, le paysage énergétique procure un cadre général dans lequel les cas particuliers peuvent être interprétés.


16

F. 1.6 – Représentation tridimensionnelle du paysage énergétique d’une protéine : (a) un paysage énergétique idéal, (b) un paysage énergétique réel, (c) deux chemins de repliement différents, (d) le paradoxe de Levinthal. [source : Gilquin Bernard : Exploration des méchanismes de repliement des protéines par dynamique moléculaire]


17

1.3 La prédiction de la structure native des protéines 1.3.1 Objectif La structure native d’une protéine est une richesse en soi. Elle est l’inconnue nécessaire à la résolution d’une foule de questions, dont la principale est : « Comment ça marche ? ». En effet, pour comprendre le fonctionnement d’une protéine il faut impérativement en connaître la structure native. Par exemple, les sites actifs des protéines ne sont généralement fonctionnels qu’une fois la protéine complètement repliée. Outre cet aspect primordial, de multiples possibilités s’offrent à nous : – l’étude des interactions possibles entre une protéine et diverses molécules permet par exemple la création de médicaments qui agiront précisément là où c’est nécessaire ; – certaines protéines subissent des mutations pour diverses raisons et la compréhension fine de l’effet de ces mutations passe par la connaissance de la structure native ; – la conception de protéines dont la séquence est modifiée permettrait d’obtenir de nouvelles propriétés comme par exemple une augmentation de la stabilité de la protéine ou encore l’amélioration de réactions enzymatiques dans le milieu industriel ; – la compréhension des mécanismes de stabilisation de l’état d’une protéine pourrait conduire, dans un futur plus ou moins proche, à la création de nouvelle protéines avec un rôle bien précis et défini à l’avance. Le besoin sera alors à l’origine de la protéine... Il existe une base de données accessible à tous, la Protein Data Bank (PDB) [8], qui répertorie toutes les connaissances actuelles concernant les structures des protéines. On y retrouve de nombreuses structures natives établies expérimentalement par cristallographie aux rayons X ou par résonance magnétique nucléaire. Pour plus d’information dans ce domaine, le lecteur intéressé se dirigera vers [9, 10]. Malheureusement, ces méthodes, en plus d’être lentes, sont plutôt coûteuses. De plus, le nombre de séquences protéiques connues est bien plus grand que le nombre de structures natives établies et l’écart entre les deux ne cesse d’augmenter. Ces divers arguments ont poussé les chercheurs à se tourner vers d’autres méthodes. L’expérience ayant montré que la plupart des protéines adoptent des confor-


18

mations spécifiques, l’étude in silico des protéines va bon train. De nombreux projets tels que l’étude théorique du repliement des protéines ou la simulation des interactions protéine-protéine ont vu le jour. Parmi eux, on retrouve la prédiction de la structure native des protéines.

1.3.2 Aperçu général des méthodes de prédiction Avant de découvrir ce qu’est la prédiction, une question se pose : pourquoi ne pas simuler intégralement un repliement protéique ? Après tout, les connaissances actuelles permettent d’envisager une simulation complète, picomètres par picomètres, des mouvements atomiques d’une protéine... Actuellement, ce genre de technique existe mais il est calculatoirement impossible d’effectuer une simulation longue de plus de quelques nanosecondes (au mieux une microseconde). Or, nous le découvrions tout à l’heure, un repliement protéique nécessite quelques millisecondes (au pire secondes). L’échelle de temps n’est donc tout simplement pas accessible et par ailleurs, l’informatique a ce défaut de toujours procéder par approximation de valeurs, ce qui peut entrainer sur le long terme une dégradation fortement prononcée des résultats. Il en va d’ailleurs de même pour les évaluations énergétiques qui sont certainement imparfaites [11]. L’intérêt de la prédiction est donc évident : gagner du temps ! Deux grandes catégories de prédictions existent : celle basée sur les connaissances actuelles et l’autre, plus pure, se basant uniquement sur la structure primaire de la protéine, c’est-à-dire sa séquence. La première catégorie repose sur l’utilisation de bases de données de structures protéiques. Il faudra donc qu’il existe un certain degré de similitude de séquence ou de structure entre la protéine étudiée et une ou plusieurs protéines de la base de donnée. La modélisation par homologie, par exemple, compare les séquences de la PDB avec la séquence protéique cible [12]. Si un degré d’homologie suffisant est dégagé, alors on produit une structure dérivée des structures homologues de la PDB. Ce procédé donne de bons résultats car les protéines de séquence similaire ont généralement des structures très semblables. Un autre exemple est la reconnaissance structurale [13, 14]. On l’utilise généralement quand la modélisation par homologie a échoué. En effet, un tel échec est clairement dû à un manque de similitude entre la protéine cible et les protéines de la PDB. Cette méthode propose donc d’utiliser un squelette issu des conformations structurales connues afin d’y enfiler la séquence d’acides aminés


19

de la protéine cible. A nouveau, il existe des bases de données structurales. Par exemple, la Fold Classification based on Structure-Structure alignement of Proteins (FSSP) [15] regroupe les structures représentatives de la PDB. L’intérêt d’une telle banque provient du fait que la PDB contient des structures redondantes. Dans le cadre de la reconnaissance structurale, l’utilisation de structures représentatives est un atout significatif. L’autre catégorie de prédiction se fait ab initio, c’est-à-dire sur base de la seule séquence de la protéine. Cette catégorie est beaucoup plus générale mais est bien évidemment plus ardue à mettre en œuvre. L’hypothèse de base est celle selon laquelle la structure native correspond au minimum global d’énergie libre. Le but d’une prédiction ab initio est donc évident : trouver une structure tridimensionnelle d’énergie libre minimale.

1.3.3 Les impératifs de la prédiction de structures natives Afin de réaliser au mieux une prédiction, certains outils sont nécessaires : 1. une modélisation de la protéine étudiée ; 2. une méthode de discrétisation de l’espace dans lequel elle évolue ; 3. une méthode de recherche qui permet l’exporation de l’espace conformationnel ; 4. une fonction énergétique permettant l’évaluation d’une structure par rapport à sa séquence ; 5. dans certains cas, il sera intéressant de posséder d’un moyen d’évaluation des prédictions obtenues. Le chapitre suivant propose donc un tour d’horizon de ces divers impératifs. Nous y découvrirons plusieurs outils et méthodes permettant de procéder à une prédiction digne de ce nom.

Chapitre 2 Présentation des outils et méthodes existants Ce chapitre propose au lecteur un aperçu des quelques techniques utilisées dans le domaine de la prédiction de structures natives. Il permettra en outre de mieux se rendre compte des divers problèmes rencontrés dans ce domaine et de découvrir les solutions qui y sont apportées.

2.1 Les modèles utilisés Une difficulté récurrente en informatique est celle de la modélisation des problèmes étudiés. Dans ce cas-ci comme dans beaucoup d’autres, il faut trouver un juste compromis entre simplification du problème et qualité des résultats. Dans le cadre des prédictions de la structure native des protéines, il y a deux entités à modéliser : la protéine elle même et l’espace de recherche.

2.1.1 Modélisation de la protéine La modélisation d’une protéine peut se faire par une représentation complète au niveau atomique des acides aminés la composant. Evidemment, un tel modèle entraînerait un coût de calcul non négligeable. D’un autre côté, une représentation trop simpliste de la protéine entraînera une dégradation de la qualité des résultats. Il faut donc trouver un compromis dont on pourra tirer le meilleur parti dans les deux cas.

20

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 21 Le modèle HP Le modèle HP [16] est simple : il propose de réduire les vingt acides aminés à seulement deux groupes : les acides aminés hydrophiles (nommé P) et les acides aminés hydrophobes (nommé H). Dès lors, une protéine, plutôt que d’être représentée par une séquence composée de 20 acides aminés différents, sera simplement représentée par une suite d’acides aminé H ou P, ce qui simplifie vraiment les choses ! De plus, un résidu est réduit à un seul point dans l’espace, on ne tient donc absolument pas compte du détail atomique. Ce modèle étant basé sur la théorie de la création d’un noyau hydrophobe, il faudra, pour satisfaire aux conditions de repliement des protéines, minimiser l’interaction des résidus H avec l’environnement extérieur (généralement l’eau). Le résultat d’une prédiction effectuée à partir d’un modèle de ce genre se trouve en figure 2.1.

F. 2.1 – Résultats d’une prédiction bidimensionelle basée sur le modèle HP pour une protéine de séquence ’HHHPPHPHPHPPHPHPHPPH’ (le début de la séquence est symbolisée par ’1’ et les résidus ’H’ sont en noir). Malheureusement, une représentation aussi simpliste fourni peu de résultats intéressants dans le cadre de la prédiction de structure native de protéine. En effet, réduire les propriétés des acides àminés à un simple caractère d’hydrophobie-hydrophilie est assez limitatif. Par contre, au niveau des études théoriques, elle permet de mieux comprendre certaines étapes du repliement des protéines. Le modèle (φ, ψ, ω) Comme nous l’avons vu dans l’introduction, les résidus peuvent adopter diverses conformations en fonction des angles de torsion adoptés autour de leurs liaisons chimiques. Bien que l’angle ω puisse sembler nouveau au lecteur,

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 22 il s’agit en fait du nom donné à celui de la liaison inter-résidu C − N. Pour rappel, ce lien à un caractère partiel de double liaison, ce qui limite à deux les valeurs prises par l’angle ω : 0◦ et 180◦ .

F. 2.2 – Illustration de l’angle ω. La molécule entourée est un résidu quelconque venu s’attacher au résidu que nous étudions. Ainsi, la conformation d’un résidu est donnée par un triplet (φ, ψ, ω) et le squelette d’une protéine peut être modélisé par une suite de triplets, de taille équivalente au nombre de résidus de la protéine. Finalement, la structure tridimensionnelle d’une protéine sera calculée à partir des conformations (φ, ψ, ω) successives adoptées pas les résidus la composant. Les modèles intermédiaires On entends par modèle intermédiaire un modèle offrant une description de la protéine plus détaillée que dans le modèle HP, mais moins détaillée qu’un modèle complet. Ainsi, il arrive de regrouper les acides aminés en plusieurs classes afin de permettre un niveau de détail plus élevé. Une fois cette classification effectuée, on peut encore se contenter de représenter le résidu par un seul point dans l’espace ou bien par quelques atomes. A ce sujet, il est important de se souvenir d’une simplification que nous avons effectuée lors de l’introduction : le remplacement de la chaîne latérale par

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 23 un atome fictif, le Cµ . Dans le cas ou la chaîne latérale serait prise en compte dans le modèle, une représentation encore plus fine de la protéine aurait lieu. En effet, cette chaîne latérale possède également des libertés de rotation autour de certaines liaisons chimiques et il serait tout à fait envisageable de les prendre en compte.

2.1.2 Discrétisation de l’espace conformationnel Nous l’avons déjà dit, l’espace conformationnel est immensément grand. Diverses techniques de discrétisation existent donc afin de permettre une diminution substanciellment importante des conformations possibles. Bien entendu, il s’agit à nouveau de trouver un bon compromis entre degré de simplification et qualité de résultat. Deux grands types de discrétisation existent : 1. le modèle réseau propose une restriction de l’espace conformationnel en maille fixe ; 2. le modèle hors réseau propose une restriction de l’espace conformationnel en limitant par exemple les rotations autour des liaisons chimiques à certaines valeurs angulaire bien définies, ce qui a pour conséquence de créer un modèle réseau dynamique. Les modèles réseau Les premiers modèles réseaux étaient basiques et consistaient en une simple grille bidimensionnelle fixe. Le lecteur comprendra aisément qu’une telle représentation protéique est bien loin de la réalité. Grâce à l’amélioration des techniques de recherche dans l’espace conformationnel, des modèles réseaux plus fins commencèrent à être utilisés. L’utilisation de grille tridimensionnelle à maille étroite offre une nette amélioration de représentation mais le caractère statique de ces modèles représente un facteur limitatif important au niveau de la qualité des solutions obtenues. En effet, dans un modèle réseau, on travaille avec des coordonnées fixes définies par la maille utilisée. Ainsi, dans le cas d’une représentation atomique de la protéine, les distances inter-atomiques seront fixes et les possibilités de rotation angulaire seront également définies (par exemple, dans un modèle bidimensionnel carré où la maille représente la distance inter-atomique, ces rotations seront toujours multiples de 90◦ ).

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 24 Toutefois, certains avantages existent de par la rigidité du modèle. Par exemple, il est tout à fait possible qu’un résidu passe d’une conformation à une autre sans que cela nous oblige à recalculer la position de tous les autres résidus (cfr. figure 2.3). C’est là l’avantage de travailler avec des coordonnées externes au modèle de protéine utilisé.

F. 2.3 – Exemple de modification conformationnelle locale : (a) l’état avant la modification conformationnelle, (b) l’état après la modification conformationnelle.

Les modèles hors réseau Les modèles hors réseau sont généralement basés sur la limitation des conformations accessibles pour un résidus, c’est-à-dire sur la limitation des valeurs prises par les angles de rotation autour de leurs liaisons chimiques. Par exemple, la liaison N − Cα sera limitée à un ensemble de quelques angles φ accessibles. Par ailleurs, la distance des liaisons inter-atomiques peut être définie pour chaque paire d’atome [17], ce qui accentue encore le réalisme de ce modèle. La discrétisation du modèle dépendra donc du nombre de valeurs disponibles pour chaque résidu. Ce système simple est facilement transposable dans le cas où la chaine latérale des protéines est prise en compte. Son avantage est évident : il est basé sur la nature même des protéines et utilise donc un système de coordonnée internes pour la représentation de cellesci. De plus, il permet de travailler dans un espace conformationnel accessible à la prédiction tout en préservant une certaine qualité de résultat. Cependant, pour un même niveau de discrétisation, le modèle hors réseau demande généralement plus de calculs que le modèle réseau. Pour reprendre l’exemple du modèle réseau, si nous modifions la conformation d’un seul résidu de la protéine dans le cadre du modèle hors réseau, il nous faudra recalculer

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 25 toutes les coordonnées tridimensionnelles de la protéine à partir des triplets (φ, ψ, ω) de chaque résidu. Il s’agit en fait d’effectuer la traduction des coordonnées internes en coordonnées réelles.

2.2 L’exploration de l’espace conformationnel De par son paysage énergétique, une protéine ne peut parfois atteindre son minimum global d’énergie libre qu’en procédant à certaines modifications conformationnelles défavorables. Autrement dit, il arrive parfois que pour passer d’une état X à un état Y plus stable, la protéine se retrouve dans un état intermédiaire de stabilité moindre que son état d’origine X. C’est là le grand problème des minima locaux. Précisément, le problème est de détecter que la protéine se trouve dans une conformation d’énergie libre minimale locale qui n’est pas native. Le problème sous-jacent est évidemment de trouver un moyen de sortir de ce minimum local, c’est-à-dire d’accepter une modification conformationnelle défavorable. Diverses solutions ont été apportées à ce problème, mais aucune d’entre-elles n’est optimale. On parle donc couramment d’heuristiques et de métaheuristiques. Une heuristique est une méthode de production de solution adaptée à un problème particulier alors qu’une métaheuristique constiste en une stratégie de choix pouvant piloter une heuristique. En résumé, un bon compromis entre l’exploration de l’espace conformationnel et l’intensification de la recherce d’une solution minimale doit être trouvé. Si la priorité est donnée à l’exploration, il est fort probable qu’on ne converge jamais vers une solution minimale (même locale) alors que si l’on privilégie l’intensification d’une solution, on aura de grandes chances de rester bloqué dans un minimum local.

2.2.1 Monte Carlo Monte Carlo est une méthode de prédiction probabiliste basée sur la modification aléatoire de la conformation de certains résidus de la protéine étudiée [18]. A chaque itération de la méthode, une nouvelle structure protéique peut donc être évaluée par la fonction énergétique. La suite de l’algorithme dépendra de l’énergie calculée. Dans les première versions de Monte Carlo, le choix était simple :

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 26 – si l’énergie calculée est inférieure à l’énergie de la conformation précédente, alors c’est que la nouvelle conformation offre à la protéine un état plus stable. On conservera donc cette conformation qui servira de base à la prochaine itération ; – dans l’autre cas, on procède simplement au rejet de la nouvelle conformation. Ce genre de rejet automatique entraîne fatalement la découverte de solutions locale, ce qui n’est pas vraiment l’objectif recherché. Ainsi, certains critères d’acceptation ont été développés. Le concept est commun : permettre l’acceptation de modifications conformationnelles défavorables afin de diminuer les chances de tomber dans un minimum local. Un critère régulièrement utilisé est celui dit de Metropolis [19] où ∆E représente la différence d’énergie entre la conformation évaluée et la dernière conformation acceptée : ≤ e−(∆E)(kT) Le facteur kT représente la multiplication entre la constante de Boltzman et la température à laquelle s’effectue la prédiction et est un nombre aléatoire compris entre 0 et 1 [19]. Ainsi, si le critère de Métropolis est vérifié, alors une modification conformationnelle défavorable sera acceptée et servira de base à la prochaine itération. On constate donc qu’une première approche d’exploration de l’espace conformationnel passe simplement par la considération d’un critère d’acceptation stochastique. En effet, plus la modification conformationnelle est défavorable, plus ∆E est grand, ce qui implique que l’exponentielle est plus petite que 1 et donc que les chances d’acceptations sont réduites. Par contre, dans le cas d’une modification conformationnelle favorable, ∆E sera négatif et l’exponentielle sera toujours plus grande que 1. La conformation sera donc toujours acceptée.

2.2.2 Le Recuit Simulé La méthode de Monte Carlo incluant un critère de Métropolis est attrayante mais peut encore être améliorée. En effet, il y a un paramètre intéressant sur lequel nous pouvons jouer : la température. C’est donc en 1983 qu’apparaît le Recuit Simulé. Il propose simplement de commencer la prédiction à haute température et de la faire baisser progressivement [20]. Dans la vie courrante, on constate par exemple que la qualité d’un rail de chemin de fer est fortement

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 27 dépendante de son refroidissement lors de sa confection. S’il est refroidi trop vite, la structure métallique sera imparfaite et le rail présentera des défauts. Il faut donc prendre le temps, ce qui rend généralement cette méthode moins rapide que les autres. L’introduction d’une variation progressive de la température lors d’une prédiction de structure native introduit donc un facteur dynamique intéressant dans le critère d’acceptation. Celui-ci permettra un degré d’acceptation de modification conformationnelle défavorable plus élevé au début de la prédiction et plus faible à la fin. Il favorise donc une phase d’exploration dans un premier temps pour permettre ensuite une intensification de la recherche. Il arrive parfois d’effectuer un Recuit Simulé sur une protéine déjà partiellement repliée. Il est donc important de remarquer que le choix de la température initiale est un facteur très important : si une température trop élevée est choisie, on pourrait assister à un dépliement de la protéine, ce qui n’est pas du tout le but recherché ! Enfin, signalons que sous certaines conditions (dont par exemple un schéma de décroissance particulier), des preuves de convergence asymptotique existent pour le Recuit Simulé. Ceci constitue un avantage mathématique inexistant dans le cas des méthodes heuristique [21, 22]. Malheureusement, il est très difficile d’appliquer au mieux la théorie sans consacrer beaucoup de temps au paramétrage de l’algorithme.

2.2.3 Recherche Tabou En 1986, Fred Glover propose une nouvelle méthode concernant le problème des minima locaux : l’utilisation d’une liste taboue simulant en fait le principe de mémoire [23]. L’application de cette méthode donne une heuristique en ce sens qu’elle ne converge pas spécialement vers le minimum global. Elle est par contre un excellent guide ! L’idée de la recherche Tabou est de permettre la poursuite d’une recherche même quand un minimum local est trouvé. La mémoire sera utilisée pour éviter les retours en arrière et donc les cycles. Il est important de signaler qu’une mémoire élevée deviendra trop restrictive tout comme une mémoire trop petite sera probablement inutile. L’intérêt typique de cette mémoire est de forcer la recherche à explorer une autre zone que celle dans laquelle se trouve la solution la plus récente. Formellement, définissons s comme étant la solution actuelle, sG comme étant

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 28 la solution golbale et V(s) comme étant l’espace de solutions voisines à s. Une fois qu’une solution s est trouvée, elle est mise en mémoire (dans la liste Tabou) pour un nombre d’itération i. Il sera dés lors impossible de revenir à cette solution tant que le nombre i d’itération ne sera pas dépassé. Toutefois, une exception existe : lorsqu’un mouvement tabou permet d’améliorer la solution s, il sera autorisé. C’est ce que l’on appelle le critère d’aspiration. Le grand avantage de la recherche Tabou est que lorsqu’elle devra choisir un voisin dans V(s), elle choisira toujours le meilleur alors que pour le Recuit Simulé, on assistait à une modification aléatoire de la solution. On découvre ainsi une autre approche de résolution du problème des minima locaux : plutôt que de les éviter, on les utilise Enfin, notons qu’il existe diverses variantes de la recherche Tabou utilisant des modification sur la fonction d’évaluation des voisins. L’idée générale est d’ajouter des pénalités ou des avantages sur base de certains critères remplis par le voisin évalué. Une autre variante de la recherche Tabou alterne volontairement les phases d’intensification de la recherche avec les phases d’exploration de l’espace en changeant dynamiquement le nombre de mouvements tabous stockés en mémoire.

2.2.4 Algorithmes évolutionnaires Les diverses méthodes que nous venons d’explorer ont toutes un point commun : leurs itérations reposent sur l’amélioration d’une solution unique. En effet, on travaillera toujours à partir d’une solution initiale en essayant de l’améliorer au maximum en fonction des contraintes du problème étudié. Les méthodes que nous allons explorer à partir de maintenant sont des méthodes à population de solutions. Le therme algorithme évolutionnaire regroupe en fait une multitude de métaheuristique différentes. Cependant, toutes ont un point commun : elles utilisent le concept d’algorithme génétique, c’est-à-dire l’application informatique simplifiée de la théorie de l’évolution de Darwin. L’actuel succès des algorithmes génétiques est évidemment dû aux ouvrages de David Goldberg [24] mais leur date de naissance réelle provient de la publication en 1975 du livre de John Holland intitulé Adaptation in Natural and Artificial System. Les algorithmes évolutionnaires sont de par leur nature de métaheuristique extrèmement flexibles. Ils peuvent donc être appliqués à de nombreux problèmes. Toutefois, leur utilisation en tant que « magic box »est à proscrire ! En

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 29 effet, certains mécanismes internes nécessitent un coût calculatoire important et les performances d’une implémentation irréfléchie de ce genre de méthode seraient plus que médiocres. De manière générale, il est important de spécialiser l’algorithme en fonction du problème étudié. En effet, l’inclusion de connaissance au processus de recherche ne fera qu’améliorer la qualité des résultats obtenus. Toutefois une chose est certaine, un algorithme évolutionnaire ne fournira jamais de meilleures performances qu’une méthode de recherche exacte (si elle existe). Avant d’examiner le fonctionnement de ce type de méthode, définissons-en le vocabulaire. L’espace de recherche sera appelé l’environnement et les solutions seront les individus regroupés en population. Le but de l’algorithme est d’évaluer la capacité d’une population à s’adapter à son environnement. Pour mesurer la qualité des individus, l’algorithme utilise une fonction d’évaluation appellée fonction de fitness. Le but devient donc la recherche d’un individu maximisant la fonction de fitness. Pour représenter un individu, les algorithmes génétiques utilisent des chromosomes. On dira que le génotype d’un individu représente son codage quand son phénotype sera l’expression du code génétique dans l’environnement. Les différent types d’algorithmes évolutionnaires sont en fait caractérisés par le moyen de codage des individus utilisé. Maintenant que ce vocabulaire est acquis, penchons-nous sur le fonctionnement de l’algorithme avant d’approfondir les méthodes d’exploration de l’espace et d’intensification des solutions. Comme le montre la figure 2.4, tout commence au départ d’une population initiale. L’idée des algorithmes évolutionnaires est de fonctionner par génération de population. Il faut donc créer de nouveaux individus sur base de la population courante. Pour ce faire, une partie de la population sera sélectionnée pour donner naissance à de nouveaux individus. Ces individus enfants seront créés par croisement et/ou mutation afin former la nouvelle population. Par ailleurs, afin de conserver une certaine mémoire des individus visités, un facteur d’élitisme sera utilisé pour inclure certains parents dans la nouvelle génération d’individus. On répetera ces opérations jusqu’à ce qu’un critère de terminaison soit atteint. Bien entendu, s’il existe de multiples algorithmes évolutionnaires, c’est non seulement parce que le codage des individus diffère de méthode en méthode mais également parce que les opérations de sélection, de croisement et de mutations connaissent beaucoup de variantes. Une foule de possibilités existe mais nous ne les explorerons pas ici. Le lecteur intéressé en trouvera quelques

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 30

F. 2.4 – Fonctionnement d’un algorithme évolutionnaire ordinaire.

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 31 exemples particulièrements intéressants dans [25]. C’est par contre le moment de pointer du doigt les mécanismes qui nous intéressent ! En effet, ce sont ces trois opérateurs (sélection, croisement et mutation) qui constituent la base du compromis entre exploration et intensification que le présent ouvrage s’attèle à déterminer ! Par ailleurs, la méthode de reproduction ainsi que le remplacement de la population sont deux autres facteurs intéressants. On peut dire que l’opérateur de sélection a une influence sur la convergence de l’algorithme. Ainsi, afin de garder un dergé d’exploration suffisant, une certaine diversité doit être conservée parmis les individus. Il faudra donc que la population de parents sélectionnés ne soit pas trop restreinte. L’opérateur de croisement quant à lui est le principal outil mettant en œuvre le compromis entre exploration et intensification. Cet opérateur plus que n’importe quel autre doit être spécifique au problère étudié ainsi qu’au codage utilisé dans la représentation des individus. En ce qui concerne l’opérateur de mutation, on peut dire qu’il garanti une certaine exploration aléatoire de l’espace. Il faudra donc veiller à ne pas trop l’utiliser afin de ne pas entrer dans un état de stagnation dû à l’éparpillement constant des individus. Vient ensuite la méthode de reproduction. On entend par là le nombre d’enfants généré par un couple d’individus parents. De multiples possibilités existent, qui influe généralement sur la rapidité de concentration de la population autour des bons individus. Pour finir, il reste le facteur de remplacement de la population. Si une population est entièrement renouvelée à chaque itération, les chances de convergence se verront réduites. Il faut adopter un certain niveau d’élitisme qui permet de conserver une mémoire des bons individus rencontrés. Comme le montrent les trois pararaphes précédents, le nombre de choix paramétriques à faire est assez important. C’est là le grand problème des algorithmes évolutionnaires : il faut du temps pour dégager une configuration paramétrique intéressante. On constate donc que la famille des algorithmes évolutionnaires dispose de plusieurs outils intéressants jouant directement sur les individus de la population. Remarquons que la qualité de la population initiale n’est pas sans importance ! En effet, elle est un facteur déterminant pour le déroulement et la convergence de l’algorithme. Il est d’ailleurs fréquent d’utiliser comme population initiale les résultats de diverses recherches locales effectuée préalablement. Cela donne naissance à des algorithmes hybrides que nous aurons l’occasion de

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 32 découvrir plus loin. Pour terminer cette section sur les algorithmes évolutionnaires, notons qu’ils offrent facilement la possibilité d’effectuer une exécution parallèle sur un cluster de machines ou sur une machine multi-processeur. En effet, dès lors que chaque individu doit être évalué individuellement, une évaluation simultanée de plusieurs individus est tout à fait envisageable ! Une autre possibilité facilement exploitable serait de conserver la population d’origine et de procéder par sélection inverse : un individu fils viendrait remplacer le plus mauvais individu de la population.

2.2.5 Colonies de fourmis Les fourmis dans la nature Tout comme les algorithmes évolutionnaires, la méthode des colonies de fourmis est basée sur une observation naturelle : le comportement des fourmis dans leur cohabitation en colonies. En 1983, des chercheurs découvrent que les fourmis ont la capacité de toujours trouver le chemin le plus court entre leur nid et une source de nourriture (par exemple) [26], ce qui est assez surprenant car la grande majorité des fourmis est aveugle. Par ailleurs, les fourmis sont également aptes à réagir aux modifications de l’environnement. Ainsi, lorsqu’un ancien chemin est obstrué par un obstacle quelconque, les fourmis trouveront un autre chemin, toujours plus court que tous les autres possibles. Ces particularités surprenantes ont donc fait l’objet de nombreuses recherches et finalement, un moyen de communication chimique entre les fourmis fut découvert : la phéromone. Concrètement, les fourmis déposent régulièrement une certaine quantité de phéromone sur le chemin qu’elles parcourent. Lorsqu’elles arrivent à un point où plusieurs chemins s’offrent à elles, elles font un choix probabiliste sur base de la quantité de phéromone présente sur chaque chemin. Ainsi, lorsque plusieurs fourmis sillonnent le même espace, on assiste à l’émergence du meilleur chemin entre le nid et la source de nourriture. Pour mieux comprendre cette particularité, examinons ce qu’il se passe dans le cas d’une modification environnementale. Initialement, les fourmis utilisent un chemin donné (cfr. figure 2.5(a)) entre le nid (A) et la source de nourriture (E). Subitement, une bûche bloque le chemin d’origine. Les foumis sont donc contraintes de trouver un autre chemin. Dans le sens A vers E, les premières fourmis à devoir faire ce choix auront

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 33 la même probabilité de se diriger à gauche ou à droite de l’obstacle car aucune autre fourmi ne sera passée avant elles. Il n’y a donc aucune trace phéromonale (cfr. figure 2.5(b)). Toutefois, une fourmi ayant choisi d’aller à droite (chemin BCD) arrivera en D bien avant une fourmi ayant choisi la gauche. De plus la situation est identique en ce qui concerne le sens E vers A (cfr. figure 2.5(c)). Finalement, la quantité de fourmis utilisant l’axe BCD est plus élevée que pour l’axe BHD ce qui entraîne une quantité de phéromone croissante sur l’axe BCD. Après un certain laps de temps, toutes les fourmis devraient (statistiquement) utiliser l’axe BCD (cfr. figure 2.5(c)). C’est ce phénomène que l’on appelle la stigmergie.

F. 2.5 – Evolution du comportement des fourmis en fonction d’une modification de l’environnement : (a) situation initiale, (b) un ostacle s’interpose et les fourmis choisissent d’aller d’un côté ou de l’autre avec la même probabilité, (c) l’axe BCD étant plus court, les fourmis arrivent plus vite de l’autre côté de l’obstacle, (d) le nombre croissant de phéromones fait émerger le chemin BCD.

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 34 La transposition informatique des fourmis Comme nous venons de le voir, les fourmis possèdent naturellement un outil leur permettant de communiquer de façon générale afin de faire ressortir un comportement global de la colonie. Les premières publications de cette idée datent de 1991 [27, 28] et la première application fut faite sur le problème du voyageur de commerce. Cette transposition était évidente vu que l’objectif de ce problème est de trouver un chemin minimal répondant à certains critères. Il est toutefois possible d’utiliser un système de fourmis pour résoudre d’autres problèmes. Les méthodes de recherche par colonie de fourmis on une particularité : la recherche de solutions s’effectue sur base de la totalité de la population et non plus a partir de quelques meilleurs individus (comme c’est le cas dans les algorithmes évolutionnaires). Cette particularité non négligeable a donné naissance la création de nombreux types d’algorithme de recherche, nous en découvrirons ici les grandes lignes. Le concept général est de simuler le comportement de plusieurs agents collaborants dans la recherche d’une meilleure solution en utilisant un moyen de communication simple : la phéromone. Dans le cadre de la métaheuristique d’optimisation par colonie de fourmi, les agents sont ... des fourmis. Bien entendu, il y aura quelques différences entre les fourmis naturelles et celles utilisée dans le domaine informatique que nous nommerons désormais fourmis artificielles. Ainsi, une fourmi artificielle : 1. n’est pas totalement aveugle et peut donc voir un peu plus loin que son entourage direct ; 2. évolue dans un univers ou le temps est discret ; 3. dispose d’une certaine mémoire lui permettant de retenir la solution qu’elle a construit ; 4. gère la quantité de phéromone déposée en fonction de la qualité de la solution et peut procéder à différents types de dépôts phéromonaux. Ces quelques améliorations ne sont pas sans conséquences, nous allons le découvrir ci-dessous.. Pour commencer, découvrons l’impact engendré par l’acuité visuelle de nos fourmis artificielles. Très clairement, cela permet à la fourmi de se guider autrement que par son seul environnement direct. Ce guide est en fait une valeur que nous appellerons valeur heuristique, c’est-à-dire une valeur pouvant influencer le

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 35 choix fait par notre fourmi. Par exemple, dans le cadre de la recherche d’un plus court chemin, lorsqu’une fourmi devra faire un choix sur le chemin à suivre, une valeur heuristique pouvant lui simplifier la tâche serait, par exemple, la longueur des chemins accessibles. Cette valeur heuristique est donc une donnée propre au problème qui permet de spécialiser la métaheuristique. Il s’agit là d’un atout important car, nous l’avons déjà vu, utiliser une métaheuristique comme méthode de recherche universelle n’est pas une bonne idée et conduirait certainement à de médiocres résultats. C’est par l’inclusion de connaissances propres au problème étudié que la méthode de recherche fournira les meilleurs résultats. Revenons en donc à nos fourmis. Chacune tente de trouver une solution minimisant la fonction d’évaluation utilisée. Cette solution est donc construite pas-à-pas en fonction de divers paramètres parmis lesquels nous retrouvons la valeur heuristique. Bien évidemment, les phéromones sont toujours de la partie. Pratiquement, le choix fait par une fourmi artificielle dépends de ces deux données. De manière plus formelle, on dit qu’à chaque point nécessitant un choix, une table de décision, résultat d’une opération mathématique particulière, produit les règles probabilistes sur lesquelles sont basés les mouvements des fourmis. Afin de découvrir cette opération mathématique, notons par τi j la quantité de phéromone présente sur un chemine reliant i à j et par ηi j la valeur heuristique associée à ce même chemin. En considérant que la fourmi se trouve sur un point i, la probabilité pi j qu’elle se dirige vers le point j sera donnée par : pi j = P

[τi j ]α · [ηi j ]β k∈V(i) [τik ]

α

· [ηik ]β

où V(i) représente le voisinage de i, c’est-à-dire tous les points accessibles à partir de i et où α et β sont des paramètres permettant de contrôler l’importance relative donnée aux phéromones et aux valeurs heuristiques. Pour être complet, rappelons nous que nos fourmis artificielles évoluent dans un univers ou le temps est discret. La formule devient donc : pi j = P

[τi j (t)]α · [ηi j (t)]β β α k∈V(i) [τik (t)] · [ηik (t)]

De ces formules, nous retirons q’une fourmi procède à la construction d’une solution étape par étape en effectuant un choix stochastique en fonction des phéromones et des valeurs heuristiques présentes sur son chemin. Nous venons

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 36 donc de mettre le doigt sur deux des paramètres les plus fondamentaux de la méthode des colonies de fourmis : α et β. En effet, ces paramètres ont un degré d’influence important par rapport à l’exploration de l’espace conformationnel et à l’intensification des solutions. Nous aurons l’occasion d’y revenir au chapitre suivant lorsque nous décrirons la manière dont nous avons adapté la méthode des colonies de fourmis au problème de la prédiction de structures natives. Maintenant que nous savons comment s’orientent les fourmis, nous allons découvrir comment se gère le dépôt de phéromones. En effet, c’est en déposant des phéromones qu’une fourmi augmente la probabilité de fréquentation du chemin qu’elle utilise. Pour en revenir à notre petite énumération des différences entre fourmi réelle et fourmi artificielle, nous allons découvrir simultanément l’intéret de la mémoire de la fourmi artificielle ainsi que les diverses façons dont elle dépose ses phéromones. Grâce à sa mémoire, une fourmi peut construire une solution valable répondant à certains critères dépendant du problème étudié. De plus, une fois la solution élaborée, elle peut l’évaluer dans sa totalité ! C’est grâce à ces capacités de souvenir et d’évaluation qu’une fourmi peut procéder à différents types de dépôts de phéromones : – lorsqu’une fourmi dépose des phéromones tout au long de l’élaboration de sa solution, on parlera de dépôt en ligne, pas-à-pas (online step by step) ; – par contre, si une fourmi procède à un dépôt phéromonal à la fin de la construction de sa solution, on parlera de dépôt en ligne, retardé (online delayed). De plus, la quantité de phéromone déposée sera directement proportionnelle à la qualité de la solution (partielle) élaborée. Si cette solution est intéressante, la fourmi déposera une grande quantité de phéromone afin de maximiser les chances de passage par les autres fourmis. Nous venons donc de découvrir que l’intensification d’une solution se fait suite à la quantité phéromonale présente sur un chemin. Bien entendu, la valeur donnée à α reste prépondérante. Dans le cas ou celle-ci serait de 0, les phéromones ne seraient tout simplement pas pirses en compte. Il est par ailleurs intéressant de constater qu’une réduction néfaste de l’espace de recherche aura lieu si la quantité phéromonale déposée est mal gérée. Dans ce cas, on pourrait assister à l’émergence rapide d’un chemin minimal qui engendrerait un état de stagnation de la recherche. En effet, ce chemin étant utilisé par un nombre croissant de fourmis, la quantité phéromonale ne cesserait d’augmenter et plus aucune exploration de l’espace n’aurait lieu.

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 37 Heureusement, une phéromone n’est pas éternelle ! En effet, dans la nature tout comme dans la transposition informatique, les phéromones s’évaporent. Nous introduisons donc un facteur de persistance qui sera appliqué à tout dépôt phéromonal. En nottant par ρ le facteur de persistance, nous pouvons dire que (1 − ρ) correspond à l’évaporation. Ainsi, un moyen de contrôle de la convergence de l’algorithme nous est donné. Ce facteur de persistance évoque d’une certaine manière la mémoire de la colonie. Pour faire la comparaison avec les algorithmes génétiques, la persistance est à la phéromone ce que l’élitisme était à la population. Attention toute fois à ne pas se méprendre, l’exemple d’émergence rapide d’un chemin minimal sera toujours possible ! La métaheuristique des colonies de fourmis, comme bien d’autres, demande un temps considérable dans la recherche d’une paramétrisation satisfaisante. Une solution intéressante à ce problème consiste à modifier le paramètre α pour quelques tours, afin de relancer l’exploration de l’espace. Un autre problème potentiel résolu par le facteur de persistance est ce que l’on appelle en anglais le phénomène d’attrition. Ce problème provient de l’accumulation de phéromones due à un tronçon partiel commun dans plusieurs solutions différentes. Toutefois, il faut remarquer que ce phénomène se produit plus facilement dans un espace bidimentionnel fortement discrétisé que dans un espace tridimentionel tel que ceux dont nous avons parlé ci-dessus. Les cinq derniers paragraphes nous ont donc appris que chaque fourmi fait un dépot de phéromone en fonction de la solution qu’elle a calculé et qu’un facteur de persistance ρ permet de gérer l’évaporation des dépots de phéromones. Il est donc intéressant de formaliser cela. Ainsi, en notant par ∆τkij (t, t + n) la quantité de phéromone déposée par la fourmi k sur le chemin i j dans l’interalle de temps (t, t + n) nous obtenons que : τi j (t + n) = ρτi j (t) +

m X

∆τkij (t, t + n)

k=1

où m représente le nombre total de fourmis, n est l’intervalle unitaire de temps discret et τi j (t) est toujours la quantité de phéromone présente sur le chemin reliant i à j. On peut traduire cette formule en disant que la quantité de phéromones présentes sur le chemin i j à l’instant (t + n) est égale à la persistance de la quantité de phéromone présente à l’instant précédent (t) à laquelle on ajoute les nouveaux dépots de toute les fourmis ayant parcouru le chemin i j depuis l’instant précédent (t).

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 38 Enfin, pour clôturer cette section, notons que lorsqu’une fourmi a procédé à l’élaboration d’une solution ainsi qu’à la mise à jour phéromonale correspondante, elle meurt. Cela nous permet de bien nous rendre compte de l’indépendance des fourmis par rapport à la colonie. En effet, tout au long de sa vie, la fourmi n’aura utilisé que les informations locales qui se trouvaient sur son chemin, sans oublier la valeur heuristique. Les particularités offertes par la colonie artificielle Les avantages présentés par la métaheuristique des colonies de fourmis ne s’arrêtent pas là ! En effet, en plus des améliorations effectuées sur les fourmis, il est possible d’effectuer des actions centralisées qu’une fourmi réelle serait incapable de réaliser. Par exemple, il est possible de procéder à un dépot phéromonal supplémentaire pour certaines solutions jugées particulièrement intéressantes. On parlera alors de dépôt hors-ligne (offline), ce qui correspond à une certaine forme d’élitisme anticipatif. Nous pouvons donc récapituler le fonctionnement d’une colonie de fourmi par le pseudocode suivant : while ( c r i t e r e de t e r m i n a i s o n non a t t e i n t ) { planification { fourmis_en_action ( ) ; evaporation_des_pheromones ( ) ; actions_generales ( ) ; } }

Il est important de remarquer que les opérations de planification ne sont pas définies, ce qui permet une multitude d’implémentation de cette méthode [29]. De plus, la métaheuristique, de par son fonctionnement, permet des implémentations faciles pour des systèmes informatiques en cluster [30, 29]. Enfin, notons que d’autres méthodes heuristiques basées sur des observations naturelles existent (comme par exemple les réseaux de neurones). Le lecteur intéressé sera certainement comblé par [22] qui propose un excellent tour d’horison en la matière.

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 39

2.2.6 Méthodes hybrides Il arrive fréquemment que des solutions hybrides soient apportées au problème d’exploration de l’espace. L’idée est de combiner deux méthodes complémentaires : une de recherche globale et une de recherche locale. En effet, les effets combinés obtenus par ce genre d’hybridation sont souvent encourageants. A titre d’exemple, examinons l’hybridation d’un algoritme évolutionnaire [31] avec, comme méthode de recherche locale la méthode Tabou (mais nous aurions pu choisir n’importe quelle autre autre). Typiquement, la phase d’initialisation de l’algorithme évolutionnaire sera modifiée en ce sens qu’une recherche locale sera effectuée sur chaque individu de départ. Ensuite commencera seulement la reherche proprement dite. De plus, plutôt que de procéder simplement à un croisement lors du renouvellement de la population, une intensification locale sera à nouveau effectuée sur les individus résultant du croisement. Ces modifications permettent de concentrer la recherche globale autour de solutions à priori locales. Bien entendu, tout autre mélange est possible. Il faudra toutefois procéder à l’hybridation de manière intelligente afin de ne pas engendrer de surcoût de calcul ni de dégradation de compromis entre exploration et intensification. D’autres exemples d’hybridation existent en masse, [32] propose par exemple l’hybridation d’un réseau de neurone (métaheuristique non explorée dans le cadre de cet ouvrage) avec la recherche Tabou.

2.3 Les fonctions énergétiques Quelle que soit la recherche effectuée in silico, un besoin commun est la fonction énergétique capable d’évaluer la correspondance qualitative et quantitative entre une séquence d’acides aminés et une conformation spatiale donnée. En effet, la performance atteinte par le matériel et la qualité des algorithmes utilisés ne servent à rien sans une fonction d’évaluation pertinente. Les fonctions énergétiques fournissent des résultats sur base de potentiels qui sont représentatifs des interactions présentes au sein d’une structure tridimensionnelle. Lorsqu’une interaction interne est native, son potentiel énergétique sera faible. Bien entendu, le nombre de types d’interactions dans une protéine est élevé, ce qui engendre une multitude de potentiels calculables. A titre d’exemple, le potentiel d’hydrophobicité est calculé sur base des différences entre résidus hydrophobes et résidus hydrophiles et est totalement indépendant de la conformation locale adoptée par la protéine. Son utilisation

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 40 influencera la prédiction vers un arrangement spécifique des résidus hydrophobes. Autre exemple, le potentiel de distance qui est calculé sur base de la propension d’un couple de résidus distants dans la séquence, à être séparés dans l’espace par une certaine distance. Nous aurons l’occasion de revenir sur ces potentiels lors de l’explication détaillée du calcul de l’énergie libre de nos prédictions.

2.4 L’évaluation des résultats d’une prédiction Afin de développer des algorithmes toujours plus performant, de nombreux tests sont effectués. Afin de procéder à ces mesures, il est fréquent de faire de la prédiction de structure native sur des protéines dont on connaît déjà la conformation d’énergie libre minimale. Pour procéder à l’évaluation, on effectue généralement un calcul du RMSD (Root Mean square Deviation) entre la conformation native et celle obtenue par prédiction. Cette opération mesure simplement les différences de position entre les atomes de Cα représentatifs du squelette des structures comparées. Bien entendu, il serait tout à fait envisageable d’effectuer un calcul RMSD sur chaque atome constituant la protéine. Notons également qu’un calcul RMSD retourne le plus petit résultat possible. En effet, le calcul s’effectue sur le meilleur alignement des deux structures comparées, ce qui a pour effet d’en minimiser le résultat. En notant par (xin , yin , zin ) les coordonnées de l’atome Cα du résidu i de la structure native et par (xip , yip , zip ) les coordonnées de l’atome Cα du résidu n de la structure prédite, le calcul du RMSD s’effectue en calculant : v t  r  1 X   i i 2 i i 2 i i 2  min  (xn − xp ) + (yn − yp ) + (zn − zp )   r  i=1 où r est le nombres de résidus de la protéine. Une seconde méthode d’évaluation serait de comparer les résultats énergétiques obtenus entre diverses méthodes de prédiction. Malheureusement, l’imperfection des multiples fonctions énergétiques existantes nous oblige à regarder ces résultats avec circonspection. Il n’est pas rare que l’évaluation énergétique d’une structure native connue offre un score supérieur à celle d’une prédiction ! D’autres comparaisons peuvent avoir lieu entre la protéine à l’état natif et la protéine prédite. Citons par exemple l’étude du pourcentage de contacts natifs,

CHAPITRE 2. PRÉSENTATION DES OUTILS ET MÉTHODES EXISTANTS 41 c’est-à-dire de contact inter-résidus semblables entre la structure native connue et la structure prédite. Un autre facteur fréquemment évalué est celui du temps requis par une prédiction. Au vu de la grande diversité de machines et de méthodes utilisées dans ce domaine, il va de soi qu’une unité d’évaluation temporelle simple serait un non sens. Il est donc fréquent d’utiliser comme source de prédiction une base de données contenant des protéines-type présentant des propriétés particulières (benchmarks). Ainsi, en spécifiant les caractéristiques de la machine utilisée, on peut se faire une idée de la qualité de la méthode de prédiction utilisée.

2.5 Récapitulatif La prédiction de la structure native d’une protéine a un intérêt scientifique important. Outre le gain de temps et d’argent que procurent les méthodes in silico, une multitude d’applications dérivées sont accessibles dès lors que la struture native d’une protéine est connue. Pour effectuer une prédiction de la structure native d’une protéine, certains outils sont indispensables, nous venons d’en faire le tour. Il est important de remarquer que quelles que soient les méthodes utilisées, la situation initiale de la recherche a une grande importance. De plus, la qualité des résultats obtenus est souvent suffisante pour des problèmes basiques mais ce n’est pas toujours le cas lors de l’étude de problèmes plus importants. La recherche du minimum global est un problème difficile pour lequel il n’existe actuellement aucun algorithme efficace garantissant un résultat optimal. Dans le chapitre suivant, nous présenterons la totalité des choix effectués pour le développement de notre logiciel de prédiction de structure native de protéine. Cela fait, nous nous lancerons dans la découverte de son fonctionnement pour ensuite en révéler les résultats.

Chapitre 3 Nos outils et méthodes Le problème de prédiction de la structure native d’une protéine est en fait un problème d’optimisation combinatoire particulier. En effet, il existe de nombreux problèmes de ce type. Chacun peut être défini par une paire (S, f ) où S représente l’ensemble fini de toutes les solutions existantes et f est une fonction de type f : S → R permettant d’évaluer la qualité d’une solution i particulière [22]. La valeur optimale de f est : f0 = min{ f (i) : i ∈ S} et l’ensemble des solutions optimales est : S0 = {i ∈ S : f (i) = f0 } La résolution du problème se fait donc en trouvant certaine(s) solution(s) i0 ∈ S0 . Afin de bien comprendre tout cela, transposons ces données mathématiques au cas qui nous intéresse. Pour cela, rappelons que nous nous basons sur l’hypothèse largement admise selon laquelle la structure native d’une protéine correspond à son minimum d’énergie libre. Dans notre cas, S est l’espace conformationnel accessible et f sera la fonction énergétique utilisée. Si nous appliquons ces formules à une protéine dont nous connaissons déjà la structure native (que nous noterons n) on aura que f (n) = f0 . D’autre part, on sait que n ∈ S et donc que n ∈ S0 car l’équation de S0 est bien vérifiée. Ces problèmes d’optimisation combinatoire sont difficiles à résoudre car pour la plupart, une exploration exhaustive des solutions possibles sera requise afin de déterminer la ou les meilleures solutions. Pour les problèmes de grande taille, une telle exploration n’est actuellement pas réalisable en un temps raisonnable. Il n’existe donc pas de méthode déterministe garantissant un résultat 42

CHAPITRE 3. NOS OUTILS ET MÉTHODES

43

optimum. Ainsi, des métaheuristiques proposent des stratégies générales d’exploration de l’espace qui peuvent être appliquées à un problème particulier. Dans notre cas, nous avons choisi de concevoir un algorithme de prédiction de la structure native de protéines à partir de la métaheuristique des colonies de fourmis que nous avons explorée dans le chapitre précédent. Ce chapitre-ci propose donc de découvrir les divers choix que nous avons effectués quant aux diverses exigences requises par la prédiction de structure native.

3.1 Modélisation de nos protéines Nous avons opté pour une représentation atomique simplifiée des acides aminés (cfr. figure 3.1). Ainsi, comme nous le voyions dans l’introduction, nous ne prenons pas en compte la totalité de la chaîne latérale. Celle-ci sera simplifiée au moyen des atomes Cβ et Cµ . Cette représentation offre l’avantage d’être assez réaliste tout en ne demandant pas une quantité de calcul trop importante.

F. 3.1 – Représentation simplifiée d’un acide aminé, attention le cas particulier de la Glycine verra ses atomes Cβ et Cµ se trouver au même endroit que l’atome Cα . Une protéine sera donc représentée par sa chaîne principale (son squelette), ainsi que par les atomes Cbeta et Cµ de chaque résidu. Par ailleurs, nos protéines évolueront dans l’espace (φ, ψ, ω), c’est-à-dire que les angles de rotation autour des liaisons chimiques N − Cα et Cα − C ainsi que la liaison C − N entre un résidu et son voisin, prendront respectivement les valeurs φ, ψ et ω. Rappelonsnous également que l’angle ω résulte d’un lien peptidique particulier ayant un caractère de double liaison, ce qui limite ses valeurs à 0◦ ou 180◦ . Le lecteur se


44

posant des questions sur les valeurs initiales des angles φ et ψ vera sa curiosité assouvie en page 55. Il nous faut encore préciser que les distances de liaison inter-atomique auront des valeurs fixes variant en fonction des atomes liés. Cela n’entraîne aucune surcharge de calcul mais nous permet d’approcher encore un peu mieux la réalité. Enfin, dans notre cas de représentation simplifiée des acides aminés, les angles de valences (cfr. figure 3.2 ) seront également fixés pour chaque atome.

F. 3.2 – Exemple d’angle de valence : en pointillé, l’angle de valence associé à l’atome O. On dira qu’une prédiction de protéine sera valable si elle ne contient pas de clash stériques. Cette appellation dénote en fait une collision entre deux atomes. En effet, il est tout a fait possible de prédire des protéines présentant une structure tridimensionnelle irréalisable dans la nature. Il devra donc être possible d’évaluer la présence de telles collisions. Dans notre modèle, nous veillerons à l’absence de collisions entre chaînes latérales et donc uniquement entre les atomes Cµ . Pour effectuer ces vérifications, un rayon atomique est fixé pour l’atome Cµ de chaque résidu, on l’appelle le rayon de van der Waals [17]. En notant par rVDW(Cµ ) le rayon de Vanderwaals de l’atome Cµ , on dira qu’il y a collision si rVDW(C1µ ) + rVDW(C2µ ) 2

· a > distanceSpatiale((C1µ ), (C2µ ))

où distanceSpatiale calcule la distance réelle séparant deux atomes et a représente un facteur d’aération permettant donc d’éviter des prédictions trop compactes. Il faudra donc évaluer l’absence de collision entre tous les atomes Cµ de la protéine prédite afin que celle-ci soit réalisable et puisse participer à la recherche


45

de la structure native.

3.2 Discrétisation de notre espace de recherche Comme notre protéine évolue dans l’espace (φ, ψ, ω), il va de soi que la discrétisation de l’espace conformationnel passera par une limitation des valeurs accessibles par chaque triplet angulaire. Au sein d’une protéine réelle, certaines valeurs de (φ, ψ, ω) ne sont pas accessibles aux résidus car elles correspondraient à des collisions atomiques. Il a été observé que l’ensemble des conformations possibles pour chaque résidu peut être réduit à sept domaines dans l’espace des triplets (φ, ψ, ω) [33] : A, C, B, P, G, E et O. La discrétisation de l’espace conformationnel s’opère donc en dégageant un certain nombre de représentants par domaine, calculés à partir de conformations natives réelles observées expérimentalement. Ces représentants peuvent être calculés par résidu car la conformation de chacun d’entre eux est sensiblement différente dans la réalité. Néanmois, nous avons initialement décidé d’utiliser des représentants généraux pouvant être appliqués sur n’importe quel résidu, cela facilitant les opérations à effectuer (nous verrons plus loins que finalement, il eut été aisé de prendre des représentants particuliers pour chaque résidu). Dans le cas ou seulement un représentant est disponible, il n’y aura aucun choix stochastique à effectuer. Par contre, dès que plusieurs représentants sont accessibles, il faut permettre une sélection proche de la réalité. C’est pourquoi en plus des trois données angulaires nécessaires, un représentant contient une quatrième donnée : son nombre d’occurences. En effet, diverses conformations sont observées dans la nature et certaines d’entre-elles sont plus fréquentes que d’autres. Recherchant un maximum de réalité lors de la conception de notre modèle, il nous a paru évident d’opérer un choix stochastique dépendant du nombre d’occurence de chaque représentant. Cela se fait simplement en associant à chaque triplet angulaire de chaque domaine une probabilité d’apparition. Cette probabilité intra-domaine se calcule simplement en divisant le nombre d’occurences du représentant considéré par la somme du nombre d’occurences de chaque représentant. En effectuant ce calcul sur chaque domaine, on obtient donc une discrétisation de l’espace conformationnel accessible. Il est important de remarquer que cette discrétisation s’opère par un système


46

de coordonnées interne à la protéine. Ainsi, la représentation tridimensionnelle d’une prédiction s’obtiendra en traduisant les coordonnées (φ, ψ, ω) en coordonnées (x, y, z). Pour effectuer cette traduction, les données propres de chaque atome seront nécessaires (angle de valence et distance inter-atomique).

3.3 Exploration de l’espace conformationnel Comme nous l’écrivions en début de chapitre, nous avons opté pour une méthode de recherche dérivée de la métaheuristique des colonies de fourmis. Le chapitre précédent nous proposait un aperçu général du fonctionnement de ces colonies, examinons donc ensemble les spécificités de fonctionnement de notre colonie. L’univers d’évolution de nos fourmis artificielles correspond concrètement à l’espace conformationnel accessible par la protéine, celui-ci étant discrétisé en sept domaines. La tâche de nos fourmis consistera donc à trouver une suite de domaines offrant une structure tridimensionnelle minimisant la fonction énergétique que nous découvrirons plus loin. Afin de réaliser au mieux sa tâche, une fourmi opérera une suite de choix probabilistes basés sur la quantité locale de phéromone τ ainsi que sur les valeurs heuristiques disponibles η. Une question se pose ici : à quoi faire correspondre cette valeur η ? L’idée est de trouver un indicateur permettant d’aider la fourmi dans son choix. Il faudrait donc disposer d’une valeur indiquant un domaine de conformation à priori favorable pour un résidu donné. Heureusement, depuis quelques années, il existe à l’ULB certains projets de recherche et de prédiction correspondant exactement à nos besoins. Ces projets se nomment Prelude et Fugue [14, 34]. Ainsi, il est possible d’obtenir, pour une séquence donnée, une séquence correspondante de domaines favorables ainsi qu’un facteur exprimant la valeur de cette correspondance. Ce facteur prend des valeurs entières comprises entre 0 et 9 et lorsqu’il vaut 0, c’est qu’il n’y a pas de prédiction pour le résidu concerné (voir exemple ci-dessous). Il est important de signaler que ces prédictions se basent exclusivement sur les interactions locales le long de la séquence et ne prennent pas en compte les interactions tertiaires ! La séquence de domaines fournie par ces programmes ne constitue donc généralement pas une solution convenable pour la structure tertiaire complète de la protéine mais elle est utile car elle permet de connaître les préférences individuelles de petites portions de la séquence.

CHAPITRE 3. NOS OUTILS ET MÉTHODES Séquence: Domaine favorable: Facteur de correspondance:

47

QQRLIFAGKQLEDGRTLSDYNIQKESTLHLVLRLRGG CPBBBBGGBPPBACPBAAxCGPBPGPBBBBBBAxxxx 1145789996211222220467776412333320000

Grâce à ces projets, nous disposons donc d’un biais inter-domaines en plus du biais intra-domaine (dû à la probabilité d’apparition de chaque représentant d’un domaine). Fort de cet avantage, nous avons donc décidé d’utiliser comme valeur heuristique pour un domaine, le facteur proposé par Prelude et Fugue, quand il existe. Cela nous permet de procéder au calcul que nous proposait l’équation du chapitre précédent : pi j = P

[τi j (t)]α · [ηi j (t)]β k∈V(i) [τik (t)]

α

· [ηik (t)]β

où pi j représente donc la probabilité d’associer le domaine j au résidu i. Il sera donc très important de paramétrer au mieux les facteurs α et β. Maintenant que nous avons une stratégie de décision, il nous faut fixer les caractéristiques du dépôt phéromonal de nos fourmis. Cela signifie qu’il nous faut décider quand déposer des phéromones et en quelle quantité. Avant d’expliquer nos choix, quelques définitions sont nécessaires. Appelons itération une recherche complète effectuée par une fourmi, c’est-à-dire qu’après une itération, une fourmi aura construit et évalué une prédiction complète et mourra. De plus, nommons cycle un nombre m d’itérations où m représente le nombre total de fourmis de la colonie. Ainsi, après un cycle, toutes les fourmis auront procédé à l’élaboration et l’évaluation d’une prédiction complète. Le nombre de cycles sera une variable à définir lors de l’exécution de l’algorithme. Enfin, en se souvenant que nous travaillons dans un univers où le temps est discrétisé, nous pouvons également dire qu’un cycle se passera en une unité de temps. Cela signifie qu’un cycle commençant en (t) se terminera en (t + 1). Revenons maintenant au dépôt de phéromones. Nous avons décidé d’effectuer un dépôt offline dans le but de permettre à chaque fourmi de disposer de la même base phéromonale lors d’un cycle de recherche. En effet, le dépôt offline s’effectue après un cycle entier, ce qui permet de ne pas modifier la situation initiale du cycle de recherche. Nous utilisons donc les propriétés de la discrétisation du temps, nous permettant d’effectuer une simulation d’activité simultanée de la colonie sans aucune interférence inter-fourmi. Bien entendu, certaines structures de données seront nécessaires à l’application de ce choix, nous aurons l’occasion de les découvrir au chapitre suivant.


48

Outre ce dépôt phéromonal offline, nous avons également pensé à exécuter un dépôt phéromonal élitiste, c’est-à-dire de déposer des phéromones supplémentaires pour les x meilleures conformations protéiques. Une fois de plus, x sera un paramètre à définir, de même que le facteur d’élitisme spécifiant la quantité de phéromone additionelle à déposer. Formellement, nous obtenons donc que : τi j (t + 1) = ρτi j (t) +

m X

∆τkij (t, t

+ 1) +

e X

k=1

∆τli j (t, t + 1)

l=1

où e est le nombre de fourmis élitistes et m le nombre total de fourmis. Nous n’en avions pas encore parlé mais il est évident qu’un taux de persistance ρ sera utilisé afin de permettre le renouvellement phéromonal et donc l’exporation de l’espace conformationnel. En résumé, notre prédiction se fera au moyen d’une colonie de fourmis contenant un certain pourcentage de fourmis élitistes et procédant d’une part à un dépôt phéromonal offline et d’autre part à un dépôt phéromonal élitiste. Les particularités d’implémentations seront explorées au prochain chapitre.

3.4 Fonction énergétique utilisée Dans le cadre de ce mémoire, nous utiliserons une classe de fonctions énergétiques appelée potentiels statistiques. Ces potentiels sont dérivés de bases de données de protéines dont les structures natives sont connues [35]. Dans notre étude, nous avons décidé d’utiliser deux types de potentiels : les potentiels de torsion entre petits éléments de structure et les potentiels de distance entre petits éléments de séquence.

3.4.1 Les potentiels de torsion Les potentiels de torsion décrivent exclusivement les interactions locales le long de la chaîne principale et sont donc lié à la prédiction de structures secondaires telles que des hélices-α. Chaque potentiel porte un nom différent où s signifie résidu de la séquence et t signifie domaine de torsion. Ces potentiels mesurent l’adéquation des données qui leurs sont fournies. Si cette adéquation est bonne, l’énergie correspondante sera basse. Par exemple, le potentiel st0 mesure l’adéquation entre la nature d’un résidu donné et le domaine prédit pour ce résidu ou pour un résidu voisin. Prenons le


49

cas de l’acide aminé nommé alanine. S’il se trouve en position 5 dans la séquence et qu’en position 7, le domaine prédit est P, alors l’énergie fournie par le potentiel st0 vaudra 0.114887. Par contre si le domaine prédit avait été A, l’énergie aurait été de -0.111486. Il faudra donc effectuer ces calculs le long de toute la séquence de la protéine. Toutefois, au vu du caractère local des potentiels de torsion, la fenêtre d’exploration autour du résidu de base sera de 8 résidus dans la séquence. Dans le cas du potentiel st0, on évaluera donc l’adéquation entre un résidu en position i de la séquence et une prédiction de domaine en position j de la séquence, avec −8 ≤ ( j − i) ≤ 8. Les potentiels de torsion que nous utiliserons sont les suivants : – st0, qui mesure l’adéquation entre le résidu positionné en i et le domaine prédit en position j ; – stt1, qui mesure l’adéquation entre le résidu positionné en i, le domaine prédit en position j et le domaine prédit en position k ; – tss1, qui mesure l’adéquation entre le domaine prédit en position i, le résidu positionné en j et le résidu positionné en k ; – tt0, qui mesure l’adéquation entre les domaines prédits en position i et j ; – ttt1, qui mesure l’adéquation entre les domaines prédits en position i, j et k. Le lecteur attentif aura remarqué une redondance de calcul dans les potentiels tt0 et ttt1. En effet, pour le potentiel tt0 (mais cela est transposable pour le potentiel ttt1), qu’on aie les domaines P et B respectivement en position 2 et 4 ou l’inverse revient exactement au même calcul ! Ainsi, pour ces deux potentiels, le calcul ne s’effectuera que sur une fenêtre positive, c’est-à-dire pour les 8 résidus voisins de droite dans le sens de lecture de la séquence de la protéine. Il sera encore temps de répartir ce potentiel énergétique intelligemment, nous y reviendrons au chapitre suivant .

3.4.2 Les potentiels de distance Nous le découvrions dans l’introduction, la création de la structure tertiaire d’une protéine est principalement due aux interactions entre résidus distants dans la séquence mais proches dans l’espace. Les potentiels de distance sont donc un moyen de mesurer la contribution de telles interactions. Dans ce cas-ci, il est tout à fait possible que des résidus espacés par plus de 8 autres résidus dans la séquence se retrouvent rapprochés dans l’espace. Nous


50

effectuerons donc ce calcul pour chaque résidu par rapport à toute la séquence à l’exception de ses deux voisins directs dans la séquence, qui lui seront forcément très proche dans l’espace. Par contre, si les résidus considérés sont distants dans l’espace de plus de 8Å, on considèrera qu’il n’y a tout simplement aucune interaction. L’énergie calculée dépend donc de la nature des deux résidus ainsi que de leur distance spatiale, mais aussi de leurs positions relatives dans la séquence. Par exemple, si nous avons une alanine en position 5 et une proline en position 8, l’énergie ne sera pas la même que pour une alanine en position 5 et une proline en position 9 présentant la même distance spatiale. Remarquons tout de même que lorsque les deux résidus considérés sont éloignés de plus de 8 résidus dans la séquence, le potentiel calculera une valeur énergétique dépendant uniquement du type de résidu considéré et de la distance spatiale. La distance dans la séquence ne sera donc plus prise en compte. Par exemple, admettons que le résidus 5 soit toujours une alanine et que les résidus 24 et 65 soient des leucine. Si l’on considère les résidus 8 et 24 spatialement distants de 3Å et les résidus 8 et 65 également distants de 3Å, le potentiel énergétique calculé sera le même. Nous utiliserons deux potentiels de distance : ds0 et ds1. Le potentiel ds0 mesure l’adéquation entre la présence d’un résidu particulier en position i de la séquence et un résidu quelconque en position j présentant une certaine distance spatiale avec le résidu i. D’un autre côté, le potentiel ds1 effectue le même calcul mais avec un résidu défini en position j. Par exemple, prenons le cas où le résidu 5 est une alanine et le résidu 8 une proline, ces deux résidus étant séparé par 4Å. Le potentiel ds1 calculerait une valeur énergétique de 0.031414 et le potentiel ds0 une valeur énergétique de 0.397226. Si par contre le résidu 8 était une valine, la valeur prise par le potentiel ds1 serait différente, mais pas celle prise par le potentiel ds0. A nouveau, comme pour les potentiels tt0 et ttt1 remarquons que le potentiel ds1 offre un certain niveau de redondance, ce qui permet de réduire les calculs.

3.4.3 Utilisation combinée des potentiels Notre espoir est que l’utilisation combinée de ces deux types de potentiels permette la prédiction de structures secondaires tout en prenant en compte les interactions tertiaires. Ainsi, la protéine prédite devrait avoir un degré de similitude plus important avec la structure tidimensionnelle réelle d’une protéine


51

que dans le cas ou seuls les potentiels de torsion seraient pris en compte.

3.5 Evaluation de nos solutions Afin d’évaluer nos solutions, nous procéderons au calcul du RMSD (Root Mean Square Deviation) sur les atomes Cα de nos protéines prédites. Pour rappel, ce calcul permet de mesurer les différences de position entre deux structures données. Nous procèderons donc à un jeu de test sur un fragment de protéine dont la structure native est déjà connue. Nous découvrirons également dans quelle mesure notre algorithme trouve des conformations d’énergie libre minimale et profiterons de l’occasion pour faire une comparaison entre les énergies obtenues et les RMSD correspondants.

3.6 Récapitulatif Nous venons de découvrir l’ensemble des outils que nous utiliserons pour procéder à la prédiction de la structure native de protéines. Dans le chapitre suivant, nous allons découvrir les spécificités de l’algorithme que nous avons développé avant d’en découvrir les tests et les résultats. Ce chapitre présentera les diverses particularités d’implémentation des outils utilisés.

Chapitre 4 Notre algorithme : AC-ProPre Dans ce chapitre, nous allons passer en revue les spécificités d’implémentation des outils découverts au chapitre précédent. Notre algorithme, AC-ProPre (Ant Colony for Protein Prediction), ayant été réalisé en C++, nous structurerons ce chapitre à l’image du code source, c’est-à-dire classe par classe. Nous commencerons donc par explorer l’implémentation de la discrétisation de l’espace conformationnel pour ensuite découvrir la manière dont les potentiels énergétiques sont codés et utilisés. Enfin, nous nous pencherons sur la représentation des protéines ainsi que sur les divers calculs nécessaires à l’élaboration de la prédiction. Pour finir nous découvrirons évidemment la manière de fonctionner de nos fourmis artificielles.

4.1 L’espace conformationnel Le code source de cette section se trouve dans les fichiers domaine.h et domaine.cc. Comme nous l’avons maintes fois signalé, nous travaillons dans l’espace (φ, ψ, ω) pouvant être représenté en sept domaines : A, C, B, P, G, E et O. Il y aura donc une instance de la classe domaine pour chaque domaine utilisé. Un domaine contient des représentants. Comme nous l’expliquions au chapitre précédent, nous avons choisi de travailler avec des représentants communs, c’est-à-dire appliquables à chaque résidu. Néanmoins, une amélioration relativement simple à implémenter serait de considérer des représentants propres à la nature du résidu, cela pourra être fait dans une version ultérieure de notre algorithme. Les représentants proviennent de fichiers source nommés cah_rama_141_x 52

CHAPITRE 4. NOTRE ALGORITHME : AC-PROPRE

53

où x représente le nombre de représentants effectifs par domaine. Il est donc possible de travailler avec un nombre de représentants variable. Il est évident qu’un nombre de représentants élevé permettra une prédiction plus proche de la réalité. Chaque représentant est caractérisé par un triplet angulaire et un nombre d’occurences. C’est au moyen de ce nombre d’occurence que nous calculons la probabilité de sélection d’un représentant particulier. Pour calculer cette probabilité, il nous suffit de diviser le nombre d’occurences de chaque représentant par le nombre d’occurences totale d’un domaine. Ainsi, nous obtenons un modèle probabiliste cohérent où la somme des probabilités d’apparition de chaque représentant vaut 1. Afin de permettre un choix aléatoire rapide pour la sélection de représentants, nous construisons un vecteur d’intervalle de probabilité proportionnel à la probabilité d’apparition de chaque représentant. Ainsi, à chaque fois qu’un représentant devra être sélectionné, il suffira de tirer aléatoirement un nombre entre 0 et 1 et de choisir le représentant dont l’intervalle contient ce nombre aléatoire. Nous voilà donc avec un espace conformationnel discrétisé. Pour obtenir un représentant, il nous suffira d’appeler la fonction getRepresentant().

4.2 Les potentiels énergétiques Chaque potentiel énergétique utilisé est fourni sur fichier texte. Sachant qu’il existe vingt types d’acides aminés et sept domaines, ce n’est pas une surprise d’apprendre que le fichier potentiel le plus lourd (tss1) fait 6,1 Mo. Au vu du nombre d’accès nécessaires au calcul de l’énergie libre d’une conformation, il est évident qu’il faut trouver un moyen de stockage en mémoire vive pour ces données vitales. Malheureusement, les potentiels sont tels qu’un gaspillage significatif de la mémoire vive aurait lieu si nous utilisions simplement des matrices multidimensionnelles pour le stockage de leurs données. De plus, chaque potentiel a des caractéristiques propres qui empêchent la création d’une méthode générale automatique de stockage et d’accès. Toutefois, un concept commun peut-être dégagé : l’idée est de stocker les valeurs énergétiques d’un potentiel dans un simple vecteur de données et de développer parallèllement une méthode d’accès à ce vecteur basée sur les paramètres nécessaires au calcul du potentiel énergétique.


54

Nous avons donc développé, pour chaque potentiel, un calcul de l’espace nécessaire pour le stockage des énergies ainsi que la méthode d’accès à cet espace. Le problème était simple : trouver l’index dans le vecteur correspondant aux paramètres de calcul du potentiel afin de retourner la bonne valeur énergétique. Il était donc nécessaire d’attribuer à chaque domaine ainsi qu’à chaque acide aminé un numéro d’identification unique. Ainsi, afin d’éviter toute erreur de manipulation ou de conversion, les fonctions communes ddTOint() et aaTOint() sont un point de passage obligé convertissant un domaine (respectivement un acide aminé) en nombre entier. Il est évident que ces nombres sont une suite croissante commencant à 0. Dans le cas des potentiels de distances, une fonction similaire a été prévue afin de convertir la distances spatiale entre deux résidus. Par ailleurs, afin de généraliser l’utilisation de notre algorithme, nous avons tenu compte de l’éventuelle modification de la taille de la fenêtre. Pour rappel, cette fenêtre délimite l’espace possible entre deux positions de la séquence pour le calcul d’un potentiel. Dans notre cas, comme présenté au chapitre précédent, cette fenêtre vaut 8. Finalement, nous avons obtenu des vecteurs de taille optimale par rapport aux spécificités des fonctions potentiels et ce pour un coût de calcul supplémentaire négligeable. En effet, si nous avions travaillé avec des matrices multidimensionnelles, le calcul d’index automatique aurait été comparable, mais nous aurions gaspillé beaucoup d’espace mémoire. Nous pouvons donc conclure en nous réjouissant de disposer de fonctions de calcul de potentiel simples. Par exemple, un simple appel à la fonction getDs0() avec les paramètres nécessaires nous retournera le potentiel énergétique associé. Bien entendu, cet exemple se transpose à chaque potentiel utilisé.

4.3 La structure et les calculs effectués sur les protéines Maintenant que nous savons comment obtenir un représentant de domaine ainsi qu’un potentiel énergétique, nous pouvons nous attaquer à la structure de la protéine proprement dite.


55

4.3.1 La structure protéique tridimensionnelle Avant toute chose, il est important de faire une petite précision en ce qui concerne les angles de torsion φ, ψ et ω. Comme nous le savons, ce sont eux qui détermineront la conformation d’un résidu. Plus précisément, la conformation du résidu sera calculée sur base de ces trois angles ! En effet, au repos, ces angles de torsion ne sont pas nuls. C’est pourquoi nous avons besoin des données primaires sur ces angles afin d’en calculer la valeur réelle. Par exemple, considérons l’angle φ. La valeur prise par cet angle influe directement sur la position des atomes C, Cβ et Cµ . Mais concrètement, l’angle φ proposé par les représentants nous donne la valeur de l’angle de rotation pour l’atome C. Il nous faudra donc ensuite mettre à jour les donnée des atomes Cβ et Cµ . Cela se fait assez facilement grâce aux formules suivantes où TX0 représente l’angle de torsion initial de l’atome X du résidu et TX représente son angle réel : TC − TC0 = TCβ − TC0 β et donc :

(TC − TC0 ) + TC0 β = TCβ

A partir de l’angle φ, qui n’est autre que l’angle de torsion réel de C, nous obtenons donc un écart qu’il nous suffira d’ajouter à l’angle de torsion initial de Cβ afin d’obtenir son angle de torsion réel . En ce qui concerne l’atome Cµ , le calcul est comparable à cette différence près que de par son statut particulier, le Cµ possède un angle de torsion différent pour chaque résidu. C’est d’ailleurs le moment de rappeler l’exception du résidu glycine qui n’a pas de chaîne latérale conséquente et pour lequel l’angle de rotation du Cµ vaut forcément 0 ! En conclusion, certaines données extérieures seront nécessaires au calcul de la structure tridimensionnelle de la protéine : – les angles de valence et de torsion de tous les atomes de la chaîne principale ; – les distances inter-atomiques de chaque atome de la caîne principale ; – les données particulières concernant l’atome Cµ de chaque résidu. Nous aurons également besoin du fameux rayon de van der Waals des atomes Cµ afin de détecter la présence de collisions dans nos prédictions. Ces données se trouvent dans les fichiers levitt.in pour le Cµ et dans donneeAtome.h pour les autres atomes.


56

Maintenant que ces subtilités sont acquises, nous pouvons examiner la représentation de la structure d’une protéine utilisée par notre programme. Notre idée générale était de permettre un accès direct et rapide aux divers éléments nécessaires à l’obtention d’une structure protéique particulière. Nous avons donc stocké toutes ces informations dans des vecteurs similaires, basés sur le nombre de résidus composant la protéine. Ainsi, avec un numéro d’index unique par résidu, nous avons la possiblité d’accèder directement à l’information nous intéressant dans chaque structure de données concernant ce résidu. Ces vecteurs nous permettent donc de retrouver instantanément et pour chaque résidu, la nature du résidu, le domaine de conformation dans lequel il se trouve, le représentant choisis dans ce domaine et les angles φ, ψ et ω proposé par ce représentant. De plus, lorsque les calculs énergétiques et tridimensionnels auront eu lieu, nous pourrons accéder à la valeur énergétique de la protéine prédite ainsi qu’à sa structure tridimensionelle en coordonnées x, y et z. La section suivante donnes quelques explications à propos de ces calculs.

4.3.2 Les méthodes de calcul Afin de pouvoir procéder à la comparaison des diverses prédictions réalisées par nos fourmis artificielles, il sera nécessaire de connaître la structure tridimensionnelle prédite pour la protéine et l’énergie libre qui lui est associée. Il y aura donc principalement deux calculs à effectuer : 1. la construction de la structure tridimensionnelle prédite pour la protéine ; 2. le calcul de l’énergie libre de cette conformation. D’autres calculs tels que la distance spatiale séparant deux points ou la présence de collisions inter-atomiques dans la séquence seront à effectuer mais de par leur simplicité, nous ne nous y étendrons pas. La construction de la structure tridimensionnelle Le but du calcul est très simple : construire la structure tridimensionnelle de la protéine à partir des angles φ, ψ et ω de chacun de ses résidus. Ce calcul s’effectue au moyen de matrices de rotation permettant de calculer les coordonnées spatiales d’un point en fonction : – des coordonnées du point le précédant ; – de la distance séparant ces deux points ;


57

– de deux angles caractérisant les mouvements à effectuer (dans notre cas, il s’agit des angles de torsion et de valence de l’atome considéré que nous avons défini respectivement aux pages 11 et 44). Dans notre cas, nous fixerons toujours les coordonnées spatiales du premier atome du premier résidu à (0,0,0). De plus, nous supposerons que la structure évolue sur l’axe des x pour le calcul du second point. Cette situation initiale étant définie, la suite des calculs pourra se passer itérativement. Il est intéressant de remarquer qu’une fois les calculs effectués pour un atome, il est possible de dégager une matrice contenant toutes les informations géométriques nécessaires concernant cet atome. En sauvegardant cette matrice, il sera donc possible de recommencer le calcul à partir d’un certain point de la structure, ce qui évitera de devoir recalculer toutes les coordonnées de tous les atomes dans le cas où un résidu viendrait à changer de conformation. Nous effectuerons donc une sauvegarde de cette matrice exclusivement pour les atomes C de chaque résidu. En effet, lorsqu’un résidu change de conformation, c’est chacun de ses angles de torsion qui prend une nouvelle valeur. Il est donc inutile de conserver les matrices internes d’un résidu. Nous avons donc ajouté un vecteur supplémentaire dans la représentation de nos protéines, permettant de recommencer une construction de structure tridimentionnelle au départ d’un point particulier. Cela sera particulièrement intéressant dans le cas où une collision interatomique sera découverte. En effet, lorsqu’une nouvelle conformation réaliste sera trouvée pour ce résidu, la reconstruction tridimentionnelle pourra recommencer non pas au début mais bien sur le résidu précédant directement celui nouvellement modifié ! Le gain de temps sera donc significatif, d’autant plus que l’espace mémoire supplémentaire requis par cette opération est minime ! Le calcul de l’énergie libre d’une conformation prédite Ce calcul s’effectue sur base des valeurs retournées par les potentiels énergétiques. Il nous faudra donc faire appel aux bonnes fonctions avec les bons paramètres, ce qui est aisé vu la structure de stockage des données utilisée. Par contre, un calcul mons évident est celui de la répartition des énergies. En effet, lors du chapitre précédent, nous avons soulevé le point de la répartition des énergies potentielles calculées lorsqu’un potentiel présente des redondances. Pour reprendre l’exemple de la page 49 dans le cas du potentiel tt0 (mais cela est transposable pour le potentiel ttt1), qu’on aie les domaines P et B respectivement


58

en position 2 et 4 ou l’inverse revient exactement au même calcul ! Dés lors, il suffira en fait de reporter ce potentiel sur chaque résidus concerné, dans notre exemple les résius 2 et 4. Ce type de redondance se présente concrètement pour les potentiels tt0, ttt1 et ds1. Toutefois, un autre type de répartition, quelque peu plus subtil, doit également être effectué ! Pour bien comprendre ce nouveau problème, quelques précisions sont nécessaires. Lors de la prédiction de structure native de protéine, la séquence source est bien évidemment une donnée immuable. Par contre, la prédiction de la conformation adoptée par chaque résidu changera tout au long de la recherche. De plus, pour définir si une conformation est intéressante ou non, nous procédons à l’évaluation de celle-ci au travers des potentiels de torsion, entre-autres. Ces potentiels expriment dans quelle mesure il est intéressant d’avoir une prédiction de domaine particulière pour un résidu particulier. Ainsi, lorsque nous procédons par exemple au calcul du potentiel st0 qui mesure l’adéquation entre le résidu positionné en i et le domaine prédit pour un résidu en position j, une question se pose : à quel résidu attribuer l’énergie calculée ? C’est de cette question que provient le nouveau type de répartition que nous venons de soulever. Toujours dans le cas du potentiel st0, l’énergie sera attribuée au résidu de position j car c’est bien le domaine prédit qui doit être évalué et non la nature du résidu en position i de la séquence, qui de toute façon ne changera jamais. Le même raisonnement est applicable dans le cas des potentiels stt1 et tss1. Pour rappel, le potentiel stt1 mesure l’adéquation entre le résidu positionné en i, le domaine prédit en position j et le domaine prédit en position k. L’énergie calculée sera répartie sur les deux résidus j et k qui présente une prédiction particulière. En ce qui concerne le potentiel tss1, qui mesure l’adéquation entre le domaine prédit en position i, le résidu positionné en j et le résidu positionné en k, l’entièreté de l’énergie calculée sera reportée sur le résidu positionné en i dans la séquence. Une dernière remarque concerne la valeur des énergies retournées par les différents potentiels. En effet, certains retourneront facilement des valeurs fortement négative lors même que d’autres seront plus difficiles à minimiser. Nous avons donc décidé de sous pondérer les valeurs retournées par les potentiels tt0 et ttt1. En effet, ces derniers sont facilement minimisables et la valeur énergétique qu’ils calculent pourrait être telle que l’impact des autres calculs de potentiels pourrait être négligeable. Nous nous retrouverions donc dans le cas d’une prédiction où seules les interactions locales seraient prise en compte, ce


59

que nous souhaitons éviter.

4.4 La colonie de fourmis Nous voici arrivé à la classe maîtresse de notre algorithme : la colonie de fourmis. C’est en effet ici que les différentes pièces du puzzle s’assemblent afin de permettre la prédiction de la structure native de la protéine étudiée. Pour commencer, examinons les paramètres requis par la colonie : – le nombre d’itération maximale, critère de terminaison principal ; – le nombre de fourmis qui participeront à la recherche ; – le nombre de fourmis d’élite et le facteur d’élitisme indiquant le degré d’élitisme souhaité ; – les paramètres de calcul probabilistes fondamentaux α, β et ρ ; – enfin, le nombre maximum d’itérations acceptées sans modification substancielle de la meilleure solution et une valeur quantitative représentant cette amélioration minimale. Ces deux paramètres constituent un critère de terminaison anticipée. Le lecteur s’en sera certainement rendu compte, le nombre de paramètres attendu est assez important. Par ailleurs, la question se pose de savoir quelle valeur donner à chacun de ces paramètres. Nous voici donc confronté au problème de la paramétrisation dont nous parlions au chapitre deux. Concrètement, il n’existe aucune configuration prédéfine générale, il faut donc procéder à une multitudes de tests pour espérer dégager un comportement particulier de l’algorithme en fonction du problème à solutionner. Nous reviendrons sur ces tests au chapitre suivant. Examinons maintenant les structures de données internes. La mémoire de la colonie, c’est à dire l’ensemble des phéromones déposées par toutes les fourmis de la colonie, se trouve stockée sous forme de matrice tridimensionnelle que nous appellerons matrice phéromonale globale. Ainsi, la mise à jour et l’évaporation phéromonale peuvent s’effectuer assez simplement. Par ailleurs, l’intégralité de l’espace matriciel est potentiellement utilisé : chaque résidu peut adopter une conformation propre à un des 7 domaines et chaque domaine aura le même nombre fixe de représentant. Nous obtenons donc une matrice de r ∗ 7 ∗ R où r représente le nombre de résidus de la séquence et R donne le nombre de représentant par domaine. Parallèlement à cette matrice, nous avons stocké les valeurs heuristiques dans un vecteur. Pour rappel, ces valeurs heuristiques permettront


60

aux fourmis de s’orienter préférablement vers un domaine en particulier. Dès lors, au départ de ces deux sources de données, il nous est possible de calculer la matrice de probabilité générale qui sera utilisée par les fourmi. En effet, nous effectuerons ce calcul à chaque début de cycle (pour rappel, un cycle est constitué de m itérations où m représente le nombre total de fourmi de la colonie et où une itération consiste en la prédiction et l’évaluation d’une conformation complète par une fourmi). De plus, comme nous le disions précédament, nous avons décidé de procéder à une mise à jour phéromonale offline. Nous auront donc besoin d’une matrice de stockage temporaire afin de simuler la simultanéité de recherche des fourmis. Ainsi, à la fin d’un cycle, nous disposerons d’une matrice contenant la totalité des phéromones déposées pas nos fourmis. Par ailleurs, afin de procéder à la mise à jour phéromonale élitiste, les meilleurs solutions seront conservées jusqu’à la fin d’un cycle. Il sera ensuite temps d’effectuer un dépôt phéromonal additionnel dans la matrice phéromonale globale. Par la suite, seule la meilleure conformation prédite depuis la genèse de la colonie sera conservée. Il y aura donc de nouvelles conformations d’élite différentes à chaque fin de cycle. A la fin de l’algorithme, la colonie renverra tout simplement la meilleure conformation protéique prédite. Généralement, cette conformation sera aussi écrite sur fichier PDB [8] afin de pouvoir procéder ensuite à une visualisation de la protéine prédite. On pourra également, dans le cas où la structure native de la protéine est connue, procéder à un calcul RMSD de la prédiction avec la structure native réelle. Les différentes étapes de notre algorithme sont schématisées à la figure 4.1.

4.5 Récapitulatif Ce chapitre nous a permis de passer en revue les diverses spécificités de notre algorithme. Ainsi, au fil des différents chapitres de cet ouvrage, nous avons découvert toujours plus en profondeur l’application et le développement de la métaheuristique des colonies de fourmis au problème de la prédiction ab initio de la structure native d’une protéine. Le chapitre suivant propose donc de découvrir quelques réultats obtenus par AC-ProPre. Nous pourons ensuite tirer les conclusions de notre travail et développer quelques perspectives de développement futur.


F. 4.1 – AC-ProPre : déroulement schématique.

61

Chapitre 5 Tests et résultats Afin d’évaluer les performances de notre algorithme, nous avons effectué une batterie de tests sur les fragment 1 à 35 de la protéine 1UBQ, l’ubiquitine. La structure native de cette protéine étant connue, il nous est possible de calculer le RMSD entre notre prédiction et la conformation réelle des fragments 1 à 35 de l’ubiquitine. Afin d’obtenir un maximum d’informations intéressantes, nous avons évidemment exécuté notre algorithme dans de multiples configurations paramétriques. Nous allons donc dans un premier temps découvrir ces configurations pour ensuite en étudier les résultats.

5.1 Les configurations d’exécution de notre algorithme Nous avons découvert au chapitre précédent que notre algorithme possède beaucoup de paramètres pouvant chacun prendre de nombreuses valeurs différentes. Il a donc été nécessaire de faire quelques choix afin de procéder à nos tests. Ainsi, nous avons décidé de dégager une douzaine de configurations paramétriques statiques permettant de tester dynamiquement d’autres paramètres. Le premier paramètre statique est α, indicateur de l’importance relative accordée aux phéromones locales lors de la prise de décision d’une fourmi. A l’instar de quelques études paramétriques réalisées dans le cadre de l’application de la métaheuristique de type colonie de fourmis (comme par exemple dans [28]), ce paramètre sera toujours fixé à 1. Par contre, β, qui donne l’importance relative de la valeur heuristique, variera entre 1, 2 et 5. D’un autre côté, nous avons le facteur de persistance phéromonal, ρ, dont la valeur sera de 0,6 ou 0,8. 62

CHAPITRE 5. TESTS ET RÉSULTATS

63

Enfin, dans le but de tester l’impact des potentiels de distances sur la prédiction, nous avons décidé de tantôt utiliser les potentiels ds0 et ds1, tantôt de ne pas le faire. Pratiquement, nous avons donc créé 12 fichiers exécutables différents proposant une première configuration paramétrique statique définie. En effet, les paramètres α, β et ρ seront invariable tout au long de l’exécution d’un des fichiers exécutables, de même que les potentiels ds0 et ds1 seront ou non pris en compte. En ce qui concerne les autres paramètres, nous les avons fait varier dynamiquement et de la même manière dans chaque exécutable. Ainsi : – le nombre de fourmis varie de 50 à 500 par intervalle de 50 ; – le nombre de fourmis d’élite prends les valeurs de 0.00%, 0.01%, 0.10% et 0.50% du nombre de fourmis de la colonie ; – le facteur d’élitisme prends les valeurs de 2, 5 et 10 ; – et le nombre d’itération est de 100, 200 ou 500. De plus, nous avons fixé le nombre autorisé d’itérations sans modifications à 10% du nombre d’itérations total et à 0.10 l’amélioration énergétique minimale. Au final, nous obtenons donc 360 configurations dynamiques différentes. Dès lors, comme il y a 12 configurations statiques d’exécution, nous effectuerons 4320 prédictions. De plus, afin d’avoir un minimum d’aspect statistique pour nos résultats, nous avons exécuté chaque configuration 5 fois. En résumé, nous avons effectués 21600 exécutions de notre algorithme. Cela dit, découvrons-en les résultats.

5.2 Résultats obtenus 5.2.1 Evaluation énergétique de notre algorithme Le but à atteindre par notre algorithme est, rappelons le, la minimisation de la fonction énergétique utlisée. Cela signifie donc que si la configuration prédite présente une énergie libre minimale, l’objectif primaire sera atteint. Nous avons donc établi des graphiques reprennant l’énergie libre moyenne calculée dans une configurations algorithmique bien définie. En ne considérant dans un premier temps que les configurations ne prenant pas en compte les potentiels de distance, les graphiques de la figure 5.1 nous montrent très clairement que l’énergie minimale est généralement trouvée avec un nombre important de


64

fourmis. Malheureusement, aucun de ces graphiques ne fait ressortir de configuration préférentielle. Notons toutefois que le nombre d’itérations influe de manière générale sur la qualité des conformations prédites. En effet, il apparaît visiblement que pour tout autre paramètre équivalent, la moyenne d’énergie est plus basse quand le nombre d’itération est à 500. Ces observations sont valables pour toutes les configurations algorithmiques ne prenant pas en compte les potentiel de distance à l’exception toutefois de la configuration (α = 1, β = 1, ρ = 0.6) (cfr. figure 5.2). En effet, dans ce cas-ci, les courbes sont plutôt anarchiques et ne présentent aucune pente particulière. Cela pourrait s’expliquer par un comportement exploratoire trop poussé. En effet, dans cette configuration-ci, la valeur heuristique n’est pas mise en avant (β = 1) et la persistence phéromonale est faible. En ce qui concerne les prédictions effectuées avec la prise en compte des potentiels de distance, un comportement similaire se dégage avec néanmoins quelques nuances (cfr. figur 5.3). Premièrement, il est flagrant que l’énergie libre globale moyenne des prédictions effectuées dans ce cas-ci sont toutes suppérieures à celle calculée sans la prise en compte des potentiels de distance. Cela s’explique simplement par le fait que les potentiels de distance prennent en compte les distances spatiales et sont donc difficiles à minimiser. A nouveau, certaines configurations (celles avec β = 1) présentent des courbes plutôt anarchique (cfr. figure 5.4). A nouveau, cela s’explique par un comportement probablement trop exploratoire.

5.2.2 Calcul du RMSD En ce qui concerne l’évaluation par RMSD de nos prédictions, les résultats ont fort probablement souffert suite au calcul de la moyenne des résultats d’une configuration équivalente. En effet, lors des cinq exécutions à configuration identique, l’écart entre le RMSD le plus petit et le RMSD le plus élevé était généralement élevé. Ainsi, aucune moyenne intéressante ne ressort de notre analyse. Par exemple, nous constatons que la figure 5.5 offre des RMSD moyen compris entre 8Å et 16Å, ce qui est assez élevé. Par ailleurs, on aurait pu s’attendre à ce que la prise en compte des potentiels de distance améliore sensiblement le résultat du calcul RMSD. En effet, il serait logique que de par la prise en compte des interactions spatiales, une structure prédite présente plus de similitudes avec la structure réelle de la protéine étudiée. Malheureusement ce n’est pas du tout le cas, les courbes RMSD de toutes


65

F. 5.1 – Graphiques de l’énergie moyenne des prédictions effectuées avec une configuration de (α = 1, β = 1, ρ = 0.8, ds0 et ds1 ne sont pas pris en compte). Les graphique allignés se passent à un nombre d’itérations équivalent et les graphiques en colonne ont le même facteur d’élitisme. Les quatres courbes représentent chacune un nombre de fourmis d’élite différent.


66

F. 5.2 – Graphiques de l’énergie moyenne des prédictions effectuées avec une configuration de (α = 1, β = 1, ρ = 0.6, ds0 et ds1 ne sont pas pris en compte). Dans ce cas-ci, les diverses courbe montrent un comportement tout à fait disparate.


67

F. 5.3 – Graphiques de l’énergie moyenne des prédictions effectuées avec une configuration de (α = 1, β = 2, ρ = 0.8, ds0 et ds1 sont pris en compte). Les graphique allignés se passent à un nombre d’itérations équivalent et les graphiques en colonne ont le même facteur d’élitisme. Les quatres courbes représentent chacune un nombre de fourmis d’élite différent.


68

F. 5.4 – Graphiques de l’énergie moyenne des prédictions effectuées avec une configuration de (α = 1, β = 1, ρ = 0.6, ds0 et ds1 sont pris en compte). Dans ce cas-ci, les diverses courbe montrent un comportement tout à fait disparate.


69

les configurations utilisées présentent la même allure que celle montrée en figure 5.5. Par contre, l’observation des résultats obtenus individuellement montre que 59% des RMSD inférieurs à 6Å sont obtenus lors de la prise en considération des potentiels de distance. Dès lors, au vu des résultats de cette section et de la section précédente, il est légitime de se demander dans quelle mesure il ne faudrait pas accorder un poids prépondérant aux potentiels de distance. En effet, il se pourrait que dans notre algorithme, les valeurs énergétiques calculées par les potentiels de distance n’aient qu’une toute petite influence sur l’énergie libre globale de la protéine. Rappelons nous d’ailleurs que l’impact des potentiels de tortions tt0 et ttt1 est volontairement déiminué dans notre algorithme, afin justement de ne pas favoriser les interactions locales au détriment des interactions spatiales ! Pour conclure cette section, il est tout de même important de relever le meilleur RMSD calculé sur une de nos prédiction. En effet, nous avons obtenus une conformation protéique présentant un RMSD de 3,52Å avec la structure native correspondante, ce qui est assez encourageant ! A titre comparatif, [18] obtient, pour le même fragment de protéine, un RMSD minimal de 3,1Å. Bien entendus, les outlis et méthodes utilisés sont différents mais le résultat est là.

5.2.3 RMSD Vs énergie libre Une question intéressante serait de se demander quelle est la correspondance entre les meilleurs RMSD et les meilleurs énergies libre. Cette question permet en fait d’évaluer la qualité des fonctions énergétiques utilisées. En effet, dans le cas ou ces fonctions sont parfaites, obtenir une structure d’énergie minimale garantirait un RMSD faible. Les schémas ci-dessous (cfr. figure 5.6) nous montrent ce que donne une telle comparaison. Le premier reprend l’intégralité des résultats obtenus pour les 10800 exécutions ne prenant pas en compte les potentiels de distance. Néanmoins, le schéma correspondant dans le cas où ces potentiels serait actifs et tout à fait semblable. Afin d’y voir plus clair, nous avons isolé les RMSD inférieurs à 7Å ce qui donne le schéma de la figure 5.6(b). Nous constatons donc que pour une même énergie libre, les RMSD peuvent être différent ! Cela montre clairement l’imperfection des fonctions énergétiques utilisées ainsi que l’impact de la limitation de notre algorithme à 8 potentiels de deux types différents. En effet, d’autres types de potentiels existent et pourraient peut-être améliorer ces résultats.


70

F. 5.5 – Graphiques du RMSD moyen des prédictions effectuées avec une configuration de (α = 1, β = 5, ρ = 0.8, ds0 et ds1 sont pris en compte). L’allure de ces courbes n’évoque rien de particulier.


71

Les observations effectuées ci-dessus nous permettent par contre de mieux comprendre la diversité des résultats RMSD et donc les valeurs moyennes élevées obtenue pour ces résultats. A nouveau, la question d’une pondération particulière accordée à certains potentiel peut être soulevée.

F. 5.6 – Graphiques de comparaison entre le RMSD et l’énergie libre d’une conformation prédite : (a) prise en compte de toute les prédictions effectuées sans considération des potentiels de distance, (b) limitation aux prédictions présentant un RMSD inférieur à 7Å. A titre quantitatif, voici les cinq meilleurs scores énergétiques et leurs RMSD associé dans le cas où les potentiels de torsions ne sont pas pris en compte : Energie libre -34,8138 -34,5582 -34,3296 -33,8029 -33,6585

RMSD 11,5475 15,9341 15,2262 10,0224 15,58

et le meilleur RMSD (3,9989Å) possède un score énergétique de -20,5775 ! En considérant les potentiels de distances nous obtenons le tableau suivant :


72

Energie libre RMSD -26,0189 9,9723 -25,5462 9,6658 -24,7493 14,4044 -24,0495 14,8154 -23,8619 13,8676 et le meilleur RMSD (3,5222Å) possède un score énergétique de -10,6868 !

5.3 Notre meilleure prédiction : une RMSD de 3,5222Å Pour terminer en beauté, découvrons les étapes intermédiaires de notre meilleure prédiction. Le tableaux ci-dessous reprends l’évolution de notre prédiction en ce qui concerne les RMSD et le score énergétique. Notons que nous avons pu procéder à la réexécution de cette simulation grâce à la sauvegarde du numéro aléatoire généré par une fonction interne au angage C++. La prédiction a donc fourni les résultats suivants : Cycle Energie libre RMSD 1 29,182 9,51864 11 7,71167 10,5208 21 3,49845 9,18431 31 3,49845 9,18431 41 -9,41256 10,4728 51 -9,41256 10,4728 61 -9,41256 10,4728 71 -9,41256 10,4728 81 -9,41256 10,4728 91 -10,6507 5,55853 101 -10,6507 5,55853 111 -10,6507 5,55853 121 -10,6507 5,55853 131 -10,6868 3,52219 140 -10,6868 3,52219 Pour remettre cette prédiction dans son contexte, signalons qu’elle a été obtenue par une colonie de 100 fourmis dont la moité étaient des fourmis d’élite. Le facteur d’élitisme était de 5, ce qui implique une sérieuse prise en compte des


73

meilleurs résultat de chaque cycle ! En ce qui concerne les autres paramètre, α et β étaient à 1, ρ à 0,6. Enfin, le nombre d’itération maximal était fixé à 500. En image, la figure 5.7 nous montre la prédiction effectuée par AC-ProPre et la figure 5.8 la structure native réelle de l’ubiquitine (résidu 1-35). Enfin, la figure 5.9 nous montre ce que donne la meilleure superposition de ces deux conformations tridimentionnelles.

F. 5.7 – Représentation de la structure tertiaire de la prédiction effectuée par AC-ProPre.


74

F. 5.8 – Représentation de la structure tertiaire de l’ubiquitine (résidu 1-35).


75

F. 5.9 – Supperposition de la conformation simulée par AC-ProPre avec la conformation native connue.

Conclusion De nombreuses méthodes de prédiction de la structure native d’une protéine existent à ce jour. Parmi elles, la métaheuristique d’optimisation par colonie de fourmis propose un moyen particulier pour résoudre ce problème. Basé sur la communication indirecte entre les agents de recheche (les fourmis artificielles), cette méthode a déjà connu diverses implémentations, les plus fréquentes ayant été réalisées sur un modèle réseau utilisant comme représentation de la protéine le modèle HP [16, 36]. Cet ouvrage a présenté une nouvelle implémentation de cette métaheuristique faisant usage d’un modèle hors réseaux et d’une représentation détaillée de la protéines et des acides aminés la composant. Le but recherché était d’appliquer une méthode de recherche performante sur une modélisation fine du problème. Malgrés ses nombreux avantages, la métaheuristique d’optimisation par colonie de fourmi nécessite un temps considérable en tests afin d’en déterminer la configuration paramétrique. Nous avons donc effectué quelques choix afin de procéder à des test qui nous semblaient intéressants. Les résultats de ces tests sont extrèmement encourageant pour un premier essai. En effet, des prédictions tout à fait valables ont eu lieu. Malheureusement, nous n’imaginions pas l’étendue des tests faisables ainsi que l’influence des méthodes statistiques sur les résultats obtenus. Ainsi, de nomreux essais et expériences sont encore à faire. Une perspective intéressante serait de permettre l’autoadaptativité de nos paramètres. Cette méthode propose de faire de la configuration paramétrique une série de variable directement modifiable par l’algorithme en cours, en fonction des résultats obtenus. Une autre idée serait de procéder à la paramétrisation de l’algorithme sur base des résultats fournis par le calcul RMSD. Bien entendu, cela ne sera faisable que dans le cas d’une protéine dont la structure native est connue. Une fois que des paramètres satisfaisants auront été dégagé, les performances atteintes 76

Conclusion

77

pour une protéine du même ordre de grandeur que celle ayant servi à la paramétrisation seraient probablement élevées. Il sera ainsi intéressant d’étudier le comportement de AC-ProPre dans le cas de protéines plus importantes ou présentant des caractéristiques particulières. Il serait également possible de modifier le fonctionement de notre algorithme en faisant des choix stratégiques différents. En effet, il est certain que le comportement des fourmis sera modifié en fonction de la stratégie de dépôt phéromonal utilisée. De plus, le calcul de la valeur heuristique pourrait se baser sur d’autres données permettant ainsi une exploration différente de l’espace conformationnel. Un dernier point intéressant serait d’utiliser ce genre de méthode pour déterminer la robustesse des fonctions énergétiques. En effet, comme nous le découvrions dans le chapitre précédent, une fonction énergétique idéale devrait permettre d’obtenir, en fonction de sa minimisation, un RMSD également minimal. Notre travail constite donc une première étape fortement intéressante dans l’implémentation de métaheuristiques de recherche utilisant une représentation détaille de la protéine ainsi que de l’espace conformationnel dans lequel elle évolue.

Bibliographie [1] T.E. Creighton. Proteins : structure and molecular properties. W.H. Freeman and Company, New York, 1993. [2] J-PA. Kocher, MJ. Rooman, and SJ. Wodack. Factors influencing the ability of knowledge-based potentials to identify native sequence-structure matches. J. Mol. Biol., 234 :1598–1613, 1994. [3] C. Levinthal. Mossbauer spectroscopy in biological systems. In How to fold graciously, pages 22–24. Debrunner, P. and Tsibris J. and Munk, E. editors, 1969. [4] Nora Benhabilès, Annick Thomas, and Robert Brasseur. Les méchanismes de repliements des protéines solubles. Biotechnol. Agron. Soc. Environ., 4(2) :71–81, March 2000. [5] J.N. Onuchic, H. Nymeyer, A.E. Garcia, J. Chahine, and N.D. Socci. The energy landscape theory of protein folding : insights into folding mechanisms and scenarios. Advances in protein Chemistry, 53 :87–152, 2000. [6] HS. Chan and KA. Dill. Protein folding in the landscape perspective : chevron plots and non-arrhenius kinetics. Proteins : Struct. Funct. Genet., 30 :2–33, 1998. [7] HS. Chan and KA. Dill. From levinthal to pathways to tunnels. Nat. Stru. Biol., 4 :1O–19, 1997. [8] H.M. Berman, J. Westbrook, Z. Freng, G Gililand, H. Bhat, T.N. ans Weissig, I.N. Shindyalov, and P.E. Bourne. The protein data bank. Nucleic Acids Research, 28 :235–242, 2000. [9] G. Wagner, S.G. Hyberts, and T.F. Havel. Nmr structure determination in solution :a critique and comparison with x-ray crystallography. Annual Review of Biophysics and Biomolecular Structure, 21 :167–198, 1992.

78

BIBLIOGRAPHIE

79

[10] G. Wider. Structure determination of biological macromolecules in solution using nuclear magnetic resonance spectroscopy. Biotechniques, 29 :1278– 1290, 2000. [11] J. Kubelka, J. Hofrichter, and WA. Easton. The protein folding speed limit. Current Opinion in Structural Biology, 14 :76–88, 2004. [12] AC. May and TL. Blundell. Automated comparative modelling of protein structures. Current Opinion in Biotechnology, 5 :355–360, 1994. [13] A. Godzik. Fold recognition methods. Methods in Biochemical Analysis, 44 :525–546, 2003. [14] MJ. Rooman, J-PA. Kocher, and SJ. Wodak. Prediction of protein backbone conformation based on 7 structure assignments : influence of local interaction. J. Mol. Biol., 221 :961–979, 1991. [15] L. Holm and C. Sander. Dali/fssp classification of three-dimensional protein folds. Nucleic Acids Res., 25(1) :231–234, 1997. [16] KA. Dill, S. Bornberg, K. Yue, Y. Fiebig, D. Yee, P. Thomas, and H. Chan. Principles of protein folding - a perspective from simple exact models. Protein Science, 4 :561–602, 1995. [17] M. Levitt. A simplified representation of protein conformation for rapid simulation of protein folding. J MOL Biol, 104 :59–107, 1976. [18] D. Gilis and M. Rooman. Identification and ab initio simlations of early folding units in protins. Proteins : Struct. Funct. Genet., 42 :164–176, 2001. [19] N. Metropolis, AW. Rosenbluth, MN. Rosenbluth, and AH. Teller. Equation of state calculation by fast computing machines. J. Chem. Phys., 21 :21087– 1092, 1953. [20] S. Kirkpatrick, CD. Jr Gelatt, and MP. Vecchi. Optimization by simulated annealing. Science, 220 :671–680, 1983. [21] V. Granville, M. Krivanek, and J.P. Rasson. Simlated annealing : A proof of convergence. IEEE Transactions on Pattern Analysis and Machine Intelligence, 16(6) :652–656, 1994. [22] A. Colorni, M. Dorigo, F. Maffioli, V. Maniezzo, G. Righini, and M. Trubian. Heuristics from nature for hard combinatorial optimization problems. International Transaction in Operational Research, 3(2) :1–21, 1996. [23] Glover F. Tabu search : a tutorial. Special Issue on the practice of Mathematical Programming - Interfaces, 20(1) :74–94, 1990.

BIBLIOGRAPHIE

80

[24] DE. Goldberg. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison Wesley, 1989. [25] S. Ben Hamida. Algorithmes évolutionnaires : prise en compte des contraintes et applications réelles. Paris-Sud XI, 2001. [26] JL. Deneubourg, JM. Pasteels, and Verhaeghe JC. Probabilistic behaviour in ants : a strategy of errors ? Journal of Theoretical Biology, 105 :259–271, 1983. [27] A. Colorni, M. Dorigo, and V. Maniezzo. Distributed optimization by ant colonies. In Proceedings of European Conference on Artificial Life, pages 134– 142, 1991. [28] A. Colorni, M. Dorigo, and V. Maniezzo. An investigation of some properties of an ant algorithm. In Proceedings of the Parallel Problem Solving from Nature Conference, pages 509–520, 1992. [29] M. Dorigo, G. Di Caro, and LM. Gambardella. Ant algorithms for discrete optimization. Artificial Life, 5(3) :137–172, 1999. [30] M. Dorigo. Parallel ant system : An experimental study. Unpublished manuscript, 1993. [31] MG. Norman and P. Moscato. A competitive-cooperative approach to complex combinatorial search. In Proceedings of the 20th Joint Conference on Information and O.R., pages 315–329, 1991. [32] D. Beyer and R. Ogier. Tabu learning : a neural network search method for solving nonconvex optimization problem. In Proceedings of the International Joint Conference on Neural Network, 1991. [33] G. Ramachandran and V. Sasilekharan. Conformation of peptides and proteins. Adv. Protein Chem., 23 :283–438, 1968. [34] MJ. Rooman, J-PA. Kocher, and SJ. Wodak. Extracting information on folding from the amino acid sequence : accurate prediction for protein regions with preferred conformation in the abscence of tertiary interactions. Biochemistry, 31 :10226–10238, 1992. [35] Y. Dehouck, D. Gilis, and M. Rooman. A new generation of statistical potentials for proteins. Biophys J., 90(11) :4010–4017, June 2006. [36] A. Shmygelska and HH. Hoos. An ant colony optimisation algorithme for the 2d and 3d hydrophobic polar protein folding problem. BMC Bioinformatics, 6 :30–51, 2005. [37] K. Braden. A simple approach to protein structure prediction using genetic algorithms. Stanford University, CS 426 :36–44, 2002.

Algorithme d'optimisation par colonie de fourmis

Recommend Documents