Simulations de dynamique moléculaire ab initio

Simulations de dynamique moléculaire ab
initio

Durant ces dernières années, de nombreuses approches théoriques et expérimentales ont permis d’explorer la solvatation d’ions lanthanides, mais très peu d’études ont porté sur des complexes de gadolinium utilisés en imagerie médicale (IRM). Yazyev et al. ont étudié par dynamique moléculaire ab initio de type Car-Parrinello la structure et la dynamique de l’ion gadolinium dans l’eau. Cette étude a mis en évidence une distance moyenne Gd-O de 2.37 Å, ce qui est tout à fait en accord avec les valeurs obtenues lors d’expériences de double résonance électronique et nucléaire (ENDOR) (2.4 ± 0.5 Å), et une distance moyenne Gd-H de 3.04 Å, pour une valeur issue des expériences ENDOR de l’ordre de 3.1 ± 0.1 Å. La distance moyenne Gd-H estimée par Yazyev et al. est aussi en accord avec les observations expérimentales de Caravan   sur des agents de contraste, celui-ci estime la distance moyenne Gd-H à 3.1 ± 0.1 Å. Les résultats issus de l’étude de Yazyev par dynamique moléculaire ab initio sont proches de ceux obtenus par Ikeda et al. par métadynamique biaisée sur l’ion diamagnétique Y+3 (où une distance moyenne Y-O avait été estimée à 2.38 Å). Les ions yttrium et gadolinium ont les mêmes charges et des rayons ioniques similaires , leur comportement en solvatation dans l’eau peut donc être comparé. Malgré ces premières études sur l’ion libre, peu d’études de dynamique (classique ou ab initio) ont été menées sur des complexes de gadolinium utilisés en tant qu’agents de contraste en IRM. Pourtant l’ion libre est très toxique, et il est donc nécessaire de le complexer dans le cas d’une utilisation médicale. De plus, l’agent étant solvaté dans une boîte d’eau, il est possible d’avoir accès à des contributions de la relaxation dites de sphère interne et de sphère externe, ce qui constitue un avantage pour la recherche dans le domaine [149]. Malgré les espoirs apportés par la méthode ab initio dans l’étude de systèmes biologiques [150, 151], l’augmentation de taille du système pose de nouvelles difficultées pour sa modélisation. Yazyev et al. [37] ont étudié (par une approche de dynamique classique) un agent de contraste commercial connu sous le nom de Dotarem qui est un complexe de gadolinium de ligand chélatant DOTA (voir Chapitre 1). Les distances moyennes Gd-O et Gd-H obtenues par simulation sont respectivement de 2.56 ± 0.06 Å et 3.27 ± 0.14 Å. Nous comparerons ces distances, qui ne sont pas très éloignées de celles que nous avons obtenues, dans le cas d’un autre ligand chélatant. Dans ce 69 Simulations de dynamique moléculaire ab initio type de système, il est important de tenir compte de la polarisation, comme l’ont fait Clavaguéra et al. [152] ou Villa et al. [153], qui ont respectivement obtenu par dynamique moléculaire classique, en utilisant un champ de force polarisable, une distance Gd-O de 2.44 Å et de 2.55 Å. Dans ce chapitre, le recours à la dynamique classique ne permet pas d’avoir accès à des informations relatives à la structure électronique du métal, aux liaisons faibles se formant avec son environnement. Comme nous l’avons vu dans le chapitre précédent, afin d’avoir accès à ces différentes propriétés électroniques, nous avons choisi d’avoir recours à la dynamique moléculaire ab initio pour suivre le comportement de complexes de gadolinium d’intérêt pour l’IRM. Je rappelle ici que les dynamiques ont été réalisées à l’aide du logiciel CPMD , ce qui implique l’utilisation de pseudopotentiels, standards pour les atomes présents dans la boîte de simulation. La description du gadolinium pose quelques questions quant à la définition électronique du coeur et de la valence. Les électrons 4f des lanthanides sont généralement inclus dans le coeur des pseudopotentiels du fait de leur forte localisation dans la région du noyau [154, 155]. Un des membres de notre groupe de recherche a montré que les électrons 4f pouvaient être insérés dans le coeur du pseudopotentiel, mais en gardant bien en mémoire que certains problèmes de structures peuvent en découler [32]. Ces résultats avaient conduit à l’implémentation de pseudopotentiels de type Vanderbilt dans CPMD, et nous avons donc pu en profiter. La construction du pseudopotentiel du gadolinium est abordée dans le Chapitre 2. Dans ce chapitre, la procédure mise en place afin d’équilibrer les dynamiques que j’ai simulées sera détaillée. Un agent de contraste commercial, le ProHance sera examiné, ainsi que deux dérivés de ce système. Chacune des dynamiques sera étudiée, et des propriétés structurales et dynamiques de ces trois systèmes seront discutées et comparées. Les différentes dynamiques moléculaires ab initio présentées dans la suite ont été réalisées à l’aide de la fonctionnelle PBE, et de pseudopotentiels de type Vanderbilt. Les autres détails de calculs seront précisés dans chacun des cas. 

Dynamique moléculaire ab initio d’un agent de contraste solvaté 

La première dynamique moléculaire ab initio (DMAI) qui a été réalisée par notre groupe concerne l’agent ProHance solvaté [34, 156]. Pour cela, une structure cristallographique de l’agent a été utilisée afin que la géométrie de départ de la dynamique concorde avec la géométrie expérimentale du ProHance. La géométrie cristalline récupérée a été solvatée dans une boîte d’eau, et une dynamique classique a été réalisée à l’aide du logiciel GROMACS afin de relaxer le système, et d’obtenir une première répartition des couches de solvatation. Dans un premier temps la dynamique classique a été effectuée dans l’ensemble canonique NPT afin d’identifier le volume à l’équilibre du système. Puis la simulation s’est poursuivie dans l’ensemble NVT au  volume choisi. Cette première étape permet d’obtenir à moindre coût une géométrie équilibrée de l’agent de contraste ProHance solvaté. C’est à partir de cette configuration que la géométrie de départ de la DMAI a été construite. Nous avons vérifié qu’une molécule d’eau était bien coordinée au centre métallique gadolinium dès le début de la DMAI, et 98 autres molécules d’eau sont contenues dans la boîte d’eau. Le système ainsi construit a permis de réaliser une dynamique moléculaire ab initio, dont l’équilibre a été atteint et contrôlé au cours de différentes étapes qui seront détaillées par la suite. Pour déterminer si le système est à l’équilibre thermodynamique, il s’agit de suivre les fluctuations de l’énergie potentielle du système (énergie de Kohn-Sham), mais aussi les fluctuations des distances dans le complexe. Enfin, la solvatation des groupements carboxylates constitue aussi un indice de l’atteinte de l’équilibre thermodynamique. Deux autres dynamiques moléculaires ab initio ont été obtenues. Les géométries de départ ont été construites à partir d’une configuration équilibrée du ProHance, et par changement de la molécule coordinante ou modification du ligand. Des précautions particulières ont été prises afin de relaxer de manière progressive ces nouveaux systèmes, elles seront détaillées dans les paragraphes correspondants.

Etude du ProHance : Gd(HP-DO3A) 

Description du système 

L’agent de contraste ProHance est un complexe de gadolinium, dont le ligand est un dérivé du cyclène [157] comportant trois bras acétate et un bras hydroxy-propyl. La géométrie est de type anti-prisme à base carrée (voir Figure 3.1). Figure 3.1 – Géométrie du ProHance : anti-prisme à base carrée Il existe différents isomères du ProHance dont les géométries diffèrent par le positionnement des bras (voir Figure 3.2 ∗ ). Dans le cas du ligand DOTA (voir Chapitre 1), des études expérimentales ont permis d’identifier l’isomère majoritaire en solution [158]. Pour notre ligand HP-DO3A, nous nous sommes intéressés à l’analogue de l’isomère majoritaire du DOTA : l’isomère que nous appellerons SA. ∗. Figure originale issue de l’article de Caravan [76] 71 Chapitre 3 : Simulations de dynamique moléculaire ab initio Figure 3.2 – Différents isomères du ligand type HP-DO3A (SA et TSA). Ce chélate se caractérise par huit sites de coordination sur le gadolinium, laissant un site vacant pour une molécule d’eau coordinante. En effet, les 4 atomes d’azote du complexe, ainsi que les 4 atomes d’oxygène appartenant aux groupements hydroxyls et acétates, sont chélatants (Figure 3.3). Cet agent de contraste est neutre, la triple ionisation du gadolinium étant contrée par les charges négatives des trois oxygènes des groupements acétates. Dans la suite de notre étude, nous nous intéresserons à la relaxivité du ProHance qui dépend de l’isomère considéré. Dans notre cas, les temps nécessaires pour passer d’un isomère à un autre sont au minimum de l’ordre de la dizaine de millisecondes (10 à 50 ms selon l’isomérisation). Nous pouvons donc étudier la relaxivité d’un unique complexe de manière approfondie, sans risque d’observer au cours de la dynamique une isomérisation, étant donnés les temps accessibles en DMAI. La molécule d’eau coordinée possède un temps de résidence relativement long (de l’ordre de 220 ns à 310 K [2]), ce qui nous permet également, à l’échelle de la dynamique ab initio, d’échantillonner les positions dans l’espace d’une molécule d’eau coordinée, sans observer d’échange avec l’environnement. 

Equilibration de la dynamique

 La géométrie de départ du complexe est tirée d’une structure cristalline à partir de laquelle la dynamique classique a été lancée (voir Table 3.1). Pour la partie DMAI, le système a directement été relaxé à partir d’une géométrie issue de la dynamique classique, sans contrainte particulière sur le ligand ou la molécule d’eau coordinée [34]. Afin de vérifier l’équilibration du système, l’énergie de Kohn-Sham est observée le long de la dynamique. Les temps d’équilibration nécessaires pour chaque dynamique sont détaillés dans 72 3.2 Etude du ProHance : Gd(HP-DO3A) Figure 3.3 – Agent de contraste ProHance et sa molécule d’eau coordinée, à gauche vue de dessus et à droite vue de côté. Les atomes de carbone sont représentés en noir, ceux d’oxygène en rouge, les atomes d’azote en bleu, et les hydrogènes en blanc. le tableau 3.1 et une géométrie équilibrée du système solvaté est représentée en figure 3.4. Au cours de la DMAI, différents pas de temps, différentes masses fictives, ainsi que quelques autres options de calcul ont été utilisées, afin de contrôler l’adiabaticité du système et d’établir les conditions permettant d’optimiser la précision et le temps de calcul. Ainsi, pour cette dynamique, une fois le système équilibré, la propagation a été effectuée avec une masse fictive de 700 u.a et un pas de temps de 6 u.a ( ≈ 0.15 fs). Ces deux paramètres sont importants pour une dynamique de type Car-Parrinello (CP), le but étant de rester au plus près de la surface de Born-Oppenheimer (BO). Ces paramètres ont été correctement choisis puisque la conservation de l’énergie cinétique fictive des électrons a été vérifiée, ce qui indique que le système reste proche de sa surface de Born-Oppenheimer. J’ai également utilisé l’astuce de calcul qui consiste à remplacer les hydrogènes par des deutériums (voir Chapitre 2), ce qui nous a permis d’utiliser une grande masse fictive et un grand pas de temps. Nous avons donc pu avoir une dynamique rapide et stable du point de vue de la surface de Born-Oppenheimer. La simulation de dynamique moléculaire a été réalisée avec un thermostat de Nosé pour les électrons et les ions, ce qui a permis de contrôler la température thermodynamique du système à 300 K pour une fréquence du thermostat de 2500 cm−1 . L’énergie cinétique fictive des électrons a été fixée à 0.03255 u.a pour une fréquence de 15000 cm−1 , permettant d’éviter un transfert d’énergie entre les sous-systèmes. La boîte de simulation est périodique cubique, de côté 15.40 Å. L’énergie de coupure est de 30 Rydberg, ce qui est faible pour ce type de système (l’énergie de coupure pour un système carboné étant en général d’environ 90 Rydberg), grâce à l’utilisation de pseudopotentiels de Vanderbilt.