1 Introduction

L'utilité de l'analyse de la dynamique interne des macromolécules biologiques pour la compréhension des processus impliqués dans leurs mécanismes d'action n'est plus à démontrer. De ce fait, les mesures de relaxation hétéronucléaire — plus spécialement celles des noyaux ¹⁵N — se sont particulièrement développées au cours de cette dernière décennie (pour une revue récente, voir : [1]). L'analyse de ces données en termes de dynamique interne peut se faire par différentes approches, allant de l'interprétation directe de l'allure des fonctions de densités spectrales (Spectral density mapping) [2–6], jusqu'à une analyse plus élaborée par divers formalismes, permettant l'obtention des paramètres du mouvement. Dans cette dernière approche, le formalisme de Lipari–Szabo (Model Free) [7] est probablement le modèle le plus populaire et le plus utilisé. La relevance des résultats obtenus par cette approche repose néanmoins sur la qualité et la quantité des données expérimentales recueillies (R₁, R₂ et NOE hétéronucléaires, en général). Bien que relativement simple, l'ajustement du modèle aux données expérimentales sera d'autant plus précis et fiable que les données expérimentales fournies seront plus nombreuses, surtout dans le cas de mouvements complexes nécessitant l'ajustement de nombreux paramètres (modèle de Lipari–Szabo étendu) [8,9]. De ce fait, une analyse des données de relaxation mesurées à plusieurs intensités de champ magnétique paraît être une approche souhaitable [10,11]. Nous présentons ici un programme MATLAB (DYNAMOF, pour DYNamics Analysis with MOdel-Free) permettant l'analyse de telles données. Ce programme permet de calculer les densités spectrales réduites J(0), J(ω_N) et < J(ω_H) > à partir des données expérimentales précitées, puis de les ajuster avec le(s) formalisme(s) de Lipari–Szabo. Nous présenterons un exemple d'application sur des données de relaxation mesurées à trois inductions B₀ différentes sur la protéine oncogène P13^MTCP1 [12,13].

2 Matériels et méthodes

Théorie : quand la relaxation d'un hétéronoyau X dépend essentiellement de l'interaction dipolaire avec un proton directement lié et de l'anisotropie de son déplacement chimique, les données de relaxation peuvent être interprétées en termes de mouvement du vecteur X–¹H [14]. Sachant que les trois paramètres de relaxation mesurés expérimentalement (R_X(X_z) ou R₁, R_X(X_xy) ou R₂, et X{¹H}NOE) dépendent des valeurs des fonctions de densité spectrale à cinq fréquences différentes [2,3], le calcul des densités spectrales peut être approché par l'approximation des hautes fréquences permettant un échantillonnage « réduit » de la fonction de densité spectrale. L'approximation des hautes fréquences implique une quasi égalité des valeurs des densités spectrales à haute fréquence (J(ω_H – ω_N) ≈ J(ω_H + ω_N) ≈ J(ω_H) = < J(ω_H) >, ce qui impose un temps de corrélation relativement important (> 1 ns) et des valeurs d'inductions de champ magnétique relativement élevées pour réaliser les mesures de relaxation (en général ≥ 9,4 T). Ces conditions sont généralement réalisées lors de l'étude de la dynamique interne d'une protéine en solution. Si ces conditions sont réunies, les vitesses de relaxation sont directement traduites en densités spectrales à trois fréquences [4–6] par la relation :

Cette approche, dite de « cartographie », de la fonction de densité spectrale permet de juger rapidement de la dynamique interne de la protéine. En effet, une protéine en solution est soumise aux lois de la thermodynamique : elle représente un système à l'équilibre dont l'énergie est conservée. Cela se traduit par la conservation de la quantité de mouvement pour chacun des vecteurs (¹⁵N–H ou ¹³C–H) caractéristiques de chacun de ses résidus, et donc la conservation de l'aire sous la courbe des fonctions de densité spectrale propres à chacun de ces vecteurs (voir Fig. 2). Ainsi, les zones flexibles de la protéine sont caractérisées par des valeurs faibles de J(0) (fréquence représentative des mouvements lents), contrebalancées par des valeurs élevées de < J(ω_H) > (fréquence représentative des mouvements rapides). À l'inverse, les zones rigides présentent des valeurs élevées de J(0) et faibles de < J(ω_H) >. Certains résidus peuvent déroger à cette règle à cause d'une contribution adiabatique à leur relaxation, généralement apportée par des phénomènes d'échange chimique ou conformationnel s'établissant sur une échelle de temps supérieure à celle du temps de corrélation de la molécule (μs–ms). Bien que parfaitement rigoureuse, cette analyse des fonctions de densité spectrale n'apporte qu'une information qualitative sur la dynamique de la protéine étudiée : aucune information sur les échelles de temps caractéristiques des mouvements qui animent chaque résidu ne peut être obtenue. Néanmoins, elle peut être facilement complétée par l'utilisation d'un modèle de mouvement permettant d'ajuster aux mieux les densités spectrales des différents vecteurs. L'approche model-free de Lipari et Szabo [7] peut être ainsi utilisée pour ajuster les valeurs expérimentales des fonctions de densités spectrales et obtenir les paramètres du mouvement. Du fait de la différence importante de leur domaine de temps respectif, ce formalisme fait l'hypothèse d'une contribution indépendante du mouvement global (ou des mouvements globaux) et des mouvements internes à la fonction d'autocorrélation des vecteurs X–¹H. Pour une protéine qui se réoriente de façon isotrope, on obtient :

Pour certains résidus, l'équation [2] s'avère insuffisante pour ajuster les valeurs expérimentales des densités spectrales : cela notamment pour des résidus qui présentent des mouvements internes sur un domaine de temps proche de 1 ns. Dans ce cas, l'expression de la fonction de densité spectrale est étendue à [8,9] :

Enfin, certains résidus présentent des valeurs de J(0) anormalement élevées, signant la présence d'une contribution adiabatique à la fonction de densité spectrale. Cette contribution est due à la présence de mouvements sur une échelle de temps allant de la micro- à la milliseconde (plus lent que le temps de réorientation de la molécule), généralement reliés à un échange conformationnel et essentiellement perçus par la mesure de R₂. Pour de tels résidus, la valeur de J(0) doit être corrigée de cette contribution d'échange (R_éch) avant d'ajuster la fonction de densité spectrale par le modèle de Lipari–Szabo. Cette correction peut être apportée grâce à une relation linéaire existant entre la valeur de 2 R₂ – R₁ et les fréquences de Larmor hétéronucléaires ω_X² [17,23] :

La complexité croissante de ces différents modèles entraîne une augmentation très significative du nombre de paramètres permettant leur description. Ainsi, généralement deux paramètres sont à ajuster pour une analyse avec le modèle de Lipari–Szabo simple, contre trois — voire quatre ! — avec le modèle de Lipari–Szabo étendu. Or, les données de relaxation classiquement enregistrées à un seul champ magnétique ne permettent d'obtenir que seulement trois valeurs de densité spectrale : J(0), J(ω_N) et <J(ω_H)>.

Le problème devient très vite sous-estimé (du moins partiellement), et le recours à un meilleur échantillonage de la fonction densité spectrale nécessaire. La simulation reportée sur la Fig. 1 montre bien qu'en effet, si l'analyse à une fréquence est suffisante pour détecter la contribution d'un mouvement interne au mouvement global d'un vecteur, seul le recours à une analyse à plusieurs fréquences est capable de révéler la relative complexité de ce mouvement interne.

Le programme DYNAMOF présenté dans ce manuscrit utilise cette stratégie pour analyser les données de relaxation enregistrées à plusieurs intensités de champ magnétique. Une procédure de type SIMPLEX a été adoptée pour ajuster les valeurs de densités spectrales expérimentales aux différentes variantes du modèle de Lipari–Szabo, dont la convergence est assurée par la minimisation d'un χ² calculé entre les valeurs expérimentales et les valeurs théoriques obtenues avec les différents modèles. Le choix du modèle pertinent se fera ensuite en comparant les valeurs de χ² obtenues pour les différents ajustements grâce à un test de Fischer–Snedecor (F-test).

Mesures expérimentales : Les expériences de RMN ont été réalisées à 9,4, 11,75 et 14,1 Tesla sur des spectromètres Bruker Avance 400, 500 et 600 équipés de sondes 5 mm triple résonance ¹H–¹³C–¹⁵N munies de gradient-z, et sur un échantillon de protéine P13^MTCP1 uniformément enrichie en ¹⁵N dissoute à une concentration de 1 mM dans 500 μl de tampon Tris-HCl 10 mM pH 7,0 (5% ²H₂O pour le lock). La température de mesure a été soigneusement calibrée à 20 °C sur chacun des spectromètres. Les séquences d'impulsions utilisées pour déterminer les valeurs des vitesses de relaxation ¹⁵N R_N(N_z) (R₁), R_N(N_xy) (R₂), et de ¹⁵N{¹H}NOE sont similaires à celles déjà décrites [2,3], et les paramètres expérimentaux utilisés ainsi que le traitement des données ont déjà été reportés en détail pour d'autres protéines étudiées au laboratoire [13,18]. La vitesse de relaxation longitudinale ¹⁵N (R_N(N_z)) a été obtenue à partir de dix expériences d'inversion–récupération, avec des délais de relaxation allant de 18 à 1026 ms. La vitesse de relaxation transverse ¹⁵N (R_N(N_xy)) a été obtenue à partir de dix expériences CPMG, avec des délais de relaxation allant de 16 à 144 ms. La valeur du NOE hétéronucléaire ¹⁵N{¹H} a été déterminée à partir du rapport de deux expériences réalisées avec et sans présaturation des protons. Les spectres RMN ont été traités à l'aide du logiciel Gifa (version 4.4) [21]. Le résultat de ces mesures expérimentales est reporté sur la Fig. 2.

3 Résultats

La Fig. 3 présente l'organigramme général du programme DYNAMOF. Le programme débute en demandant un suffixe (par exemple out), qui sera attribué à tous les fichiers de sortie (fichier out) créés lors de son utilisation. Après avoir offert la possibilité de changer en cours d'utilisation certaines valeurs de constantes entrées en valeurs par défaut (Δσ, valeurs de départ utilisées dans les procédures SIMPLEX...), l'utilisateur doit entrer le type de noyau utilisé pour l'analyse (¹⁵N ou ¹³C), puis le nombre de champs magnétiques utilisés. Après avoir laissé le choix entre l'analyse classique (Spectral density mapping) et celle récemment proposée par Roumestand et al. (Fast spectral density mapping, soumis à publication) que nous ne développerons pas ici, l'utilisateur doit entrer les différentes fréquences de mesures protons en ordre croissant (400, 500, puis 600) ainsi que le nom générique (gnam) des fichiers d'entrée correspondant à ces fréquences de mesure. Ces fichiers d'entrée consistent en trois fichiers textes pour chaque fréquence de mesure (R₁, R₂ et NOE (I/I₀)), correspondant aux fichiers directement obtenus par l'analyse des données RMN par le module « Relaxation » du programme de traitement des données Gifa 4 [21]. Ces fichiers sont des tableaux à trois colonnes, chacune contenant respectivement le numéro du résidu (position dans la séquence), la valeur du paramètre de relaxation, et l'erreur mesurée sur ce paramètre. À chacun de ces fichiers, on attribue un suffixe (en général la fréquence de mesure) qui constituera son nom générique (par exemple : R₁_400, R₂_400 et NOE_400 correspondent aux mesures de relaxation effectuées sur P13^MTCP1 à 400 MHz (generic name = 400)). Pour chaque fréquence, le programme calcule les valeurs des densités spectrales réduites en utilisant l'équation (1) ainsi que les erreurs correspondantes, qui seront stockées dans deux fichiers sortie différents (par exemple : J 400.out et dJ 400.out). Une sortie graphique reportant, pour chaque champ de mesure, les valeurs de densités spectrales réduites en fonction de la séquence de la protéine permet d'apprécier rapidement la qualité des données (Fig. 4).

À partir de cette étape, le programme propose de corriger les valeurs de J(0) d'une éventuelle contribution d'échange. Si l'utilisateur accepte ce choix, le facteur d'échange ϕ sera calculé par régression linéaire en utilisant l'équation [5] : un affichage graphique des données pour chaque résidu est proposé (Fig. 5a), puis le programme va afficher la valeur de ϕ en fonction de la séquence (Fig. 5b) et créer deux fichiers (PHI_R.out et Rex_R.out), correspondant l'un aux valeurs de ϕ par résidu, l'autre aux valeurs de R_éch calculées à partir de ϕ [éq. (6)] pour chaque champ de mesure. Les erreurs sur ϕ seront calculées par une procédure Monte-Carlo (nombre d'itérations choisies par l'utilisateur). Le fichier PHI_R.out contient également les valeurs corrigées de J(0) (et dJ(0)) : ces valeurs sont introduites dans deux nouveaux fichiers, contenant les valeurs de densités spectrales et les erreurs correspondantes : J_gnamcorR.out et dJ_gnamcorR.out (ex J_400corR.out et dJ_400corR.out pour la fréquence 400 MHz). Ce sont dès lors ces fichiers qui seront utilisés par défaut pour l'ajustement avec les différents modèles de mouvement.

Cette étape de calcul (et éventuellement de correction) des fonctions densités spectrales effectuées, le programme va débuter l'ajustement des données avec les différents modèles de mouvement décrits dans la partie Matériel et méthodes. En premier lieu, le programme va demander à l'utilisateur de choisir entre un ajustement « local » (un τ_c par résidu) ou « global » (le même τ_c pour tous les résidus de la protéine) des données : l'ajustement local des données peut s'avérer judicieux lorsque la protéine présente des domaines avec des comportements dynamiques très différents, ou lorsque la réorientation de la protéine en solution ne peut plus être considérée comme isotrope [11,19,20].

• • Ajustement Local (une valeur de τ_c par résidu). Avant de démarrer la procédure d'optimisation, le programme offre encore à l'utilisateur la possibilité d'ajuster la contribution d'échange (ϕ : Eq. (8)) si la correction proposée à l'étape précédente n'a pas été réalisée : les paramètres τ_c, τ_f, S² et éventuellement ϕ seront alors obtenus résidu après résidu par ajustement des 2 n + 1 (n étant le nombre de champs magnétiques utilisés) valeurs expérimentales de densité spectrale (J(0), J(40), J(50), J(60), < J(400) >, < J(500) > et < J(600) > dans le cas de p13^MTCP1) par le modèle de Lipari–Szabo (Eq. (2), légèrement modifiée : la contribution haute fréquence au mouvement est supposée indépendante de la fréquence et, de ce fait, ω τ négligé [11]). Une sortie graphique est proposée pour chaque résidu, permettant d'apprécier la qualité de chacun des ajustements. Deux fichiers sont créés à la fin de la procédure, contenant le résultat de chacun des ajustements (LOCAL.out, optimisation de τ_c, τ_f et S²–LOCALREX.out, optimisation de τ_c, τ_f, S² et ϕ). Un χ² expérimental est calculée pour chaque modèle et pour chaque résidu, dont le résultat est également contenu dans le fichier correspondant : ces valeurs de χ² pourront être utilisées pour sélectionner le modèle optimum grâce à une analyse statistique. Enfin, un calcul de type Monte-Carlo va permettre de déterminer les erreurs sur chacune des variables (nombre d'itérations choisi par l'utilisateur) : ces erreurs seront contenues dans deux fichiers « définitifs », à côté des valeurs des variables τ_c, τ_f, S² et ϕ (LOCAL_MC.out et LOCALREX_MC.out). Cette étape crée également deux fichiers supplémentaires contenant pour chaque résidu les valeurs de χ² obtenues après chaque itération et ce pour chacun des deux modèles (χ2_CAN_MC.out et χ2_CANREX_MC.out). Ces fichiers pourront être utilisés pour une analyse statistique plus poussée (en cours d'implémentation dans le programme).
• • Ajustement global (une valeur de τ_c commune à tous les résidus). Dans une première étape, le programme va tenter de déterminer la valeur commune du τ_c en ajustant résidu par résidu les valeurs des densités spectrales expérimentales par le modèle de Lipari–Szabo [Eq. (2)] et en tentant d'optimiser pour chaque résidu les valeurs de τ_c, τ_f et S². Une sortie graphique est proposée pour chaque résidu, permettant d'apprécier la qualité de l'ajustement (Fig. 6a). En fin de procédure, un fichier LS3.out est créé, contenant les résultats de l'ajustement pour chaque résidu. La terminologie LS3 signifie que ce fichier a été obtenu à partir du modèle de Lipari–Szabo (LS), en optimisant trois paramètres : τ_c, τ_f et S². Une sortie graphique présente les valeurs de τ_c en fonction de la séquence et propose une valeur moyenne (± 2σ) utilisable comme valeur de τ_c « global » (Fig. 6b). Cette valeur de τ_c pourra être alors utilisée comme constante dans le reste de la procédure qui va consister à ajuster les 2 n + 1 valeurs expérimentales de densités spectrales avec le modèle de Lipari–Szabo simple [Eq. (2)] (LS2, deux variables à optimiser : τ_f et S²), le modèle de Lipari-Szabo simple avec contribution d'échange [Eq. (2) et (8)] (LSREX, trois variables à optimiser : τ_f, S² et ϕ), et les modèles de Lipari–Szabo étendus [Eq. (3) et (4)] (LSE3, trois variables à optimiser : S²_s, S²_f et τ_s – LSE4, quatre variables à optimiser : S²_s, S²_f, τ_s et τ_f). Comme précédemment, une sortie graphique est proposée résidu par résidu, permettant de juger la qualité de l'ajustement par chaque modèle (Fig. 7). Quatre fichiers sont alors créés (LS2.out, LSREX.out, LSE3.out et LSE4.out), contenant les résultats obtenus avec les modèles respectifs. Tout comme pour l'approche « locale », un χ² expérimental est calculée pour chaque modèle et pour chaque résidu, dont le résultat est également contenu dans le fichier correspondant : ces valeurs de χ² pourront être utilisées pour sélectionner le modèle optimum grâce à une analyse statistique. De même, les erreurs sur chacun des paramètres peuvent alors être calculées par une procédure Monte-Carlo (l'utilisateur choisit le nombre d'itérations) : cette étape conduira à la création de quatre nouveaux fichiers contenant à la fois les valeurs des paramètres optimisés et les erreurs calculées sur ces paramètres (LS2_MC.out, LSREX_MC.out, LSE3_MC.out et LSE4_MC.out). De même, quatre fichiers (correspondant aux quatre modèles) seront créés contenant les valeurs de χ² obtenues à chaque itération du calcul de Monte-Carlo pour chaque résidu (χ2_LS2_MC.out, χ2_LSREX_MC.out, χ2_LSE3_MC.out, χ2_LSE4.out) : ces derniers seront utilisés pour une analyse statistique plus poussée qui est en cours d'implémentation dans le programme.

On notera que chaque fois qu'un ajustement par un modèle est effectué, que ce soit dans la procédure « locale » ou « globale », un fichier contenant les valeurs des densités spectrales théoriques (et, après un calcul Monte-Carlo, un fichier contenant les erreurs sur les densités spectrales théoriques) est créé, du type Jtheo_nomdumodèle.out (ou jtheo_nomdumodele_MC.out, dJtheo_nomdumodele_MC.out). Ces fichiers pourront être également utilisés dans une analyse statistique.

4 Discussion

Le programme décrit dans ce manuscrit permet d'analyser des données de relaxation enregistrées à une ou plusieurs intensités de champ magnétique B₀. Il a été écrit en langage MATLAB (version 4.2c.1) sur un ordinateur MacIntosch (MacOS9) : la portabilité de ce langage le rend utilisable sur tout autre type de plateforme et système d'exploitation (MacOS X, Windows, Linux, UNIX). Une version en langage OCTAVE a été testé notamment sur des données de relaxation mesurées à trois intensités de champ magnétique sur la protéine P13^MTCP1 : les résultats de l'analyse dynamique sont présentés sur la Fig. 8. Pour cette protéine relativement globulaire et de symétrie quasi-sphérique, l'approche « globale » (une valeur de τ_c commune à tous les résidus) est pertinente. Notre analyse nous a permis de déterminer un temps de corrélation global de 10,4 ns, compatible avec la taille de la protéine (107 résidus, 13 kDa) et la température à laquelle ont été réalisées les mesures (20 °C). Elle montre également des mouvements internes relativement contraints dans le tonneau β (modèle LS2), alors que la longue boucle, qui relie les deux motifs de feuillet β constituant le tonneau, est animée par des mouvements internes relativement complexes, avec des temps de corrélation proches de la nanoseconde (modèle LSE3, l'utilisation du modèle LSE4 n'apportant pas une amélioration statistiquement significative de l'ajustement). De tels mouvements ont été également identifiés sur l'une des deux boucles « plissées » qui émergent sur une face du tonneau β. Comme le montrait déjà la relative dispersion des valeurs de J(0), des contributions d'échange significatives ont été mesurées pour de nombreux résidus, répartis sur des zones bien particulières de la protéine : une face du tonneau β, autour de l'hélice présente dans la longue boucle, sur l'autre boucle « plissées » émergeant du tonneau β… Des valeurs similaires pour ces contributions d'échange ont été obtenues par correction des valeurs de J(0) [17] ou par ajustement direct des valeurs expérimentales (modèle LSREX).

Si des résultats similaires peuvent être obtenus en analysant les données de relaxation obtenues à une seule fréquence de mesure (résultats non montrés), l'analyse à plusieurs fréquences permet une bien meilleure discrimination entre les différents modèles par les tests statistiques utilisés, et augmente ainsi notre confiance dans les résultats de notre analyse dynamique.

Le programme DYNAMOF est téléchargeable à partir du serveur du Centre de biochimie structurale (ftp://ftp.cbs.cnrs.fr/pub/DYNAMOF_1.0.zip).

Remerciements

Le travail de Virginie Ropars est financé par une bourse de la Ligue contre le cancer. L'étude des protéines oncogéniques de la famille TCL1 constitue un projet du laboratoire soutenu par l'Association pour la recherche contre le cancer. La conversion des données de relaxation en valeurs de densité spectrale utilise une routine MATLAB écrite par J.-F. Lefèvre : Philippe Barthe et Christian Roumestand remercient cet ami, malheureusement trop tôt disparu, qui a su leur communiquer sa passion pour l'étude de la dynamique interne des protéines et la relation entre les propriétés dynamiques des protéines et leur activité.