Champs de force, énergie potentielle et autres joyeusetés en modélisation moléculaire

Introduction

Cet article, d'une très longue série, a pour but de don­ner plus de détails et peut-être — qui sait — vous titiller suf­fi­sam­ment pour enjam­ber la bar­rière phy­sique (au sens lit­té­ral) et rejoindre le fabu­leux monde d'Oz de la bio­in­for­ma­tique struc­tu­rale.

Avant de com­men­cer, je vous conseille de prendre du papier, un crayon, un thé ou un café. Ces breu­vages ser­vi­ront à illus­trer mes exemples et faire une pause dans votre lec­ture.

La bio­in­for­ma­tique struc­tu­rale, aus­si connue sous le nom bar­bare de bio­lo­gie struc­tu­rale com­pu­ta­tion­nelle, a pour objec­tif de com­prendre les fac­teurs d'influence qui déter­minent la fonc­tion bio­lo­gique. Quoi ? Vous faites ça tous les jours sur votre bureau et avec votre petit ordi­na­teur ? Je vais vous faire déjan­ter tout de suite. Com­prendre ces fac­teurs implique dès le départ d'avoir une bonne connais­sance des élé­ments sus­cep­tibles de jouer un rôle. Ain­si et comme vous le savez très bien : fonc­tion bio­lo­gique = struc­ture. Tou­te­fois, on ne s'intéresse pas exclu­si­ve­ment à la struc­ture finale. On va regar­der com­ment cette struc­ture se forme et pour cela, on peut regar­der les acides ami­nés conser­vés au cours de l'évolution et res­pon­sables de cer­taines caté­go­ries de replie­ments struc­tu­raux ou encore com­pa­rer les molé­cules homo­logues et non-homo­logues ce qui apporte éga­le­ment beau­coup d'informations car elles ren­seignent sur la conver­gence struc­tu­rale comme c'est le cas des pro­téines appar­te­nant aux familles de replie­ments dont cer­tains membres ont des scores de simi­la­ri­tés médiocres…

Depuis long­temps, on sait que l'information du replie­ment est conte­nue dans la séquence pro­téique ou nucléique que l'on peut obte­nir faci­le­ment aujourd'hui. La base de ces replie­ments réside dans les pro­prié­tés phy­si­co-chi­miques liées aux inter­ac­tions entre atomes. Le nombre de replie­ments semble donc infi­ni aux vues des nom­breuses confor­ma­tions pos­sibles. Heu­reu­se­ment pour nous, cer­tains replie­ments agissent comme une ancre pour ne for­mer qu'un nombre limi­té de confor­ma­tions stables et leurs détec­tions faci­litent le tra­vail de déter­mi­na­tion de la fonc­tion bio­lo­gique.

Ain­si, la bio­in­for­ma­tique struc­tu­rale uti­lise les modèles 1D (séquence) et 3D (struc­ture). L'histoire pour­rait s'arrêter là mais ce serait réduc­teur de pen­ser que le tra­vail se limite à faire un blast et un RMSD. En effet, la par­ti­cu­la­ri­té de s'intéresser au Vivant est justement…qu'il est Vivant, donc modu­lable. Cette modu­la­ri­té au niveau struc­tu­ral a pour consé­quence d'avoir des modules qui adoptent un replie­ment carac­té­ris­tique indé­pen­dam­ment du reste de la séquence et de manière indé­pen­dante. De plus, cer­tains replie­ments s'associent entre eux pour for­mer un domaine struc­tu­ral fonc­tion­nel mais par­fois un module seul suf­fit pour per­mettre au replie­ment d'être fonc­tion­nel. L'exemple le plus par­lant est un ARN de trans­fert, dont les modules s'associent entre eux pour for­mer un domaine spé­ci­fique d'un acide ami­né : si bien que même si l'on recon­naît un module par une matrice poids-posi­tions, cela n'indique en aucun cas son influence sur la struc­ture 3D. La pré­dic­tion et la simu­la­tion sont donc des solu­tions alter­na­tives à l'expérimentation tout azi­mut (qui est chère, longue et sou­vent hasar­deuse dans sa réa­li­sa­tion).

La modé­li­sa­tion molé­cu­laire englobe un grand nombre de tech­niques pour pré­dire la confor­ma­tion d'un sys­tème molé­cu­laire. Elles ont en com­mun de par­tir de fichiers de struc­tures connues et dis­po­nibles dans les bases de don­nées inter­na­tio­nales comme RCSB PDB ou PDBe pro­ve­nant d'expériences de Cris­tal­lo­gra­phie RX, RMN ou CryoEM et de recher­cher une nou­velle confor­ma­tion. Très uti­li­sé en Bio­lo­gie Struc­tu­rale, et notam­ment dans l'industrie phar­ma­ceu­tique, cette approche per­met éga­le­ment de mieux com­prendre les nom­breuses inter­ac­tions les plus pro­bables au sein d'une molé­cule (avec elle-même, c'est le cas lors de replie­ment secon­daire de pro­téines) ou avec son envi­ron­ne­ment ou bien encore avec un ligand. Tout cela per­met l'étude dyna­mique des molé­cules.

Cet article s'intéresse par­ti­cu­liè­re­ment à la méca­nique molé­cu­laire, il s'agit d'une des trois méthodes que l'on peut uti­li­ser pour faire de la modé­li­sa­tion de struc­ture pro­téique. Les deux autres impliquent soit d’utiliser la méca­nique quan­tique soit d’utiliser un méthode mixte entre la méca­nique molé­cu­laire et la méca­nique quan­tique. Dans les faits, la méca­nique quan­tique s'adresse à des petits sys­tèmes car les cal­culs sont com­plexes. La méca­nique molé­cu­laire s'adresse aux grands sys­tèmes car les cal­culs sont moins com­plexes. La méthode mixte com­bine à la fois la méca­nique molé­cu­laire et quan­tique : la par­tie de la molé­cule à étu­dier (un site actif par exemple) est simu­lée via méca­nique quan­tique tan­dis que le reste de la molé­cule est simu­lé par méca­nique molé­cu­laire.

La pre­mière étape est de bien com­prendre ce que l’on modé­lise. Repre­nons les bases et répé­tez avec moi : on modé­lise des molé­cules ! Bien…mais ce n'est pas fini ! Une molé­cule est faite d'atomes et chaque atome inter­agit avec son voi­sin. Au lieu de s'échanger de la farine ou du sel, des atomes proches mettent en com­mun des élec­trons et bim : ça fait une inter­ac­tion. Et voi­là, le pre­mier nœud à pro­blèmes : modé­li­ser une inter­ac­tion entre atomes. La solu­tion : uti­li­ser un champ de force.

Champs de force

Pour bien com­prendre ce qu'est un champ de force, je vais prendre un exemple très par­lant. Faites un des­sin hyper-sim­pli­fié de vous même debout : un point pour la tête, deux points pour les mains, un point pour une ver­tèbre tho­ra­cique et un autre pour une ver­tèbre lom­baire, deux points pour vos fesses et enfin, deux points pour vos pieds. Reliez les points ensemble : chaque pied est relié à une fesse, vos fesses sont reliées à la ver­tèbre lom­baire. Celle-ci est reliée à la ver­tèbre tho­ra­cique. Et la tête et cha­cune des mains sont reliées à cette ver­tèbre. Plu­tôt que de décrire toutes les carac­té­ris­tiques internes (vos petits os, vos organes, vos cel­lules), on va par­tir du pos­tu­lat que ces points suf­fisent ample­ment. Si vous avez une bonne notion du corps humain, vous avez sans doute remar­qué que votre ver­tèbre tho­ra­cique est à mi-dis­tance de cha­cune de vos mains et que cette ver­tèbre est plus proche de la tête que de vos fesses. Sinon, je vous invite à consul­ter.

Évi­de­ment, cer­tains atomes pré­sentent des carac­té­ris­tiques propres comme leur charge (due à leurs élec­trons de valence) ou leur taille (due aux nombre de nucléons). Si on va plus loin, il y a éga­le­ment des inter­ac­tions entre atomes très carac­té­ris­tiques et qui jouent un rôle dans les forces qui s'exercent sur cha­cun des atomes liés. De même que des forces s'exercent entre atomes dis­tants, on parle alors d'interactions entre atomes non liés.

Si on applique cela à la modé­li­sa­tion du corps humain, on peut dire que le point repré­sen­tant la tête dis­pose de carac­té­ris­tiques dif­fé­rentes de celui qui repré­sente une main ou un pied. De même, on peut dire que la tête est liée à la ver­tèbre tho­ra­cique mais il existe des forces qui la lie de manière dis­tante aux pieds. Ouf pour nous !

Comme vous vous en dou­tez, les échanges qui ont lieu au sein d'une inter­ac­tion sont d'ordres éner­gé­tiques. Ain­si, les forces qui s'exercent sur les atomes ont une action papable sur des valeurs phy­siques réelles comme des dis­tances ou des angles. Autre­ment dit, une force va appor­ter une éner­gie par­ti­cu­lière pour modi­fier l'interaction elle-même : c'est l'énergie poten­tielle. Celle-ci, a, par défi­ni­tion, le poten­tiel d'être trans­for­mée en éner­gie ciné­tique donc en mou­ve­ment. Et la manière de repré­sen­ter la varia­tion est d'avoir recours à une fonc­tion ana­ly­tique. Pour les adeptes des mathé­ma­tiques, il s'agit de la fonc­tion d'une variable qui est déve­lop­pable en série au voi­si­nage d'un inter­valle don­né. C'est-à-dire que sur un inter­valle don­né, la fonc­tion géné­rale est égale à la somme de fonc­tions simples.

Repre­nons mon exemple. Votre sché­ma de base peut s'appliquer à tous les types d'humanoïdes . Tou­te­fois, de petites dif­fé­rences sont notables : votre voi­sin sera plus grand ou plus petit que vous par exemple. On a donc des variables qu'il faut prendre en compte pour décrire les inter­ac­tions entre cha­cun des points de votre sché­ma.

Pour résu­mer : un champ de force est la somme des forces (répul­sions et attrac­tions) qui s'exercent sur une inter­ac­tion, autre­ment dit sur les atomes qui com­posent l'interaction. Cha­cune de ces forces se décrit par son éner­gie poten­tielle qui dépend de constantes décrites de manières empi­riques et de variables qui dépendent du modèle struc­tu­ral (du sché­ma que vous avez fait).

Il existe de nom­breux champs de force et cha­cun dif­fère par la valeur des constantes décri­vant le type d'interaction qui lui est asso­cié. Elles sont four­nies par la spec­tro­sco­pie infra­rouge, la RMN ou par des cal­culs quan­tiques.

Pour repar­tir de mon exemple, cha­cun des traits que vous avez fait entre les points de votre sché­ma est un champ de force. Réjouis­sant !

Coordonnées cartésiennes et coordonnées internes

Main­te­nant vient la dou­lou­reuse ques­tion de la repré­sen­ta­tion. Doit-on repré­sen­ter l'interaction ou les atomes qui la com­posent ? Vous allez me dire : c'est la même chose. En fait, non ! On peut soit décrire la molé­cule par la posi­tion selon les 3 coor­don­nées x, y et z en ang­ström (10^-12m) de cha­cun des atomes, on uti­lise alors des coor­don­nées car­té­siennes, soit décrire la géo­mé­trie de la molé­cule en expri­mant les lon­gueurs de liai­sons, les angles de valances et de tor­sion, on parle alors de coor­don­nées internes ou Z‑matrice ou encore appe­lé modes de vibra­tions internes et que pour évi­ter la redon­dance des expres­sions matri­cielles, le pre­mier atome est l'origine du réfé­ren­tiel.

On a vu qu'une éner­gie poten­tielle peut engen­drer une modi­fi­ca­tion d'un des para­mètres de l'interaction. Ces modi­fi­ca­tions sont les coor­don­nées internes et sont répar­ties en deux grandes caté­go­ries : les vibra­tions d'allongements (appe­lé aus­si vibra­tions de valence ou modes d'extension ou stret­ching) et les vibra­tions de défor­ma­tions (appe­lé aus­si mode de défor­ma­tion ou ben­ding). Par­mi les stret­ching, il y a les vibra­tions symé­triques (Vs) et asy­mé­triques (Vas). Par­mi les ben­ding, il y a les défor­ma­tions dans le plan de molé­cule : la rota­tion plane (notée β et appe­lée aus­si rocking) et le cisaille­ment (noté δ et appe­lé aus­si scis­so­ring) ain­si que les défor­ma­tions hors du plan de la molé­cule : le balan­ce­ment (noté ω et appe­lé aus­si wag­ging) et la tor­sion (notée τ et appe­lée aus­si twis­ting).

Dans un sys­tème 3D, un objet peut être repré­sen­té selon sa posi­tion (x, y et z) et sa quan­ti­té de mou­ve­ment. C'est ce que l'on appelle les degrés de liber­té. Ain­si, on a 3 degrés pour la posi­tion et 3 degrés pour la rota­tion chez les molé­cules linéaires (2 chez les linéaires, car un axe est blo­qué) et un cer­tain nombre de coor­don­nées internes. On peut asso­cier une coor­don­née interne à un degré de liber­té et puisqu'un cer­tain nombre a déjà été comp­ta­bi­li­sé, on peut dénom­brer les coor­don­nées internes pour chaque type de molé­cule. On a donc 3n‑5 coor­don­nées internes (n étant le nombre d'atomes) pour une molé­cule linéaire et 3n‑6 coor­don­nées internes pour une molé­cule non-linéaire. Par exemple, le dioxyde de car­bone (CO₂) est une molé­cule linéaire, une de ces rota­tions est blo­quée par la nature des liai­sons qui engagent le Car­bone : elle n'a donc que 4 coor­don­nées internes. En revanche, le méthane (CH₄) est une molé­cule non-linéaire, elle a 9 coor­don­nées internes.

Aujourd'hui, le for­mat de fichier PDB — le plus uti­li­sé — contient exclu­si­ve­ment des don­nées sous forme de coor­don­nées car­té­siennes. Vous pour­rez retrou­ver la notice de ce for­mat ici.

Un petit exemple tiré de la struc­ture d'une pro­téase pan­créa­tite bovine. J'ai volon­tai­re­ment extrait uni­que­ment les coor­don­nées car­té­siennes d'une pro­line de cette molé­cule.

Pour les coor­don­nées internes, on peut uti­li­ser des logi­ciels (VDM) ou des packages R (Rpdb) ou encore Python (Bio.PDB). Fina­le­ment, j'ai uti­li­sé VDM. J'ai obte­nu toutes les dis­tances entre atomes (7), cha­cun des angles (8) et des angles dièdres (6). Cela cor­res­pond à la géo­mé­trie de notre petite pro­line. Du coup, comme je vous ai par­lé de 3n‑6 coor­don­nées internes utiles, on a 21 coor­don­nées internes cal­cu­lées par VMD soit au final, 6 valeurs redon­dantes. À l'échelle d'une petite molé­cule, ça fait peut mais quand on consi­dère des macro­mo­lé­cules de plu­sieurs mil­liers d'atomes, ça change la donne… et notam­ment au niveau for­mat de fichier et de la com­pres­sion de don­nées. Dans les faits, cela ne sert stric­te­ment à rien d'avoir la géo­mé­trie pour un acide ami­né, il est plus inté­res­sant d'obtenir les angles dièdres du sque­lette car­bo­né d'une pro­téine (back­bone car­bone) et de les repré­sen­ter via un dia­gramme de Rama­chan­dran. Je revien­drai des­sus pour un pro­chain article.

Calcul de l'énergie potentielle en mécanique moléculaire

En méca­nique molé­cu­laire, les noyaux et les élec­trons forment une seule enti­té : un point ponc­tuel. Ils sont tous sphé­riques, de charge nette (c'est-à-dire qu'elle est la même pour tous les atomes) et leurs tailles varient selon l'atome.

On a vu plus haut que l'énergie poten­tielle est en fait l'énergie néces­saire à appor­ter pour qu'une inter­ac­tion soit modi­fiée dans un espace tri-dimen­sion­nel. Un atome dans une molé­cule pré­sente de nom­breuses inter­ac­tions : des inter­ac­tions entre atomes liés et non-liés. Cha­cun de ces termes peut être repré­sen­té par une éner­gie poten­tielle. Au sein d’interactions entre atomes liés, on a une éner­gie poten­tielle de défor­ma­tion de liai­son (stretch), de défor­ma­tion d'angles (bend) et de tor­sion d'angles dièdres (dihe­dral angle). Dans le cas d’interactions entre atomes non-liés, on a une éner­gie poten­tielle élec­tro­sta­tique (poten­tiel de Cou­lomb), des liai­sons hydro­gènes et des inter­ac­tions de Van Der Waals (poten­tiel de Len­nard-Jones). Notez que dans cer­tains cas, une éner­gie poten­tielle pour les angles de tor­sion impropres peut être uti­li­sées en plus !
À des fins de rapi­di­té de cal­cul, il existe pour chaque champ de force des termes de cou­plage qui asso­cient deux poten­tiels pour cer­taines inter­ac­tions typiques comme par exemple : une élon­ga­tion des liai­sons quand l'angle dimi­nue (stretch-bend) ou lorsque une confor­ma­tion est éclip­sée (stretch-tor­sion).

Au final, cela donne ces équa­tions, avec en vert, les para­mètres variables et en rouge, les para­mètres empi­riques.

Avec les indi­ca­teurs :

Ka : constance de force de liai­son

l et l₀ : dis­tance inter­ato­mique mesu­rée et à l'équilibre

Kb : constante de force d'angle

θ et θ₀ : angle mesu­ré et à l'équilibre

V : constante de force de tor­sion asso­cié à n

n : ordre de séries de Fou­rier

Φ : angle mesu­ré dièdre

γ : angle de phase du dièdre

q_i et q_j : charges par­tielles de l'atome i et j

ε : constante diélec­trique

r_ij : dis­tance inter­ato­mique entre les atomes i et j

A_ij : pro­fon­deur de puits du poten­tiel appro­ché

r*_ij : dis­tance inter­ato­mique au mini­mum d'énergie

Et enfin la fameuse équa­tion repré­sen­tant l'énergie poten­tielle totale :

Le sol­vant est pris en compte par le choix de la valeur de la constante diélec­trique.

Petit détail qui a son impor­tance : il est exclu de vou­loir com­pa­rer l'énergie poten­tielle de deux molé­cules même si elles pré­sentent une forte simi­la­ri­té de confor­ma­tion car de nom­breux para­mètres d'un champ de force à l'autre ne sont pas super­po­sables. Enfin, sachez que selon le logi­ciel de cal­cul de champs de force uti­li­sé, les uni­tés uti­li­sées ne sont pas les mêmes : en kcal/​mol/​A pour Charmm ou en kJ/​mol/​nm pour GROMOS.

Si vous êtes encore vivants et curieux d'en savoir un peu plus, je vous invite au pro­chain article sur l'analyse des confor­ma­tions et de la mini­mi­sa­tion. Par ailleurs, pour les afi­cio­na­dos du code et autres adeptes de la gâchette du cla­vier, sachez que Bio­Py­thon dis­pose d'un ensemble de packages que l'on uti­lise pour l'étude des struc­tures PDB : Bio.PDB. De toute manière, les prin­ci­paux logi­ciels de visua­li­sa­tion de struc­tures sont en Python…

Voi­là, c'est tout pour le moment. Vous venez de voir com­ment fonc­tionnent glo­ba­le­ment les pro­grammes et les logi­ciels de champs de force en modé­li­sa­tion molé­cu­laire. La pro­chaine fois, je vous mon­tre­rai com­ment on peut s'en ser­vir pour recher­cher une confor­ma­tion par­ti­cu­lière.