Inférence de paramètres de modèles de Thomas par vérification de propriétés LTL 

Télécharger le fichier pdf d’un mémoire de fin d’études

État de l’art de la modélisation de GRN

Depuis les années 60, de nombreux formalismes mathématiques ou informatiques ont été proposés pour modéliser les réseaux de régulation génétique : sous la forme d’équations différentielles ordinaires, partielles, de réseaux bayésiens ou sous d’autres formes stochastiques, de réseaux de Petri, de réseaux discrets booléens ou multivalués ou encore d’automates hybrides. Dans un premier temps, nous esquisserons la modélisation sous la forme d’équations différentielles
ordinaires, avant de nous intéresser à la modélisation discrète, et en particulier au modèle dit de René Thomas. Pour un aperçu plus complet de l’état de l’art sur la modélisation des GRN, le lecteur pourra par exemple consulter [Jon02 ; SBB00].

Modélisation par système d’ODE

Classiquement, la modélisation mathématique des GRN se base sur des systèmes d’équations différentielles ordinaires (ODE) décrivant la variation des concentrations des différentes espèces en présence (voir par exemple [Goo65 ; C+99]).
Par nature, ces modèles constituent une approche quantitative des réactions sous-jacentes aux phénomènes de régulation entre gènes et protéines régulatrices. Dans ce cadre, la régulation d’un gène est modélisée par une équation de vitesse de réaction, exprimant la vitesse de synthèse de la protéine associée au gène en fonction de la concentration de ses régulateurs.
La non-linéarité des ODE4 ne permet généralement pas une résolution analytique de ces systèmes. Une alternative possible est de recourir à des techniques de résolution numérique (cf. [MS95 ; SA85] par exemple pour l’étude du phage ). Cependant, les simulations nécessitent de connaître précisément les valeurs des constantes qui apparaissent dans les ODE, correspondant aux paramètres cinétiques du système (notamment les constantes de production et de dégradation des équations des vitesses de réaction). Or ces paramètres sont difficilement évaluables à cause du manque de données expérimentales in vivo ou in vitro.

Modélisation discrète booléenne

En 1969, Kauffman propose dans [Kau69] une modélisation qualitative s’appuyant sur une approche discrète (vis-à-vis du temps et des états). Dans ce formalisme, les GRN sont représentés par un modèle booléen (dit aussi modèle binaire) dans lequel l’expression des gènes est activée ou inactivée par le biais d’un switch binaire on/off enclenché par la présence ou l’absence des protéines régulatrices correspondantes (elles-mêmes étant représentées par des gènes activés ou inactivés). Chaque gène est ainsi associé à une variable qui peut prendre la valeur 0 (pour la valeur booléenne faux, le gène est considéré inactivé et les protéines correspondantes sont absentes) ou 1 (pour vrai, gène activé et protéines présentes).
Une fonction booléenne permet de connaître de manière déterministe le prochain état atteint en fonction des gènes activés dans l’état courant. Cette fonction booléenne est à déterminer : soit au hasard, soit par tâtonnement en testant toutes les possibilités une à une, en utilisant par exemples des connaissances sur le comportement du GRN.
La transition entre les états est synchrone, i.e. plusieurs gènes peuvent évoluer simultanément. L’évolution au cours du temps des valeurs des variables associées aux gènes est appelée dynamique du GRN et caractérise son comportement.
Exemple 6. Pour le GRN introduit précédemment (figure 2.2), on peut imaginer le modèle dont les entrées et sorties de la fonction booléenne sont décrites dans le tableau 2.1.
Dans ce modèle, l’état correspondant à = 0 et = 0 (i.e. les gènes et sont inactifs) sera suivi de l’état dans lesquel = 1 et = 0 (i.e. est actif tandis que reste inactif).
La modélisation booléenne des GRN a été étudiée par d’autres auteurs, citons le travail antérieur de Sugita dans [Sug61 ; Sug63 ; SF63] et celui postérieur de Thomas dans [Tho73 ; Tho79] par exemple.

Modélisation discrète multivaluée

En 1990, Thomas et al. proposent une modélisation discrètemultivaluée des GRN dans [Td90] (voir aussi [Tho91 ; TCT93 ; TTK95 ; TT95]), connue sous le nom de modèle de Thomas (ou méthode logique généralisée). Thomas y introduit la notion de seuils dans les modèles discrets, i.e. une protéine ne peut réguler un gène cible que si elle est présente en quantité suffisante, cette quantité étant spécifique d’une régulation donnée.
Par rapport au modèle booléen, la modélisation multivaluée présente l’avantage de prendre en compte de manière plus fine le phénomène de régulation. Par exemple, dans le cas du GRN précédent, les protéines peuvent être en quantité suffisante pour réguler le gène mais insuffisante pour activer le gène : cette nuance n’est possible que dans un modèle multivalué.
La notion de seuil provient de la discrétisation des concentrations du modèle continu des ODE5. Elle est liée à la forme sigmoïdale des courbes V représentant la vitesse de synthèse d’une protéine en fonction de la concentration d’un de ses régulateurs.
En effet, expérimentalement, les biologistes déterminent qu’une protéine p est activatrice d’un gène g si en augmentant progressivement la concentration de p, ils observent une valeur à partir de laquelle la synthèse des protéines de g augmente. On considère généralement que la vitesse de synthèse des protéines de g sature à partir d’un certain seuil de concentration de p, i.e. que la vitesse de synthèse est constante au-dessus d’un seuil spécifique à la régulation de g par p [JM61 ; YY71 ; Tho98].
À l’inverse, un inhibiteur de g est expérimentalement caractérisé par une diminution du taux de synthèse de g à partir d’une valeur de concentration donnée de cet inhibiteur. La vitesse maximale est dans ce cas atteinte pour une valeur de l’inhibiteur inférieure à ce seuil. Un phénomène de saturation se produit également dans le cas d’une inhibition : quand la concentration de p est supérieure à la valeur de seuil, la vitesse de synthèse ne diminue plus.
Ces observations biologiques expliquent le profil des courbes V en sigmoïde : elles sont croissantes ou décroissantes suivant le type de la régulation associée (activation ou inhibition), la régulation est faible tant qu’une concentration minimale spécifique n’est pas atteinte, et pour une concentration bien supérieure à ce seuil, la régulation atteint un palier maximal (cas d’activation) ou minimal (cas d’inhibition).
Pour illustrer notre propos, considérons trois protéines p1, p2 et p3 telles que p1 est un activateur du gène qui synthétise p2 et un inhibiteur du gène qui synthétise p3. Notons respectivement Vp2=p1 et Vp3=p1 les courbes caractérisant la vitesse de synthèse de p2 et p3 en fonction de la concentration de p1. Une illustration possible du profil de ces deux courbes est esquissée en figure 2.3.
Pour une courbe donnée, on néglige l’effet régulatoire potentiel de p1 en-dessous du seuil de régulation, tandis qu’au-dessus, on considère que le plateau de saturation est atteint. Les deux paliers de la fonction en escalier obtenue à partir de V correspondent ainsi à deux états logiques différents : un état dans lequel il n’y pas d’effet de régulation, et un état où la régulation est effective.
Dans l’interprétation logique, la concentration de p est de plus discrétisée : on parle alors de niveaux d’expression du gène associés à p. Dans la figure 2.4, l’abscisse t correspondant au saut entre les paliers est un entier naturel découlant de la discrétisation de la concentration de p : c’est une valeur relative respectant l’ordonnancement des seuils de concentrations de l’ensemble des régulations de p. Ainsi, dans notre exemple, on note t1 et t2 les seuils discrétisés obtenus à partir de c1 et c2 ; comme c1 < c2, alors t1 < t2.
Les seuils de régulation considérés dans le modèle de Thomas correspondent aux seuils t discrétisés de l’interprétation logique des courbes de régulation. Du point de vue des gènes, ils correspondent à un niveau d’expression minimal qu’un gène doit atteindre pour que les protéines qu’il synthétise régulent un ou des gènes cibles du GRN.
Le nombre de seuils de régulation différents pour une protéine, et donc la valeur maximale du niveau d’expression du gène correspondant, est au plus égal au nombre de gènes qu’elle régule (en général il est égal au nombre de gènes régulés, i.e. chaque régulation correspond à un seuil différent).
Dans le cadre de notre exemple jouet, le niveau maximum de est au plus 2 car il régule deux gènes, et (il peut cependant être fixé à 1 si on considère que les seuils de régulation de et sont identiques, ce qui n’est généralement pas le cas, à l’exception des gènes appartenant à un même opéron) ; par contre ne régule qu’un seul gène, , et donc son niveau maximum est 1.
Le plus petit palier, 0, correspond à l’absence de la protéine associée au gène, ou à une concentration trop faible pour réguler un gène du GRN.
Cette approche est qualitative, chaque niveau d’expression est lié à un effet de régulation différent. Par exemple, on peut imaginer que le gène puisse atteindre trois niveaux d’expression différents 0, 1 et 2, caractérisés de la manière suivante :
— 0 signifie que est inactif et qu’il ne régule aucun gène,
— 1 correspond au niveau à partir duquel peut réguler un gène (le gène par exemple),
— 2 correspond au niveau à partir duquel peut réguler un autre gène ().
Le fonctionnement d’un GRN est caractérisé par sa dynamique, autrement dit par l’évolution au cours du temps des concentrations des protéines, ou des niveaux d’expression des gènes. Ces dynamiques sont assujetties à la connaissance de paramètres biologiques indicatifs des effets de compétition entre protéines activatrices et inhibitrices d’un même gène.
Par exemple, quand et ont tous les deux un niveau d’expression supérieur à leur seuil respectif de régulation sur , on ignore à priori si l’effet d’activation de sur est plus important que l’effet d’inhibition de sur et donc si le niveau d’expression de va avoir tendance à augmenter, à diminuer ou à rester stable.
Chaque dynamique peut être caractérisée par un ensemble de paramètres, dits paramètres du modèle de Thomas (ou paramètres logiques) et définis comme les points focaux (constants) représentant le niveau vers lequel tend un gène en fonction des ressources disponibles (i.e. les régulateurs qui sont au-dessus de leur seuil de régulation, que nous appellerons régulateurs effectifs dans la suite).
Dans le cadre discret, l’évolution au cours du temps d’un GRN est caractérisée par une succession d’états, chacun étant déterminé par l’ensemble de valeurs des niveaux d’expression des gènes du GRN au temps t discret correspondant. Il existe deux approches antagonistes pour modéliser l’écoulement du temps : l’approche synchrone et l’approche asynchrone.
Dans le cas synchrone, l’évolution est déterministe et plusieurs niveaux d’expression peuvent évoluer en même temps de manière à atteindre immédiatement l’état ciblé. Ce cadre est notamment utilisé dans la modélisation booléenne de Kauffman [Kau69] (cf. section 2.2.2).
Dans le cas asynchrone, un seul niveau d’expression peut évoluer à la fois, en considérant que chaque niveau d’expression tend à évoluer vers la valeur qu’il aurait dans l’état ciblé par l’approche synchrone. Une transition synchrone définit de fait plusieurs transitions asynchrones : l’état source de ces transitions asynchrones est le même que celui de la transition synchrone, mais par contre chaque transition asynchrone correspond à une augmentation ou une diminution d’un seul niveau d’expression d’une seule unité.
Alors que lors d’une transition synchrone le point focal est nécessairement atteint, cela n’est plus le cas lorsque la transition est rendue asynchrone. En effet, les paramètres caractérisant la prochaine évolution dépendent des niveaux d’expressions des gènes dans l’état courant : le point focal évolue avec chaque transition asynchrone suivie, car l’état évolue lui aussi. Par conséquent, l’approche asynchrone apporte un niveau de complexité supplémentaire pour la détermination des dynamiques par rapport à l’approche synchrone.
L’approche asynchrone présente deux avantages par rapport à l’approche synchrone. D’une part, elle permet de considérer qu’il peut exister plusieurs choix d’états atteignables à partir de l’état courant, ce qui relâche une contrainte forte. D’autre part, elle permet également de s’affranchir de la possibilité de processus de régulation simultanés, qui sont peu probables en pratique. En effet, cela suppose la simultanéité du changement de niveau d’expression de plusieurs gènes à priori distants, qui reposent elle-même sur la simultanéité des mécanismes sous-jacents, i.e. la régulation de plusieurs gènes d’une part et la variation significative des concentrations des protéines correspondantes d’autre part. C’est pour ces raisons que l’approche asynchrone est privilégiée par Thomas et al. [Td90].
Enfin, troisième possibilité, on peut également se placer dans un cadre à la fois synchrone et asynchrone comme Gonzalez et al. ([Gon+06 ; Fau+06]) qui déclinent le modèle de Thomas en autorisant des transitions synchrones en plus des transitions asynchrones.
Dans le cas asynchrone, il peut exister plusieurs possibilités d’évolution à partir de l’état courant. Signalons que plusieurs choix peuvent être pris en conséquence :
— considérer que toutes les transitions se valent et donc exploiter tous les cas possibles (cas de base du modèle de Thomas [Td90]) ;
— ajouter la notion de délai ou de priorité aux transitions (dans [Ahm+09] par exemple, Ahmad et al. proposent un modèle hybride en ajoutant des contraintes temporelles au modèle de Thomas pour simuler l’évolution continue des niveaux d’expression) ;
— se placer dans le cadre probabiliste, en attachant une probabilité de franchissement aux transitions.
Ces deux derniers choix permettent de brider le non-déterminisme dû à l’asynchronicité.
Les paramètres de Thomas capturent de manière plus abstraite les paramètres des ODE.
Il est admis (ou pour certaines propriétés, démontré) que les dynamiques discrètes capturent l’essentiel des caractéristiques qualitatives des dynamiques continues. Thomas a par exemple conjecturé dans [Tho81] deux propriétés importantes des GRN, portant sur le lien entre les circuits de rétroaction positifs (i.e. qui contiennent un nombre pair d’inhibitions) d’une part ou négatifs (i.e. qui contiennent un nombre impair d’inhibitions) d’autre part, présent dans un graphe d’interactions (qui représente les régulations entre les gènes du GRN, nous le définirons formellement dans la suite, cf. définition 2.3.1) et l’existence respective de plusieurs états stables ou d’oscillations entretenues. Ces conjonctures ont été démontrées au fil des ans dans plusieurs cadres différents : tout d’abord, continu (dans [PMO95; Sno98 ; Sou03 ; Sou06] par exemple), puis discret booléen [RRT08], et enfin discret généralisé [Ric06 ; RC07 ; Ric09 ; Ric10b].
Qu’il s’agisse de modélisation se basant sur des ODE ou de modélisation discrète, la question centrale est celle de la détermination des paramètres du modèle, qui déterminent précisément le comportement des dynamiques issues du GRN. Le modèle de Thomas a donné lieu à de nombreux travaux de modélisation, appliqués par exemple à la production de mucus de Pseudomonas aeruginosa [Ber+04], au cycle de vie du phage [TT95], à la morphogenèse d’Arabidopsis thaliana [MTA99], ou encore à la régulation de la réponse immunitaire [Mur+96]). Il a conduit en particulier à la mise au point de méthodes outillées d’aide à la détermination des paramètres (présentées en détail dans la section 2.4), cadre dans lequel nous nous plaçons également.

Le modèle multivalué de R. Thomas

Dans cette section, nous détaillons le cadre bien établi de l’approchemultivaluée de R. Thomas présentée en section 2.2. En particulier, nous introduisons ou redéfinissons des notations propres à cette thèse (proches de celles de [Ber+04] ou [Bar+12]), afin de faciliter la mise en oeuvre de notre approche dans la suite.
Remarque. Par souci de simplification, nous assimilons gène et protéine en considérant qu’un gène est associé à un seul type de protéines produites.

Graphe d’interactions

Un GRN est classiquement représenté par un graphe d’interactions (abrégé en IG pour Interaction Graph), qui est un graphe orienté et étiqueté dans lequel les sommets sont associés aux gènes du GRN et les arcs aux interactions (de régulation) entre protéines et gènes (définition 2.3.1).

Table des matières

Introduction 
I Inférence de paramètres de modèles de Thomas par vérification de propriétés LTL 
1 Préliminaires
1.1 Autour des ensembles
1.2 Contraintes sur les entiers
1.3 Éléments de la théorie des langages
1.4 Systèmes de transition
1.5 Logique Temporelle Linéaire
2 Présentation des réseaux de régulation génétique
2.1 Contexte : les réseaux de régulation génétique
2.2 État de l’art de la modélisation de GRN
2.3 Le modèle multivalué de R. Thomas
2.4 Sélection de modèles à partir d’observations temporelles
3 Méthode d’inférence de paramètres de Thomas basée sur le model-checking LTL et la résolution de contraintes
3.1 GRN Paramétré
3.2 Synchronisation d’un PGRN avec un automate de Büchi
3.3 Exécution symbolique du Produit
3.4 Inférence des paramètres
4 Implémentation
4.1 SPuTNIk
4.2 Algorithme de simplification du produit
4.3 Algorithme d’exécution des arbres
4.4 Pistes d’optimisation
5 Etudes de cas
5.1 Inductibilité de la cytotoxicité de Pseudomonas aeruginosa
5.2 Cycle de vie du phage
II Analyse de la voie de signalisation Wnt/-caténine à l’aide du modelchecking statistique HASL 
2.1 Contexte : les réseaux de régulation génétique
2.2 État de l’art de la modélisation de GRN
2.2.1 Modélisation par système d’ODE
2.2.2 Modélisation discrète booléenne
2.2.3 Modélisation discrète multivaluée
2.3 Le modèle multivalué de R. Thomas
2.3.1 Graphe d’interactions
2.3.2 Espace des états
2.3.3 Dynamique
2.3.4 Paramètres du modèle de Thomas
2.3.5 Contraintes sur les paramètres
2.4 Sélection de modèles à partir d’observations temporelles
3.1 GRN Paramétré
3.1.1 Caractérisation des états par contraintes
3.1.2 GRN Paramétré
3.1.3 Dynamiques associées à un PGRN
3.2 Synchronisation d’un PGRN avec un automate de Büchi
3.3 Exécution symbolique du Produit
3.3.1 Vue d’ensemble
3.3.2 Construction des Arbres d’Exécution
3.3.3 Recherche des cycles acceptants
3.4 Inférence des paramètres
Conclusion 
Bibliographie 
Annexes 
A Mode d’emploi de SPuTNIk

Télécharger le rapport complet