Modélisation des assemblages de polymères des parois végétales

Dans beaucoup d’applications, il est indispensable d’imposer des conditions thermodynamiques au système. Il faut donc forcer la trajectoire à être compatible avec l’ensemble statistique choisi. Un ensemble est défini par les grandeurs qui y sont conservées : le nombre d’atomes N, le volume V, la pression P et la température T. On peut citer trois ensembles thermodynamiques très courants dans les simulations de dynamique moléculaire : l’ensemble microcanonique NVE ; l’ensemble canonique NVT ; l’ensemble isotherme et isobare NPT. Les simulations MD effectuées dans le cadre de cette thèse sont basées sur l’ensemble canonique NVT.

L’application correcte (en termes de pas de temps) d’un algorithme de Verlet ou Verlet-vitesse à un nombre fixe de particules dans une boîte de simulation de volume fixe permet de remplir les conditions de l’ensemble NVE. En revanche, imposer les conditions d’autres ensembles nécessite l’utilisation de méthodes plus sophistiquées : thermostats pour fixer la température, barostats pour fixer la pression, etc… Seuls des thermostats ont été utilisés dans le cadre de cette thèse.

La température est un paramètre qui mesure l’énergie cinétique instantanée d’un système moléculaire. Celle-ci peut donc être définie à l’aide des vitesses atomiques : ∑ Équation I-17

Où N est le nombre de particules, mi et vi sont respectivement la masse et la vitesse de la particule i, et kB la constante de Boltzmann. Le thermostat est un dispositif servant à maintenir la température constante dans le système moléculaire. Quelques thermostats usuels sont présentés ici.

Le thermostat de Berendsen [Berendsen et al. 1984] fixe l’énergie cinétique du système moléculaire par renormalisation directe des vitesses. Il est extrêmement efficace pour thermaliser les systèmes hors équilibre à la température que l’on souhaite :

Où λ est le facteur de modification des vitesses des particules du système moléculaire, τ la constante de la puissance du thermostat, T0 est la température du thermostat, T est la température instantanée du système, Δt est le pas de temps.

Cependant, une fois le système équilibré, ce couplage ne permet de préserver aucun ensemble canonique car ce thermostat supprime la fluctuation de la température.

Le thermostat d’Andersen [Andersen 1980] est le thermostat le plus simple capable de produire des ensembles canoniques corrects. Il s’agit de sélectionner de manière aléatoire dans chaque pas de temps un certain nombre de particules et de modifier leurs vitesses selon la distribution de Maxwell-Boltzmann à température T (Eq. I-20) :

Avec P la distribution de probabilité de la vitesse, p la quantité de mouvement. Le thermostat Andersen n’impose pas que la température soit modifiée à chaque pas de temps, mais plutôt d’intervenir tous les 1/v pas de temps, où v est la fréquence de couplage avec le réservoir de chaleur.

Le thermostat d’Andersen peut générer des ensembles canoniques corrects. Cependant, les autocorrélations des vitesses sont détruites, ce qui interdit l’usage de ce thermostat pour calculer des coefficients de diffusion.

Le thermostat Nosé-Hoover, introduit par Nosé [Nosé 1984a, Nosé 1984b], puis amélioré par Hoover [Hoover 1985], est un thermostat déterministe produisant une statistique correcte des ensembles canoniques.

En pratique, les équations du mouvement sont modifiées (Eq. I-21 Eq. I-22, Eq. I-23) par un terme semblable à un terme de dissipation visqueuse de viscosité variable :

Où ⃗ est la quantité de mouvement de l’atome i, la masse de l’atome i, Df  est le nombre de degrés de liberté, T0 la température du thermostat et Q un paramètre ajustable, appelé masse virtuelle. Notons que le paramètre ε est identique pour toutes les particules d’un système. Lorsque la température instantanée du système est plus grande que T0, le paramètre augmente et prends une valeur positive. Le terme additionnel est alors équivalent à une dissipation visqueuse, et l’énergie du système diminue. Lorsque la température instantanée du système est moins grande que T0, diminue, et lorsque la valeur de cette variable devient négative, l’énergie du système augmente.

Le paramètre Q fixe la fréquence d’oscillation de la température. Si la valeur de Q est trop grande, les fluctuations de température sont très grandes, en revanche, si la valeur de Q est trop faible, la température est correctement fixée, mais l’échantillonnage des micro-états est insuffisant et la statistique de l’ensemble canonique n’est plus correcte.

Le bois se compose élémentairement de 50% de carbone, 42% d’oxygène, 6% d’hydrogène, 1% d’azote et 1% d’éléments divers. Les constituants structuraux du bois sont la cellulose, les hémicelluloses et la lignine. Les propriétés physiques du bois sont assurées par ces trois composants chimiques principaux. Il existe en plus de ces trois constituants principaux des substances mineures extractibles de nature organique, ainsi que des cendres.

Le bois contient 40% à 60% de cellulose selon l’essence. La cellulose est la molécule organique de polysaccharide la plus abondante sur terre. Elle est la composante principale des microfibrilles qui peuvent être considérées comme des cordes fabriquées de fils de très longues molécules de cellulose. C’est un polymère linéaire d’unité de D-glucose cyclique de deux configurations natives Iα et Iβ. La structure de la cellulose Iα est P1 triclinique avec uniquement une cellobiose dans l’unité répétitive de la chaîne dont les paramètres sont : a = 0,6717 nm, b = 0,59962 nm, c = 1,040 nm, α = 118,08º, β = 114,80º, γ = 80,37º. La structure de la cellulose Iβ est P21 monoclinique avec deux cellobioses dans l’unité répétitive de la chaîne dont les paramètres sont : a = 0,7784 nm, b = 0,8201 nm, c = 1,0380 nm, α = β =90º, γ = 96,5º Voir (Fig. II-4). L.L.
La cellulose est construite par la répétition de l’unité de dimère cellobiose (Fig. I-15) qui est formée par deux monomères de glucose cyclique β, dont l’une est retournée de 180º par rapport à l’autre. La longueur et la largeur de la cellulose dans la nature sont respectivement 100 nm et 20 nm en moyenne, mais sur certaines espèces, la longueur peut atteindre 1000 nm et la largeur 20 nm [Elazzouzi-Hafraoui et al. 2007].
Les macromolécules de cellulose associées forment des microfibrilles aux états cristalline et amorphe (Fig. I-16), qui sont elles-mêmes associées en couches pour former les parois des fibres végétales. Le monomère de la cellulose est le cellobiose. Dans le bois, le degré de polymérisation « DP » est compris entre 5000 et 10000. La structure de la cellulose permet de former des macromolécules grâce aux liaisons hydrogène qui s’établissent entre les groupements hydroxyles de la cellulose et d’autres molécules. Le taux de cristallisation de la cellulose dans le bois est entre 60% et 70% [Trouy-Triboulot and Triboulot 2012]. La figure I-16 expose la coexistence de cellulose cristalline et amorphe dans les parois végétales.

Comme la cellulose les hémicelluloses sont des polysaccharides. La conformation des chaînes moléculaires des hémicelluloses a un caractère amorphe, aléatoire et présente peu de rigidité. Les hémicelluloses sont composés de différents types de monomères de sucres. Les hémicelluloses possèdent des chaînes de longueurs plus courtes que celles de la cellulose, ainsi que des chaînes ramifiés. Il existe deux types de molécules dans la catégorie des hémicelluloses : les pentosanes et les hexosanes. La distinction entre ces deux types de chaînes moléculaires est liée au type du monomère élémentaire considéré (Fig.I-17). Les hexosanes, représentées principalement par des chaînes des glucomannanes, sont construites à partir de monomères à 6 carbones de glucose ou de mannose. Les pentosanes sont constituées de monomères à 5 carbones de xylose, ce qui permet de former les chaînes de xylane.
Autrement dit, les hémicelluloses présentées dans le bois sont principalement des xylanes et des glucomannanes. La composition des hémicelluloses est différente pour des bois feuillus et pour des bois résineux. Les hémicelluloses des bois feuillus sont majoritairement composées de xylane et les hémicelluloses des bois résineux sont composées principalement de glucomannanes.

Avec la cellulose et les hémicelluloses, la lignine est l’un des principaux composants chimiques du bois. Elle est le deuxième biopolymère renouvelable le plus abondant sur la Terre juste après la cellulose. La lignine est chimiquement très différente de la cellulose et des hémicelluloses. La fonctionnalité principale de la lignine est d’apporter à la paroi cellulaire la rigidité mécanique nécessaire lui conférant son imperméabilité à l’eau et sa résistance à décomposition. La lignine ne peut être définie avec précision car sa grande variabilité et la formation de lignine dépend fortement de l’environnement physico-chimique. La lignine est bio synthétisée dans le bois à partir de trois alcools phénoliques : les alcools paracoumarylique, coniférylique et sinapylique (Fig.I-18).
La proportion de ces trois types de monomère varie de manière importante en fonction du végétal, de l’espèce, de l’organe et du tissu. Il n’est donc pas possible de définir précisément la structure chimique de la lignine, car presque toutes les lignines présentent des variations chimiques dues à l’association des monolignols, dérivés indirects d’un acide aminé servant de source pour la biosynthèse de lignines, par différentes liaisons chimiques sans caractère ordonné ni répétitif pour former après réduction, méthylation et déshydrogénation, un polymère amorphe et hydrophobe.

Les extractibles comprennent de nombreux composants pouvant être extraits du bois par des solvants divers. Généralement, les bois tropicaux contiennent plus d’extractibles que le bois venant des zones tempérées. Les éléments extraits sont divers : des tannins, des terpènes, des lignanes, des stilbènes, des flavonoïdes et d’autres composants aromatiques. On retrouve également des lipides, des cires, et des acides gras [Trouy-Triboulot and Triboulot 2012].
Les composants minéraux sont normalement analysés à partir des cendres après combustion totale du bois à 600-800 ºC. Les principaux composés minéraux sont le potassium, le calcium, le magnésium et la silice.

Avec les composants chimiques présentés dans cette section I.3.3, quelques modèles expérimentales de la paroi végétale peuvent être citées. Dans le modèle de [Terashima et al. 2009], la nanostructure des assemblages dans la couche primaire et la couche secondaire de la paroi de cellules a été présentée. Le modèle général [Terashima et al. 1993] de la paroi plus lignifiée a été également présenté, ainsi que les modèles améliorés [Terashima et al. 1995, Terashima et al. 1996, Terashima et al. 1996]. Ces modèles de la sous-structure, existent dans les couches primaire et secondaire de cellules, constituent importantes références pour la construction de modèles numériques du réseau moléculaire de la paroi.

Le champ de force utilisé pour calculer la conformation des xylanes est Molecular Mechanic 3 ou MM3 [Allinger et al. 1989, Allinger et al. 1990]. Ce champ de force a été utilisé pour le calcul de la conformation élémentaire des dimères de xylose et les hélices des chaînes de xylanes.
Le modèle atomistique de xylane est composé par l’unité répétitive de dimère D-xylopyranose avec la liaison 1-4β. (Fig. II-1) Seule la configuration 4C1 a été considérée pour la construction des cycles. Chaque monomère a été construit à partir des paramètres standards issus de la cristallographie au niveau de la longueur de liaisons, angles interne des liaisons.
La forme géométrique d’un disaccharide est majoritairement dominée par la liaison glycosidique [French et al. 1993]. Afin de construire la configuration d’une chaîne de xylane, il est nécessaire d’établir la conformation autour de la liaison glycosidique entre deux monomères. L’orientation relative de deux unités saccharides peut être exprimée à l’aide de deux paramètres portant respectivement sur les deux angles de torsion autour de la liaison glycosidique : Φ = O5’-C1’-O1’-C4 et Ψ = C5-C4-O1’-C1’.
Les hélices de la chaîne de xylane ont été générées par le programme POLYS, abréviation de « Polysaccharide » [Engelsen et al. 1996]. POLYS nécessite en entrée les géométries de différentes unités de monosaccharide et les angles de torsion Φ et Ψ. Le monomère provient de l’étude préliminaire du disaccharide et a été implanté dans le MONOBANK [Pérez and Delage 1991].
La génération des conformations de chaînes aléatoires contenant un grand nombre de monomères a été effectuée en utilisant le programme METROPOL [Boutherin et al. 1997].
La chaine de xylane est formée à partir d’un disaccharide initial, en ajoutant un monomère suivant à l’aide de valeurs Φ et Ψ aléatoires. Le nouveau résidu, décrit par Φ et Ψ, sera accepté si l’énergie de la nouvelle conformation est inférieure ou égale à celle de la conformation précédente. Sinon, un test de « Metropolis » intervient : on génère tout d’abord un nombre aléatoire entre 0 et 1, ensuite, on compare ce nombre avec le facteur de Boltzmann e-ΔE/kBT, où ΔE est l’écart d’énergie, kB est la constante de Boltzmann et T est la température.
Si le nombre généré est inférieur au facteur, la nouvelle conformation sera acceptée, sinon, Φ et Ψ se voient affecté de nouvelles valeurs aléatoires, la génération de la chaine recommencera jusqu’à ce que la nouvelle conformation soit acceptée.
Ensuite, la structure générée fait l’objet d’une thermalisation par dynamique moléculaire. Suit la génération de quatre trajectoires dans l’ensemble NVT sous la température 300 K. Les modèles atomistiques concernant le xylane [Mazeau et al. 2005] contiennent 3 extractions de trajectoire de 500 pas de temps chacune sous forme de fichier PDB « Protein Data Bank ». Chaque trajectoire ne comporte qu’une chaine de xylane. Les degrés de polymérisation (DP) dans ces trois trajectoires (Fig. II-2) sont respectivement 5, 10, 25. La représentation visuelle des modèles atomistiques sont illustrée dans la figure II-2.

Pour modéliser l’interaction intermoléculaire xylane-xylane, un extrait de trajectoire de 500 pas a été utilisé. Ce modèle atomistique est composé de deux chaines de xylane de 5 xyloses chaque. Les deux chaînes, ont été comprimées de façon à les faire interagir. La simulation atomistique a ensuite été exécutée sous la température de 300 K. Le modèle généré de l’interaction entre deux chaines moléculaire du xylane est présenté dans la figure II-3.

Pour modéliser l’interaction intermoléculaire xylane-cellulose, un extrait de trajectoire de 500 pas a été utilisé. Ce modèle atomistique est composé d’une chaine de xylane de 5 xyloses, et de 3 couches de cellulose de 8 chaines chacune, chaque chaine comportant 8 glucoses.
La surface (110) est l’une des deux surfaces dominantes, (110) et (1-10), dans les microfibrilles de cellulose [Mazeau et al. 2012, Van Daele et al. 1992, Hanley et al. 1992, Helbert et al. 1998a, Helbert et al. 1998b, Revol 1982, Sassi and Chanzy 1995, Sassi et al. 2000, Sugiyama et al. 1985]. De plus, il a été constaté par Résonnance magnétique nucléaire (RMN, ou NMR en anglais) que les chaînes du xylane sont préférentiellement en interaction avec la surface (110), beaucoup plus qu’avec la surface (1-10) [Bergenstråhle et al. 2008, Larsson et al. 2004].
Dans le cadre de la thèse, seule une trajectoire atomistique de l’absorption d’une chaîne de xylane de DP = 5 sur la surface (110) de la cellulose a été utilisée pour la modélisation de potentiel gros grain de xylane-cellulose, les deux autres surfaces (1-10) et (100) sont également prises en compte, les intensités des potentiels d’interaction xylane-cellulose étant ajustées à l’aide du ratio entre l’énergie d’adsorption du xylane sur ces surfaces et cette même énergie dans le cas de la surface (110), ratio obtenu, dans chaque cas, par dynamique moléculaire atomistique.
Le cristal de cellulose a été construit à partir des paramètres de la structure Iβ cristal, illustré dans la figure II-4, déterminé par diffraction des Rayon X et diffusion de Neutrons [Nishiyama et al. 2002], dont les paramètres de cristal sont a = 0,77848 nm, b = 0,82018 nm, c = 1,38010 nm, α = β = 90º, γ = 96,5º, les données correspondantes à la structure géométrique sont illustrées dans la figure II-4. La cellule a été dupliquée dans l’espace pour former 8 couches de celluloses de 8 chaines chacune.
Les 8 couches ont été équilibrées dans l’ensemble NPT, sous une température de 300K et une pression de 1 atm [Mazeau 2005]. Ensuite, le paramètre de la cellule de direction b a été manuellement augmenté afin de créer une membrane de cellulose exposant deux surfaces équivalentes et parallèles à (110). Le modèle atomistique final ne contient que 3 couches de cellulose qui sont déployées dans une boîte de simulation avec PBC dans toutes les directions après le retrait de 5 couches de cellulose. Ce modèle représente une surface accessible des microfibrilles natives de cellulose. Seuls les groupes fonctionnels de la couche externe de la cellulose sont actifs. Les coordonnées du reste des atomes de la cellulose sont fixées.
Une chaîne de xylane de 5 monomères a été générée selon la description résumée dans la section II.1.1.1, elle est ensuite mise au-dessus de la surface (110) de la cellulose avec une orientation aléatoire. Une trajectoire de ce système xylane-cellulose a été ensuite générée [Mazeau and Charlier 2012] dans l’ensemble canonique NVT. La représentation visuelle du modèle xylane-cellulose est présentée dans la figure II-5.

Le champ de forces UFF [Rappé et al. 1992], abréviation des « Universal Force Field », a été utilisé pour les simulations de dynamique moléculaire atomistiques. L’UFF a été utilisé auparavant [Hanus and Mazeau 2006, Henao and Mazeau 2008] pour la description des propriétés des polysaccharides complexes. La distribution des charges électrostatiques provient de la méthode d’équilibration des charges pour la dynamique moléculaire [Rappé and Goddard 1991].
Pour les interactions intermoléculaires, deux types de modèles du potentiel sont pris en compte : la force Van der Waals et la force électrostatique.