Laurence Calzone

Nathalie Chabrier-Rivier

François Fages - Loïc

Fosse -Sylvain Soliman

Langages formels dans la machine abstraite biochimique BIOCHAM

Keywords: Biologie des systèmes, logique temporelle, apprentissage, cycle cellulaire Systems biology, temporal logics, machine learning, cell cycle 1 re soumission

Le développement de langages formels pour modéliser les systèmes biologiques ouvre la voie à la conception de nouveaux outils de raisonnement automatique destinés au biologiste modélisateur. La machine abstraite biochimique BIOCHAM est un environnement logiciel qui offre un langage simple de règles pour modéliser des interactions biomoléculaires, et un langage puissant fondé sur la logique temporelle pour formaliser les propriétés biologiques du système. En s'appuyant sur ces deux langages formels, il devient possible d'utiliser des techniques d'apprentissage automatique pour inférer de nouvelles règles de réaction, estimer les valeurs des paramètres cinétiques, et corriger ou compléter les modèles semi-automatiquement. Dans cet article, nous décrivons les langages implantés dans BIOCHAM et illustrons l'utilisation du système d'apprentissage automatique sur un modèle simple du contrôle du cycle cellulaire.

ABSTRACT. With the advent of formal languages for modeling biological systems, the design of automated reasoning tools to assist the biologist becomes possible. The biochemical abstract machine BIOCHAM software environment offers a rule-based language to model bio-molecular interactions and a powerful temporal logic based language to formalize the biological properties of the system. Building on these two formal languages, machine learning techniques can be used to infer new molecular interaction rules from temporal properties, or to estimate kinetic parameter values, in order to semi-automatically correct or complete models from observed biological properties of the system. In this article we describe the formal languages of BIOCHAM and illustrate, on a simple cell cycle control model, the use of the machine learning system.

Introduction

La production en masse de données post-génomiques, telles que l'expression des ARN, la production de protéines et les interactions protéine-protéine, donne plus d'acuité au besoin, et aux difficultés, qu'il y a à raisonner à l'échelle d'un système complexe d'interactions biomoléculaires. Les connaissances concernant les réseaux d'influence géniques et les réseaux d'interaction protéiques sont intégrées dans des bases de données comme KEGG [START_REF] Kanehisa | KEGG : Kyoto Encyclopedia of Genes and Genomes[END_REF], BioCyc [START_REF] Keseler | EcoCyc : a comprehensive database resource for Escherichia coli[END_REF], . . . et présentées généralement sous forme de diagrammes annotés. Les outils comme BioSpice [START_REF] Giaever | Functional profiling of the Saccharomyces cerevisiae genome[END_REF], Copasi (VBI n.d.), GON, E-cell [START_REF] Tomita | The E-CELL Project : Towards Integrative Simulation of Cellular Processes[END_REF] etc. ont été développés pour faire des simulations numériques lorsque les données cinétiques sont connues.

Cependant, au delà de la simulation, la possibilité de faire des calculs symboliques sur les réseaux d'interactions moléculaires ouvre la voie à la réalisation de nouveaux outils de raisonnement automatique destinés au biologiste modélisateur. Notre projet avec la Machine Abstraite Biochimique BIOCHAM 1 [START_REF] Fages | Modelling and Querying Interaction Networks in the Biochemical Abstract Machine BIOCHAM[END_REF], qui a débuté en 2002, est une tentative dans cette direction. BIOCHAM est fondé sur un langage de composés moléculaires et de règles de réaction, qui fournit des sémantiques précises aux diagrammes d'interactions biomoléculaires, à trois niveaux d'abstraction :

1) la sémantique Booléenne, qui permet de raisonner sur la présence ou l'absence des molécules, associe à un ensemble de règles de réaction un système de transitions concurrent asynchrone ;

2) la sémantique des concentrations, qui permet de raisonner sur la dynamique continue du système, associe des valeurs réelles aux molécules représentant leur concentration, et associe aux règles un système d'équations différentielles (non linéaires) du premier ordre ;

3) la sémantique des populations, qui permet de raisonner sur les petits nombres de molécules, associe aux molécules une valeur entière représentant leur quantité, et associe aux règles une chaîne de Markov en temps continu.

Dans les trois cas, nous utilisons la logique temporelle comme langage de formalisation des propriétés biologiques du système. Dans la sémantique Booléenne, la logique des arbres de calcul CTL [START_REF] Clarke | Model Checking[END_REF] permet d'exprimer un large spectre de propriétés biologiques telles que l'accessibilité, les points de passage obligé, la stabilité ou les oscillations [START_REF] Chabrier | Symbolic model cheking of biochemical networks[END_REF]. Nous avons montré le passage à l'échelle des algorithmes de model-checking dans ce contexte, sur un modèle Booléen du contrôle du cycle cellulaire des mammifères, développé d'après le diagramme de Kohn [START_REF] Kohn | Molecular Interaction Map of the Mammalian Cell Cycle Control and DNA Repair Systems[END_REF] et impliquant environ 500 variables et 800 règles de réaction (Chabrier-Rivier et al., 2004a). Dans la sémantique des concentrations, la logique LTL [START_REF] Clarke | Model Checking[END_REF] avec contraintes numériques permet d'exprimer des propriétés temporelles quantitatives. Dans la sémantique des populations, la logique PCTL [START_REF] Hansson | A Logic for Reasoning about Time and Reliability[END_REF][START_REF] Kwiatkowska | PRISM 2.0 : A Tool for Probabilistic Model Checking[END_REF] avec contraintes numériques permet d'exprimer les probabilités de réalisation de propriétés temporelles quantitatives.

La vérification automatique des propriétés d'un système biologique formalisées en logique temporelle, permet non seulement d'interroger un modèle, mais aussi d'aider à le corriger, en prenant ces propriétés comme une spécification et en cherchant automatiquement des règles et des paramètres permettant de la satisfaire [START_REF] Calzone | Machine learning biomolecular interactions from temporal logic properties[END_REF]. Dans cet article, nous décrivons les deux méthodes d'apprentissage implantées dans BIOCHAM. La première, qui utilise la sémantique Booléenne, concerne la recherche de règles de réaction à ajouter ou retirer au modèle de façon à satisfaire une spécification CTL. La seconde, qui utilise la sémantique des concentrations, concerne la recherche de paramètres cinétiques permettant de satisfaire une spécification LTL avec contraintes. Les deux méthodes sont illustrées sur un modèle simplifié du contrôle du cycle cellulaire.

Travaux reliés

Les travaux pionniers de Regev, Silverman and Shapiro (2001), sur l'utilisation du π-calcul de Milner pour modéliser les voies de signalisation dans la cellule, ont été la source d'inspiration de nombreux travaux ultérieurs dans la lignée des calculs de processus [START_REF] Cardelli | Brane Calculi -Interactions of Biological Membranes[END_REF][START_REF] Regev | BioAmbients : An Abstraction for Biological Compartments[END_REF], Danos et al., 2004) et de leurs extensions stochastiques (Phillips et al., to appear). Le langage de règles de BIOCHAM se veut une simplification de cette démarche utilisant les calculs de processus, au profit d'un formalisme simple de règles de réaction, qui d'une part est très naturel pour le biologiste, et d'autre part se prête mieux aux techniques d'analyse par model-checking. Dans le cas Booléen, ce langage est similaire aux règles de réécriture de Pathway Logic implantées dans le système Maude [START_REF] Eker | Pathway Logic : Symbolic Analysis of Biological Signaling[END_REF]. L'utilisation de la logique temporelle et des techniques de model-checking en biologie des systèmes est une idée récente, elle apparaît dans [START_REF] Eker | Pathway Logic : Symbolic Analysis of Biological Signaling[END_REF][START_REF] Chabrier | Symbolic model cheking of biochemical networks[END_REF] pour les modèles Booléens, dans [START_REF] Bernot | A Fruitful Application of Formal Methods to Biological Regulatory Networks : Extending Thomas' Asynchronous Logical Approach with Temporal Logic[END_REF][START_REF] Batt | Model Checking Genetic Regulatory Networks using GNA and CADP[END_REF][START_REF] Calder | Analysis of Signalling Pathways using the Prism Model Checker[END_REF] pour les modèles discrets et dans [START_REF] Antoniotti | Model Building and Model Checking for Biochemical Processes[END_REF] pour les modèles continus.

Les techniques d'apprentissage automatique comme la programmation logique inductive [START_REF] Muggleton | Inverse Entailment and Progol[END_REF] ou la programmation génétique, ont été utilisées pour inférer des fonctions de gènes [START_REF] Bryant | Combining Inductive Logic Programming, Active Learning and Robotics to Discover the Function of Genes[END_REF], des descriptions de voies métaboliques (Angelopoulos et al., 2002a[START_REF] Angelopoulos | Slps for Probabilistic Pathways : Modeling and Parameter estimation[END_REF][START_REF] Koza | Reverse Engineering of Metabolic Pathways from Observed Data Using Genetic Programming[END_REF] ou des interactions entre gènes [START_REF] Bernot | A Fruitful Application of Formal Methods to Biological Regulatory Networks : Extending Thomas' Asynchronous Logical Approach with Temporal Logic[END_REF]. Nos travaux se rapprochent du domaine de la découverte scientifique qualitative et numérique [START_REF] Langley | Scientific Discovery : Computational Explorations of the Creative Processes[END_REF] et de la révision de théories [START_REF] Todorovski | Theory Revision in Equation Discovery[END_REF][START_REF] De Raedt | Interactive Theory Revision, an inductive Logic Programming Approach, Knowledge-Based Systems[END_REF]. Cependant, l'apprentissage à partir de formules temporelles est un sujet nouveau, que ce soit du point de vue de l'apprentissage automatique, ou de la biologie des systèmes.

Langage de description de processus biomoléculaires

Les règles de réaction BIOCHAM représentent des réactions entre des objets formels représentant des composés chimiques ou biochimiques, qui s'étendent des petites molécules aux gènes et aux protéines. BIOCHAM possède aussi un langage de schémas avec contraintes que nous ne décrivons pas ici, et qui est utilisé pour spécifier de façon concise des ensembles de molécules et de règles de réaction, ou pour indiquer la forme des règles à apprendre.

Exemple 1 Le cycle cellulaire décrit le mécanisme bio-moléculaire par lequel une cellule est amenée à se diviser et proliférer. Il s'agit d'une séquence d'événements qui inclue la réplication de l'ADN (phase S), la mitose (phase M) et l'apparition et la disparition successives de diverses protéines comme les cyclines. Pour être actives, ces cyclines s'associent à d'autres protéines, des kinases cycline-dépendantes. Les complexes ainsi formés contrôlent le déroulement du cycle.

L'exemple BIOCHAM ci-dessous est un modèle simple du cycle cellulaire tiré de [START_REF] Tyson | Modeling the cell division cycle : cdc2 and cyclin interactions[END_REF] prmeter@kIDHFHISAF prmeter@kPDHFHISAF prmeter@kQDPHHAF prmeter@kRpDHFHIVAF prmeter@kRDIVHAF prmeter@kSDHAF prmeter@kTDIAF prmeter@kUDHFTAF prmeter@kVDIHHAF prmeter@kWDIHHAF present@gdPDIAF mkesentnotpresentF

Les dernières instructions définissent l'état initial du modèle, dans lequel gdP est présent (avec une concentration initiale égale à 1) et les autres composés sont absents (concentration 0).

La sémantique des règles BIOCHAM est définie à trois niveaux d'abstraction qui correspondent à trois formes de raisonnement sur les systèmes biologiques et à trois sémantiques formelles du langage des règles : la sémantique Booléenne, la sémantique des concentrations en nombres réels, et la sémantique des populations de molécules en nombres entiers.

Sémantique Booléenne

Dans la sémantique Booléenne, on associe à chaque objet une variable Booléenne qui représente sa présence ou absence dans le système, et les règles de réaction sont interprétées comme un système de transitions concurrent asynchrone. Une règle comme A + B => C + D définit quatre transitions possibles qui correspondent à la consommation totale ou non des réactifs A et B par la réaction. Pour que la règle soit applicable, les molécules A et B doivent être présentes dans l'état courant du système. Dans l'état suivant les molécules C et D sont présentes tandis que les molécules A et B peuvent être absentes (consommation totale) ou présentes (consommation partielle). Ceci se fait de façon non-déterministe afin de rendre compte de tous les comportements possibles. Le système de transitions est asynchrone dans le sens où une seule transition est appliquée à la fois, ce qui permet de représenter les phénomènes biologiques fondamentaux de masquage d'une réaction par une autre, comme par exemple dans le cas de certaines inhibitions. A l'inverse, l'hypothèse de synchronie (dans la- quelle toutes les règles applicables se déclenchent en même temps) ne permettrait pas de représenter correctement les phénomènes de compétition entre réactions.

Formellement, la sémantique Booléenne d'un modèle BIOCHAM est définie par une structure de Kripke K = (S, R), où S est l'ensemble des états définis par le vecteur de variables Booléennes, R ⊆ S × S est la relation de transition entre états supposée totale (c'est-à-dire que pour tout état s ∈ S, il existe un état s ∈ S tel que (s, s ) ∈ R). Lorsque la relation de transition provenant uniquement des réactions BIOCHAM n'est pas totale, elle est automatiquement complétée par une transition qui ajoute une boucle sur les états n'ayant pas de successeur. Un chemin dans K partant de l'état s 0 est une suite infinie d'états π = s 0 , s 1 , ... telle que (s i , s i+1 ) ∈ R pour tout i ≥ 0. π k désignera le chemin s k , s k+1 , • • • . L'évolution temporelle du système est modélisée par la succession des étapes de transition, et les différents comportements du système sont modélisés par le choix non-déterministe des transitions. La figure 1 montre une trace de simulation Booléenne du modèle du cycle cellulaire.

Sémantique des concentrations

Dans la sémantique des concentrations, on associe à chaque objet une variable réelle qui représente sa concentration, et les règles sont interprétées par un système d'équations différentielles. Formellement, à un ensemble de règles BIOCHAM E =

dx k /dt = n i=1 r i (x k ) * e i - n j=1 l j (x k ) * e j
où r i (x) (resp. l i ) est le coefficient stoechiométrique de x dans le membre droit (resp. gauche) de la règle i.

Étant donné un état initial défini par les concentrations initiales des objets, les méthodes d'intégration numérique permettent de simuler le comportement du système. La figure 2 illustre par exemple les oscillations des concentrations des protéines décrites dans l'exemple 1.

Sémantique des populations

Dans la sémantique des populations, on associe à chaque objet une variable entière qui représente le nombre de molécules, et les règles sont interprétées par une chaîne de Markov en temps continu dans laquelle les taux de transition sont définis par les expressions cinétiques des règles BIOCHAM.

Les méthodes de simulation stochastique [START_REF] Gillespie | General Method for Numerically Simulating Stochastic Time Evolution of Coupled Chemical-Reactions[END_REF] calculent des réalisations du processus. Ces courbes sont en général des versions bruitées des courbes obtenues dans la sémantique des concentrations. Cependant dans les modèles faisant intervenir, par exemple, des petits nombres de molécules, des comportements qualitativement différents peuvent apparaître dans les simulations stochastiques. Un exemple classique de ce phénomène est le modèle du phage lambda dans lequel un petit nombre de molécules, facteurs de promotion de deux gènes, engendrent un phénomène de multiplication explosive (lyse) après une période d'attente passive plus ou moins longue (lysogénie) [START_REF] Gibson | A probabilistic model of a prokaryotic gene and its regulation[END_REF].

Dans la sémantique stochastique, pour un volume V donné, les concentrations sont donc converties en nombre de molécules (N = C ×V ×K, K constante d'Avogadro). Les coefficients cinétiques k sont convertis en taux de transition τ par les formules suivantes :

-τ = k pour les réactions du type A =>..., -τ = k V ×K pour les réactions du type A+B=>..., -τ = 2 × k V ×K pour les réactions du type A+A=>..., -τ = k (V ×K) 2 pour les réactions du type A+B+C=>..., -etc.

Langage de description de propriétés biologiques

Le trait le plus original de BIOCHAM est l'utilisation de la logique temporelle comme langage de formalisation des propriétés biologiques d'un modèle. Les logiques CTL Booléenne, LTL et PCTL avec contraintes numériques sont utilisées respectivement dans chacun des trois niveaux sémantiques. Elles permettent de formaliser les propriétés importantes d'un modèle, qui aujourd'hui sont expliquées de façon informelle dans les articles de modélisation. Leur formalisation permet cependant de les considérer comme une spécification devant être préservée lors des opérations de raffinement, de simplification ou de composition de modèles. La possibilité de vérifier automatiquement une telle spécification permet d'envisager la comparaison et la réutilisation de différents modèles de façon beaucoup plus systématique que ce qui se pratique aujourd'hui dans ce domaine.

Logique CTL pour la sémantique Booléenne

La logique temporelle des arbres de calcul (Computation Tree Logic CTL * ) [START_REF] Clarke | Model Checking[END_REF] est une extension de la logique classique qui permet de raisonner sur un arbre de transitions d'états à l'aide d'opérateurs sur les choix de branches (nondéterminisme) et sur le temps (transitions d'états). Deux quantificateurs de chemins sont introduits pour raisonner sur le non-déterminisme : Aφ qui signifie que φ est vraie sur tous les chemins, et Eφ, pour φ est vraie sur au moins un chemin. Plusieurs opérateurs temporels sont introduits : Xφ signifiant que φ est vrai à la transition suivante, Gφ pour φ est toujours vrai, F φ pour φ finit par devenir vrai, φ U ψ pour ψ finit par devenir vrai et φ est vrai jusqu'à ce que ψ soit vrai, φ W ψ pour φ est vrai jusqu'à ce qu'éventuellement ψ devienne vrai. Dans cette logique, F φ est équivalent à true U φ, Gφ est équivalent à φ W f alse et nous avons les propriétés de dualité suivantes :

¬EF (φ) = AG(¬φ), ¬Eφ U ψ = A¬ψ W ¬φ et ¬Eφ W ψ = A¬ψ U ¬φ. s |= α ssi α est une formule propositionnelle et s |= α, s |= Eψ ssi il existe un chemin π partant de s tel que π |= ψ, s |= Aψ ssi pour tout chemin π partant de s, π |= ψ, π |= φ ssi s |= φ où s est l'état de début de π, π |= Xψ ssi π 1 |= ψ, π |= F ψ ssi il existe k ≥ 0 tel que π k |= ψ, π |= Gψ ssi pour tout k ≥ 0, π k |= ψ, π |= ψU ψ ssi il existe k ≥ 0 tel que π k |= ψ et π j |= ψ pour tout 0 ≤ j < k.
Tableau 1. Définition inductive de la relation de vérité d'une formule CTL dans un état ou un chemin, dans une structure de Kripke K.

Dans le fragment CTL qui nous intéresse, chaque opérateur temporel doit obligatoirement être précédé d'un quantificateur de chemin et un quantificateur de chemin doit obligatoirement être suivi d'un opérateur temporel. Le tableau 1 définit la relation de vérité d'une formule dans une structure de Kripke K = (S, R) dont les états sont définis par des variables Booléennes. Un modèle BIOCHAM M est constitué d'un ensemble de règles et d'une spécification (partielle) de l'ensemble des états initiaux. Les règles de M définissent la structure de Kripke associé au modèle, et on note M |= φ si la formule φ est vrai dans tous les états initiaux spécifiés par M.

Expressivité

Comme montré dans [START_REF] Chabrier | Symbolic model cheking of biochemical networks[END_REF], le fragment CTL a une expressivité suffisante pour exprimer un très large éventail de propriétés biologiques : Sur l'accessibilité. Existe-t-il un chemin pour produire (c-à-d synthétiser, activer etc.) une protéine P ? Cette requête est formalisée par la formule CTL EF (P ), abrégée en rehle@A dans BIOCHAM.

Sur les voies.

Est-ce que l'état Q est un point de passage obligé pour atteindre les états où P ∼ {s} est présent ? La formule correspondante ¬(E((¬(Q)U P ∼ {s}))) est abrégée en hekpoint@D£{s}A.

Sur la stationnarité et la stabilité. Est-ce qu'un état (ou un ensemble d'états) de la cellule décrit par une formule s est un état stationnaire (resp. stable) ? s ⇒ EG(s) (resp. s ⇒ AG(s)), abrégé en stedy@sA (resp. stle@sA).

Sur les oscillations.

Le système peut-il exhiber un comportement cyclique sur la présence de P ? EG((P ⇒ EF ¬P )∧(¬P ⇒ EF P )), abrégé en osil@A. rehle@gdP~{pI}A rehle@gylinA rehle@gylin~{pI}A rehle@gdPEgylin~{pI}A rehle@gdP~{pI}Egylin~{pI}A osil@gdPA osil@gdP~{pI}A osil@gdP~{pI}Egylin~{pI}A osil@gylinA eq@@3@gdPEgylin~{pI}AAEb hekpoint@gdP~{pI}Egylin~{pI}DgdPEgylin~{pI}AA

Algorithme de model-checking symbolique

Le langage de requête CTL pour les modèles Booléens est implanté en BIOCHAM avec une interface au vérificateur de modèles symbolique NuSMV [START_REF] Cimatti | NuSMV 2 : An OpenSource Tool for Symbolic Model Checking[END_REF]. Les performances obtenues sur un gros modèle du contrôle du cycle cellulaire suivant la carte de Kohn [START_REF] Kohn | Molecular Interaction Map of the Mammalian Cell Cycle Control and DNA Repair Systems[END_REF] comprenant 800 règles et 500 variables sont de l'ordre de quelques dizaines de secondes pour compiler le modèle ou vérifier des formules CTL simples (Chabrier-Rivier et al., 2004a). Cependant ces performances sont en deçà de celles obtenues classiquement par NuSMV avec des modèles fortement structurés pour la vérification de circuits ou de programmes. Les modèles Booléens de BIOCHAM ont la caractéristique d'être très fortement non-déterministes, ce qui peut conduire à des problèmes de performance sur des problèmes de plus petite taille comprenant un grand nombre de voies parallèles, comme par exemple dans le modèle de la voie de signalisation MAPK de [START_REF] Schoeberl | Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors[END_REF].

Équivalence CTL de modèles

Il existe, en général, de nombreuses manières de modéliser un même système, suivant le niveau de détail souhaité et l'information disponible. Même pour une simple réaction enzymatique, il existe déjà deux façons classiques de la coder en BIOCHAM :

-soit en détaillant la formation d'un complexe intermédiaire

A + B <=> D D => A + C -soit directement A + B => A + C (que l'on peut aussi noter B =[A]=> C)
Le premier choix peut sembler plus simple et naturel pour une conversion en BIO-CHAM, alors que le second correspond à la traduction directe de certains modèles numériques à base d'équations différentielles. Les deux propriétés ci-dessous montrent l'équivalence des deux modèles vis-à-vis de certaines CTL ;

Nous supposerons que le complexe D n'apparaît que dans les deux règles du premier cas, sans quoi la simplification n'est bien sûr pas correcte. Les structures de Kripke associées à ces deux choix définissent deux relations de vérité :

|= 1 et |= 2 .
Proposition 1 (Accessibilité) (Chabrier-Rivier et

I |= 1 ¬E(¬φU ψ) implique I |= 2 ¬E(¬φU ψ).

LTL avec contraintes numériques pour la sémantique des concentrations

Une extension de la logique temporelle du temps linéaire (Linear Time Logic, LTL) avec contraintes arithmétiques est utilisée pour la sémantique (déterministe) des concentrations moléculaires. Cette approche est similaire à celle décrite dans [START_REF] Antoniotti | Model Building and Model Checking for Biochemical Processes[END_REF] et utilisée dans le projet DARPA BioSpice.

LTL est un fragment de la logique temporelle, qui utilise uniquement les opérateurs temporels (pas les quantificateurs de chemin), et est donc appropriée pour raisonner sur les systèmes déterministes que sont les modèles cinétiques. Dans LTL avec contraintes arithmétiques, les formules atomiques sont des formules de logique du premier ordre avec égalité, inégalité et opérations arithmétiques. Ces formules portent sur les valeurs réelles des concentrations, par exemple > £{s}, ainsi que celles de leurs dérivées, par exemple d@AGdt < H.

Expressivité

Les propriétés d'accessibilité sont formalisées par l'opérateur F . Par exemple, la propriété que la concentration de la protéine P phosphorylée sur s peut être dépassée par celle de la protéine P , s'exprime par la formule F ([P ] > [P ∼ {s}]). La propriété que la concentration de P atteint et reste au dessus d'un certain seuil c s'exprime par F G([P ] > c).

Les propriétés d'oscillation d'une concentration [P ] un nombre n de fois, portent sur le nombre d'alternances du signe de la dérivée. Elles s'expriment par la formule

F ((d[P ]/dt > 0) ∧ F ((d[P ]/dt < 0) ∧ F ((d[P ]/dt > 0) . . .)))
abrégée par osil@DnA.

Algorithme de model-checking avec contraintes

Sous l'hypothèse que les conditions initiales sont complètement déterminées, les méthodes d'intégration numérique du système d'équations différentielles associées aux règles BIOCHAM, fournissent une trace de simulation discrète (à pas variable) qui constitue une structure de Kripke linéaire, dans laquelle on peut interpréter les formules LTL. De façon à poser des contraintes sur non seulement les valeurs des concentrations mais aussi sur leurs dérivées, on considère des traces de simulation de la forme

(< t 0 , x 0 , dx 0 /dt >, < t 1 , x 1 , dx 1 /dt >, ...)
qui contiennent pour chaque point de temps t i , la valeur des concentrations x i et de leurs dérivées dx i /dt.

Sous l'hypothèse que les propriétés temporelles peuvent être vérifiées sur un horizon fini, l'algorithme de model-checking avec contraintes pour vérifier une formule φ consiste alors à effectuer les opérations suivantes : 1) calculer une trace de simulation finie, 2) étiqueter les points de la trace par les sous-formules atomiques de φ qui sont vraies en ce point ;

3) ajouter les sous-formules de la forme F φ 1 aux points prédécesseurs des points étiquetés par φ 1 , 4) ajouter les sous-formules de la forme φ 1 U φ 2 aux points prédécesseurs des points étiquetés par φ 2 tant qu'ils satisfont φ 1 On a Gφ 1 = ¬F (¬φ 1 ), ce qui se ramène ici à dire qu'on peut étiqueter le dernier état par Gφ 1 s'il est étiqueté par φ 1 , puis étiqueter les prédécesseurs d'un état étiqueté par Gφ 1 , s'il est lui même étiqueté par φ 1 .

Pour des raisons d'intégration au système BIOCHAM, cet algorithme est implanté en Prolog, de même que les méthodes d'intégration numérique (méthode de Runge-Kutta à pas variable, et méthode implicite de Rosenbrock pour les systèmes raides).

PCTL avec contraintes arithmétiques pour la sémantique des populations

La logique PCTL [START_REF] Hansson | A Logic for Reasoning about Time and Reliability[END_REF] remplace essentiellement les opérateurs de choix des chemins E et A de la logique CTL, par un opérateur P p , désignant la probabilité de réalisation de la formule sous sa portée. Par exemple, A(ψU ψ ) devient P ≥1 (ψU ψ ) (imposant que la probabilité que ψU ψ se réalise soit égale à 1). Comme dans la section précédente, on s'intéresse à des formules atomiques du premier ordre avec contraintes arithmétiques portant, cette fois-ci, sur des nombres entiers de molécules.

Les algorithmes de model-checking probabiliste existants, comme celui du système PRISM [START_REF] Kwiatkowska | PRISM 2.0 : A Tool for Probabilistic Model Checking[END_REF], ne permettent cependant pas de traiter les exemples fortement non déterministes ou pour lesquels les variables ont un grand domaine, ce qui est le cas pour la sémantique de population de BIOCHAM. De meilleurs résultats ont été obtenus avec un algorithme de type Monte-Carlo semblable à celui utilisé dans le système APMC [START_REF] Hérault | Approximate Probabilistic Model Checking[END_REF]. Cet algorithme permet de vérifier une formule LTL avec contraintes, et fournit une probabilité approchée de réalisation de cette formule. Pour ce faire, un ensemble de simulations stochastiques est effectué, en étiquetant les traces de chaque simulation par l'algorithme pour LTL avec contraintes. La probabilité de réalisation est alors estimée par comptage sur l'échantillon. Dans ce cadre, l'opérateur P p n'est autorisé qu'en tête de la formule. Ce ne sont donc pas des formules PCTL générales qui sont vérifiées, mais des formules LTL en fournissant une probabilité de réalisation. Le résultat est une probabilité, et non pas la vérification d'une formule comprenant un ensemble de probabilités imbriquées.

Apprentissage de modèles

Les modèles biologiques sont construits à partir d'expériences réalisées sur des organismes sauvages ou mutés. Ces expériences définissent les propriétés biologiques que les modèles doivent satisfaire. En particulier, les expériences sur les organismes mutés, dans lesquels la synthèse de certaines protéines est bloquée, permettent de discriminer entre différentes hypothèses d'interaction. Dans BIOCHAM, l'utilisation de langages formels pour, à la fois décrire le système (langage de règles), et décrire ses propriétés (spécification en logique temporelle), permet d'appliquer des algorithmes d'apprentissage automatique pour aider à corriger ou compléter un modèle en fonction de nouvelles spécifications.

Recherche de règles

En ce qui concerne la recherche de règles à ajouter ou supprimer au modèle, nous utilisons la sémantique Booléenne de BIOCHAM, et une spécification CTL des propriétés attendues du modèle.

Les formules CTL sont partitionnées en trois classes de façon à guider la recherche pour l'ajout ou la suppression de règles. La classe ECTL est formée des formules CTL qui ne contiennent pas l'opérateur A (c'est-à-dire qui ne contiennent pas d'occurrence positive de A ni d'occurrence négative de E). La classe ACTL est formée des formules CTL qui ne contiennent pas l'opérateur E. On note UCTL l'ensemble des autres formules.

Proposition 3 Soit K = (S, R, L) et K = (S, R , L) deux structures de Kripke telles que R ⊆ R . Pour toute formule ECTL φ, si s |= K φ alors s |= K φ. Pour toute formule ACTL φ, si s |= K φ alors s |= K φ. Preuve 1 Soit φ une formule ECTL telle que s |= K φ. Raisonnons par l'absurde et supposons que s |= K φ, alors il existe un chemin π = ss 1 s 2 ... dans la structure K tel que π |= K φ ou s |= K φ . Comme R ⊆ R , π est aussi un chemin dans la structure K , on a donc π |= K φ ou s |= K φ, contradiction. De même pour une formule ACTL φ, par dualité ¬φ est une formule ECTL. Si s |= K φ alors s |= K ¬φ, d'où l'on déduit comme précédemment s |= K ¬φ, c'est-à-dire s |= K φ.
Cette propositon montre que l'on peut chercher à satisfaire une spécification CTL à partir d'un modèle ne les satisfaisant que partiellement, en supprimant des règles du modèle qui rendent les formules ACTL fausses, et en ajoutant au modèle des règles qui permettent de satisfaire les formules ECTL fausses.

La commande learn_one_rule(reaction_pattern, spec_CTL) énumère chaque instance du schéma de règles fourni en premier argument qui, ajoutée au modèle, permet de satisfaire la formule CTL fournie en second argument (ou implicitement avec la commande add_spec). La commande delete_one_rule(reaction_pattern, spec_CTL) énumère chaque instance de règle qu'il suffit de retirer au modèle pour satisfaire la formule CTL.

Exemple 3 Dans cet exemple, on retire la règle suivante du modèle de Tyson : deleterule@gylinCgdP~{pI}abgdP~{pI}Egylin~{pI}AF Le complexe inactif gdP~{pI}Egylin~{pI} ne pouvant plus se former, certaines molécules ne sont plus accessibles. La spécification donnée dans l'exemple 2 n'est donc plus vérifiée. L'algorithme de recherche de règle propose quatre règles parmi 587 testées. L'ajout de l'une d'elles permet de satisfaire à nouveau la spécification : biocham: learn_one_rule(elementary_interaction_rules).

(1) Cyclin+Cdc2~{p1}=

[Cdc2]=>Cdc2~{p1}-Cyclin~{p1} (2) Cyclin+Cdc2~{p1}=[Cyclin]=>Cdc2~{p1}-Cyclin~{p1} (3) Cyclin+Cdc2~{p1}=>Cdc2~{p1}-Cyclin~{p1} (4) Cyclin+Cdc2~{p1}=[Cdc2~{p1}]=>Cdc2~{p1}-Cyclin~{p1}
La troisième règle correspond à la règle retirée. les autres règles ne sont pas biologiquement plausibles et auraient pu être éliminées d'avance en fournissant plus d'information dans le schéma de règles à chercher.

On définit un algorithme de révision du modèle permettant d'ajouter et retirer plusieurs règles, en fixant une stratégie de recherche des règles à ajouter ou supprimer de façon à satisfaire la spécification. La commande revise_model(reaction_pattern1, reaction_pattern2, spec_CTL) énumère de façon heuristique des ensembles d'ajouts et de retraits de règles qui permettent de satisfaire une spécification CTL fournie en argument. L'algorithme utilisé dans BIOCHAM est défini par un système de transitions non-déterministe et backtrackant. L'état d'une recherche est représenté par un triplet < QT, Q, R > où -R est l' ensemble courant des règles de réaction.

-QT = (E, U , A) est l'ensemble des formules déjà traitées et satisfaites pour R. Ces formules sont triées en 3 groupes : E représente les formules ECTL, U les formules non classifiables (UCTL) et A les formules ACTL.

-Q = (E, U, A) est l'ensemble des formules à traiter triées de la même manière en formules ECTL, UCTL et ACTL.

La spécification I de l'état initial est fixe et on s'intéresse donc à une modification de l'ensemble R des règles du modèle de façon à satisfaire la spécification CTL, ce que l'on note R |= φ par abus de notation. Les modifications sont cherchées en considérant seulement les instances des schémas de règles fournis pour les ajouts et les retraits. De plus, la recherche des règles à retirer (dans les transitions U" et A' ci-dessous) ne se fait que dans l'ensemble des règles apparaissant dans le contre-exemple fourni par le model-checker. Exemple 4 L'algorithme de révision peut être utilisé pour raffiner un modèle, de façon par exemple à prendre en compte une nouvelle molécule donnée avec une spécification de ses propriétés biologiques. Dans l'exemple du cycle cellulaire, on peut ainsi s'intéresser à la recherche d'une modification du modèle, qui préserve les spécifications précédentes, et fasse intervenir une molécule X telle que (1) X ne peut être activée sans la présence de Cdc2-Cyclin~{p1}, et (2) sans X, Cyclin~{p1} n'est jamais inactivée ou dégradée :

État initial :

< (∅, ∅, ∅), (E, U, A), R > transition E : < (E, U , A), (E ∪ {e}, U, A), R > → < (E ∪ {e}, U , A), (E, U, A), R > si R |= e transition E' : < (E, U , A), (E ∪{e}, U, A), R > → < (E ∪ {e}, U , A), (E, U, A), R ∪ {r} > si R |= e et ∀f ∈ {e} ∪ E ∪ U ∪ A K ∪ {r} |= f transition U : < (E, U , A), (∅, U ∪{u}, A), R > → < (E, U ∪ {u}, A), (∅, U, A), R > si R |= u transition U' : < (E, U , A), (∅, U ∪ {u}, A), R > → < (E, U ∪ {u}, A), (∅, U, A), R ∪ {r} > si R |= u et ∀f ∈ {u} ∪ E ∪ U ∪ A R ∪ {r} |= f transition U" : < (E, U , A), (∅, U ∪ {u}, A), R ∪ Re > → < (E, U ∪ {u}, A), (∅, U, A), R >si K : s i |= u et ∀f ∈ {u} ∪ E ∪ U ∪ A R |= f transition A : < (E, U , A), (∅, ∅, A ∪ {a}), R > → < (E, U , A ∪ {a}), (E p , U p , A), R >si R |= a transition A' : < (E ∪ E p , U ∪ U p , A), (∅, ∅, A ∪ {a}), R ∪ Re > → < (E, U , A ∪ {a}), (E p , U p , A), R >si R |= a, ∀f ∈ {u} ∪ E ∪ U ∪ A R |= f et E p ∪ U p est l'
(1) eq@@3@AA Eb hekpoint@gdPEgylin~{pI}DAAD (2) eq@@gylin~{pI}A Eb hekpoint@D3@gylin~{pI}AAA.

A cette spécification, sont également rajoutées des formules d'accessibilité et d'oscillation relatives à X :

rehle@AD rehle@3@AAD osil@A.

Par défaut, on associe à la molécule X une dynamique arbitraire avec deux règles simples de synthèse et de dégradation : _=>X et X=>_.

Ces spécifications ne sont pas vérifiées par le modèle. La commande revise_model, limitée à des règles faisant intervenir X pour l'ajout de règles, permet de trouver des solutions pour le corriger. La première solution trouvée : heletion@sAX abF gylin~{pI}abF eddition@sAX agdPEgylin~{pI}abF Cgylin~{pI}abEgylin~{pI}F propose non seulement que la synthèse de X dépende directement de Cdc2-Cyclin~{p1} mais aussi que X se complexe à Cyclin~{p1} afin de l'inactiver. Ces modifications apportées au modèle (retraits et ajouts des nouvelles règles) sont assurées de respecter les spécifications.

D'autres solutions sont également proposées, comme par exemple :

heletion@sAX abF gylin~{pI}abF eddition@sAX agdPEgylin~{pI}abF gylin~{pI}aabF où la synthèse de X est à nouveau maintenue par Cdc2-Cyclin~{p1} mais où Cyclin~{p1} est directement dégradée par X.

Recherche de paramètres

De la même manière que pour l'apprentissage de règles d'interaction à partir de propriétés CTL, il est possible d'utiliser une spécification LTL avec contraintes arithmétiques pour l'apprentissage de valeurs de paramètres cinétiques. Nous avons développé une méthode énumérative dans laquelle l'espace de recherche, fourni par des bornes sur les valeurs des paramètres, est exploré avec une précision spécifiée par le modélisateur. Pour chaque jeu de valeurs testées, le modèle est simulé et confronté à la spécification LTL.

La recherche de paramètres peut être illustrée à partir du modèle de l'exemple 1. Deux paramètres sont choisis, k 3 et k 4 , auxquels sont assignées les valeurs arbitraires 0.01 et 0.01. La simulation montre que le système atteint un état stationnaire (figure 3) or il est censé osciller. Pour y remédier, BIOCHAM recherche, pour chaque paramètre, parmi 20 valeurs différentes et dans un intervalle [0,200], des valeurs telles que gdPEgylin~{pI} oscille au moins 3 fois en 150 unités de temps : treget@kQDkRD@HDPHHAD@HDPHHADPHDosil@gdPEgylin~{pI}DQADISHAF BIOCHAM propose une première solution (k 3 = 10, k 4 = 70) qui n'est pas satisfaisante car gdPEgylin~{pI} oscille mais avec une très faible amplitude à peine décelable et ne s'active donc pas complètement (figure 4). L'ajout de la contrainte Notons que le nombre de paramètres pour lesquels une valeur est cherchée ainsi que la précision de la recherche influencent le temps de calcul. D'une certaine manière, le processus d'apprentissage automatise une partie de ce que le modélisateur fait à la main, c'est-à-dire essayer différentes valeurs de paramètres qui correspondent à des formes de courbes d'activité de protéines qu'il souhaite reproduire, et ce dans un intervalle de valeurs de paramètres jugé réaliste. L'algorithme d'apprentissage permet de tester rapidement des espaces de paramètres après formalisation de l'aspect de la courbe en spécification LTL. Cette méthode est en fait complémentaire d'autres méthodes d'estimation de paramètres telles que l'utilisation de diagrammes de bifurcation.

Conclusion et perspectives

Le projet BIOCHAM compte parmi les premières expériences d'application des méthodes formelles de vérification et d'apprentissage automatique en modélisation de systèmes biologiques. La sémantique Booléenne est utilisée d'une part pour exprimer des propriétés d'accessibilité, de points de passage obligés, de stationnarité et d'oscillation dans les réseaux d'interaction, et d'autre part pour inférer des règles de réaction Le cadre formel développé dans BIOCHAM fournit aussi une base pour étudier les questions d'abstraction et de compositionnalité de modèles en biologie des systèmes.

En particulier nous envisageons d'explorer les techniques d'interprétation abstraite pour formaliser les problèmes de simplification et de raffinement de modèles. L'étude de la compositionnalité vise à faciliter la ré-utilisation de modèles dans différents contextes et à définir des notions d'interfaces dans les modèles ouverts. Ces questions sont actuellement abordées dans l'étude de modèles couplés du cycle cellulaire et du cycle circadien avec l'objectif de contribuer à formuler de nouvelles hypothèses biologiques par ces méthodes.

  La syntaxe des molécules et des réactions peut être donnée par la grammaire (simplifiée) suivante : objet = molecule | molecule : lieu molecule = nom | molecule-molecule |molecule∼{nom,...,nom} reaction = solution => solution | cinetique for solution => solution solution = _ | objet | nombre*objet | solution + solution L'objet de base est un composé moléculaire. Il peut être localisé avec l'opérateur " :" en indiquant un lieu, qui est un simple nom représentant un compartiment tel que le noyau, le cytoplasme, la membrane, etc... L'opérateur de liaison "-" sert à représenter la complexation ou d'autres formes de liaison de molécules. L'opérateur de modification "∼" permet d'attacher un ensemble de modifications à une molécule, tel que l'ensemble des sites phosphorylés d'une protéine. Pour les règles de réaction, les abréviations suivantes peuvent être utilisées : A = [C] => B pour la règle A+C => B + C avec le catalyseur C, et A <=> B pour les deux réactions symétriques. Par exemple, M + N => M -N représente une règle de complexation, M = [N ] => M ∼ {s} une règle de phosphorylation catalysée par une kinase N , et M : L => M : L représente une règle de transport de M , de L en L . Dans tous les cas, les règles peuvent être données avec une expression cinétique, correspondant par exemple à la loi d'action de masse.

Figure 1 .

 1 Figure 1. Exemple de trace d'une simulation Booléenne du modèle du cycle cellulaire de l'exemple 1.

Figure 2 .

 2 Figure 2. Simulation numérique du modèle du cycle cellulaire de l'exemple 1

Figure 3 .

 3 Figure 3. Simulation numérique avec k 3 = 0.01 et k 4 = 0.01

Figure 4 .Figure 5 .

 45 Figure 4. Simulation numérique avec k 3 = 10 et k 4 = 70

  Sans hypothèse d'équité forte, cette formule n'est cependant qu'une approximation car elle est impliquée par l'existence d'oscillations mais n'implique pas l'alternance infinie d'états de types P et ¬P . Le modèle du cycle cellulaire précédent vérifie, par exemple, les propriétés suivantes qui expriment que toutes les molécules sont accessibles, que leur activité oscille et qu'il est nécessaire de produire la forme doublement phosphorylée gdP~{pI}Egylin~{pI} avant l'activation de gdPEgylin~{pI} :

	Exemple 2

Point de passage obligé)

  al., 2004b) Soit φ une formule CTL atomique, si I |= 2 EF (φ) alors I |= 1 EF (φ). De plus, si A et B n'apparaissent pas négativement (i.e. sous un nombre impair de négations) dans φ et si D n'apparaît pas du tout dans φ, alors I |= 1 EF (φ) implique I |= 2 EF (φ). (Chabrier-Rivier et al., 2004b) Soit ¬E(¬φU ψ) une formule de checkpoint, i.e. φ et ψ sont des formules atomiques décrivant des états, si A et B n'apparaissent pas négativement dans ψ et D n'apparaît pas du tout dans ψ, alors I |= 2 ¬E(¬φU ψ) implique I |= 1 ¬E(¬φU ψ).

	Proposition 2 (

si A et B n'apparaissent pas négativement dans φ et D n'apparaît pas du tout dans φ, alors

  ensemble des formules rendues non-satisfaites par le retrait des règles Re. La terminaison se montre en considérant l'ordre lexicographique sur le couple < a, n > où a est le nombre de propriétés ACTL non vérifiées, et n est le nombre de propriétés ECTL ou UCTL non vérifiées. Chaque transition diminue strictement a ou bien laisse a égal et diminue strictement n. La correction de l'algorithme découle du fait que chaque transition préserve l'ensemble des propriétés en les faisant passer entre les ensembles vérifiés et non-vérifiés et ne conserve dans l'ensemble vérifié que des propriétés vraies. Cet algorithme est toutefois incomplet pour deux raisons. D'une part, la satisfaction d'une formule ECTL ou UCTL n'est recherchée qu'en ajoutant une seule règle au modèle (transitions E' et U'). D'autre part, la complétion de la relation d'accessibilité en une relation totale dans la structure de Kripke associée au modèle, ne garantit pas que les conditions de la proposition 3 sont vérifiées après l'ajout ou le retrait d'une règle. Il s'agit donc d'un algorithme de révision heuristique qui peut échouer dans des cas où la spécification CTL est satisfaisable.

	Proposition 4 L'algorithme de révision termine et si sa configuration terminale est
	de la forme < (E, U , A), (∅, ∅, ∅), R > alors le modèle R vérifie la spécification
	initiale.
	Preuve 2


Francois Fages
email: francois.fages@inria.frsylvain.soliman@inria.fr

Sylvain Soliman

François Fages

Model Revision from Temporal Logic Properties in Computational Systems Biology

Systems biologists build models of bio-molecular processes from knowledge acquired both at the gene and protein levels, and at the phenotype level through experiments done in wild-life and mutated organisms. In this chapter, we present qualitative and quantitative logic learning tools, and illustrate how they can be useful to the modeler. We focus on biochemical reaction models written in the Systems Biology Markup Language SBML, and interpreted in the Biochemical Abstract Machine BIOCHAM. We first present a model revision algorithm for inferring reaction rules from biological properties expressed in temporal logic. Then we discuss the representations of kinetic models with ordinary differential equations (ODEs) and with stochastic logic programs (SLPs), and describe a parameter search algorithm for finding parameter values satisfying quantitative temporal properties. These methods are illustrated by a simple model of the cell cycle control, and by an application to the modelling of the conditions of synchronization in period of the cell cycle by the circadian cycle.

Introduction

One promise of computational systems biology is to model biochemical processes at a sufficiently large scale so that complex system behaviors can be predicted under various conditions. The biochemical reaction systems involved in these processes may contain many cycles and exhibit complex multistationarity and oscillating behaviors. While usually neglected in metabolic networks, these characteristics are preponderant in models of signal transduction and cell control. They thus provide a challenge to representation and inference methods, and the issue of representing complex biochemical systems and their behavior at different levels of abstraction is a central one in systems biology.

The pioneering use in [START_REF] Regev | Representation and simulation of biochemical processes using the pi-calculus process algebra[END_REF] of the π-calculus process algebra for modeling cell signalling pathways, has been the source of inspiration of numerous works in the line of process calculi [START_REF] Cardelli | Brane calculi -interactions of biological membranes[END_REF][START_REF] Regev | Bioambients: An abstraction for biological compartments[END_REF][START_REF] Danos | Formal molecular biology[END_REF] and their stochastic extensions [START_REF] Phillips | A correct abstract machine for the stochastic pi-calculus[END_REF]. Recently, the question of formalizing the biological properties of the system has also been raised, and formal languages have been proposed for this task, most notably using temporal logics in either boolean [START_REF] Eker | Pathway logic: Symbolic analysis of biological signaling[END_REF][START_REF] Chabrier | Symbolic model cheking of biochemical networks[END_REF], discrete [START_REF] Bernot | A fruitful application of formal methods to biological regulatory networks: Extending thomas' asynchronous logical approach with temporal logic[END_REF][START_REF] Batt | Model checking genetic regulatory networks using gna and cadp[END_REF][START_REF] Calder | Analysis of signalling pathways using the prism model checker[END_REF] or continuous models [START_REF] Antoniotti | Model building and model checking for biochemical processes[END_REF][START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF].

The biochemical abstract machine BIOCHAM1 [START_REF] Fages | Modelling and querying interaction networks in the biochemical abstract machine BIOCHAM[END_REF][START_REF] Calzone | BIOCHAM: An environment for modeling biological systems and formalizing experimental knowledge[END_REF] has been designed as a simplification of the process calculi approach using a logic programming setting and a language of reaction rules compatible with the Systems Biology Markup Language SBML [START_REF] Hucka | The systems biology markup language (SBML): A medium for representation and exchange of biochemical network models[END_REF] (http://www.sbml.org/). This opens up the whole domain of mathematical biology, through repositories like BioModels.net (http://www.biomodels.net), CMBSlib (http://contraintes.inria.fr/CMBSlib/), PWS (http://jjj.biochem.sun.ac.za/), etc. This rule-based language is used in BIOCHAM for modeling biochemical networks at three abstraction levels :

the boolean semantics, where one reasons on the presence/absence of molecules, the differential semantics, where one reasons on molecular concentrations, the stochastic semantics, where one reasons on molecule numbers and reaction probabilities.

A second language is used to formalize the biological properties known from experiments in temporal logic (the Computation Tree Logic CTL, Linear Time Logic LTL or Probabilistic LTL with constraints, according to the qualitative, quantitative or stochastic nature of the properties). Such a formalization is a first step toward the use of logic learning tools to help the modeler in his tasks [START_REF] Fages | From syntax to semantics in systems biology -towards automated reasoning tools[END_REF]. When a model does not satisfy all the expected properties, the purpose of the machine learning system of BIOCHAM is to propose rules or kinetic parameter values in order to curate the model w.r.t. a given specification [START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF]. This novel approach to biological modeling has been applied to a data set of models about the cell cycle control in different organisms, and signal transduction network (see http://contraintes.inria.fr/APrIL2/). There has been work on the use of machine learning techniques, such as inductive logic programming (ILP, see Chapter 1 or [START_REF] Muggleton | Inverse entailment and progol[END_REF]), to infer gene functions [START_REF] Bryant | Combining inductive logic programming, active learning and robotics to discover the function of genes[END_REF], metabolic pathway descriptions [START_REF] Angelopoulos | Machine learning metabolic pathway descriptions using a probabilistic relational representation[END_REF][START_REF] Angelopoulos | Slps for probabilistic pathways: Modeling and parameter estimation[END_REF] or gene interactions [START_REF] Bernot | A fruitful application of formal methods to biological regulatory networks: Extending thomas' asynchronous logical approach with temporal logic[END_REF]. However learning biochemical reactions from temporal properties is quite new, both from the machine learning perspective and from the systems biology perspective. A precursor system of this type was the system KARDIO used in drug target discovery [START_REF] Bratko | KARDIO: A study in Deep and Qualitative Knowledge for Expert Systems[END_REF]. The novelty in our approach is the use of the temporal logic setting to express semi-qualitative semi-quantitative properties of the behaviour of the system to be captured by the model.

In the following, we present successively :

the boolean semantics of reaction models in Datalog, the representation of biological properties in temporal logic CTL, the application of ILP and model revision from temporal properties, the representation of kinetic models, of quantitative properties in temporal logic LTL with constraints, and a parameter search algorithm, the evaluation on an application: the modelling of the synchronization in period of the cell cycle by the circadian cycle.

2 Reaction rule learning from temporal properties

Biochemical reaction models in Datalog

From a syntactical point of view, SBML and BIOCHAM models basically consists in a set of reaction rules between molecules, protein complexes and modified proteins such as by phosphorylation. Each reaction rule for synthesis, degradation, complexation, phosphorylation, etc. can be given with a kinetic expression.

Example 1.

Here is for instance a simple model of the cell cycle control after Tyson (1991). Each rule is given here with an arithmetic expression (its rate) followed by the keyword for and then a list of reactants separated by + on the left side of the reaction arrow => and a list of products on the right side. The notation represents the empty list.

k1 for _=>Cyclin. k2*[Cyclin] for Cyclin=>_. k3*[Cyclin]*[Cdc2~{p1}] for Cyclin+Cdc2~{p1} => Cdc2~{p1}-Cyclin~{p1}. k4p*[Cdc2~{p1}-Cyclin~{p1}] for Cdc2~{p1}-Cyclin~{p1} => Cdc2-Cyclin~{p1}. k4*([Cdc2-Cyclin~{p1}])^2*[Cdc2~{p1}-Cyclin~{p1}] for Cdc2~{p1}-Cyclin~{p1} =[Cdc2-Cyclin~{p1}]=> Cdc2-Cyclin~{p1}. k5*[Cdc2-Cyclin~{p1}] for Cdc2-Cyclin~{p1} => Cdc2~{p1}-Cyclin~{p1}. k6*[Cdc2-Cyclin~{p1}] for Cdc2-Cyclin~{p1} => Cyclin~{p1}+Cdc2. k7*[Cyclin~{p1}] for Cyclin~{p1} =>_. k8*[Cdc2] for Cdc2 => Cdc2~{p1}. k9*[Cdc2~{p1}]
for Cdc2~{p1} => Cdc2.

The first rule represents the synthesis of a cyclin with a constant rate k1. The second rule represents the degradation of the cyclin with a reaction rate proportional to the cyclin concentration. The third rule represents the phosphorylation of the cyclin when it gets complexed with the kinase Cdc2~{p1}. The fourth rule is an autocatalyzed dephosphorylation of the complex, etc. For a more complete account of BIOCHAM syntax see for instance [START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF].

From a semantical point of view, reaction rules can be interpreted under different semantics corresponding to different abstraction levels. The most abstract semantics of BIOCHAM rules is the boolean semantics that associates to each molecule a boolean variable representing its presence or absence in the system, and ignores the kinetic expressions. Reaction rules are then interpreted as an asynchronous transition system2 over states defined by the vector of boolean variables. A rule such as A+B=>C+D defines four possible transitions corresponding to the complete or incomplete consumption of the reactants A and B. Such a rule can only be applied when both A and B are present in the current state. In the next state, C and D are then present, while A and B can either be present (partial consumption) or absent (complete consumption). Formally, the boolean semantics of a reaction model is a Kripke structure (see for instance [START_REF] Clarke | Model Checking[END_REF]) K = (S, R) where S is the set of states defined by the vector of boolean variables, and R ⊆ S × S is the transition relation between states, supposed to be total (i.e. ∀s ∈ S, ∃s ∈ S s.t. (s, s ) ∈ R). A path in K, starting from state s 0 is an infinite sequence of states π = s 0 , s 1 , • • • such that (s i , s i+1 ) ∈ R for all i ≥ 0. We denote by π k the path s k , s k+1 , • • •.

Biological properties in temporal logic CTL

In the boolean semantics of reaction models, the biological properties of interest are reachability properties, i.e. whether a particular protein can be produced from an initial state; checkpoints, i.e. whether a particular protein or state is compulsory to reach another state; stability, i.e. whether the system can (or will) always verify some property; etc.

Such properties can be expressed in the Computation Tree Logic CTL * [22] that is an extension of propositional logic for reasoning about an infinite tree of state transitions. CTL * uses operators about branches (non-deterministic choices) and time (state transitions). Two path quantifiers A and E are introduced to handle non-determinism: Aφ meaning that φ is true on all branches, and Eφ that it is true on at least one branch. The time operators are F, G, X, U and W ; Xφ meaning φ is true at the next transition, Gφ that φ is always true, F φ that φ is eventually true, φ U ψ meaning φ is always true until ψ becomes true, and φ W ψ meaning φ is either always true or until and when ψ becomes true.

, π |= ψ, s |=!ψ iff s |= ψ, s |= ψ & ψ iff s |= ψ and s |= ψ , s |= ψ | ψ iff s |= ψ or s |= ψ , s |= ψ ⇒ ψ iff s |= ψ or s |= ψ, π |= φ iff s |= φ where s is the first state of π, π |= Xψ iff π 1 |= ψ, π |= ψ U ψ iff there exists k ≥ 0 s.t. π k |= ψ and π j |= ψ for all 0 ≤ j < k. π |= ψ W ψ iff either for all k ≥ 0, π k |= ψ.
or there exists k ≥ 0 s.t. π k |= ψ&ψ and for all 0

≤ j < k, π j |= ψ. π |=!ψ iff π |= ψ, π |= ψ & ψ iff π |= ψ and π |= ψ , π |= ψ | ψ iff π |= ψ or π |= ψ , π |= ψ ⇒ ψ iff π |= ψ or π |= ψ,
Table 1. Inductive definition of the truth value of a CTL * formula in a state s or a path π, in a given Kripke structure K.

In this logic, F φ is equivalent to true U φ, Gφ to φ W f alse, and the following duality properties hold:

!(Eφ) = A(!φ), !(Xφ) = X(!φ), !(F φ) = G(!φ), !(φ U ψ) =!ψ W !φ and !(φ W ψ) =!ψ U !φ,
where ! denotes negation. The following abbreviation are used in BIOCHAM :

reachable(P) stands for EF (P ); steady(P) stands for EG(P ); stable(P) stands for AG(P ); checkpoint(Q,P) stands for !E(!Q U P ); oscillates(P) stands for EG((F !P ) ∧ (F P )).

These temporal properties can be checked in the Prolog representation of reaction rules, by using a symbolic model-checker written in Prolog. The BIOCHAM model checker in Prolog proceeds by computing both backward and forward frontiers of states, starting from the initial states (resp. the goal states) leading to a goal state (resp. an initial state). These sets of states are represented by Prolog facts with variables. Their cardinalities are reduced by subsumption checks in this representation. In its simplest form, the forward reachability analysis proceeds by computing the transitive closure of the transition relation, starting from the initial state, up to the reaching of a state in the query. The simplest case in such a model checker is thus a standard transitive closure algorithm in Prolog.

For performance reasons in large reaction models however, the symbolic model checker NuSMV [START_REF] Cimatti | Nusmv 2: An opensource tool for symbolic model checking[END_REF] based on ordered binary decision diagram (OBDD) is preferred and is used by default in BIOCHAM, through an interface. NuSMV is restricted to the fragment CTL of CTL * in which each time operator must be immediately preceded by a path quantifier. This restriction causes a difficulty for the oscillation properties only, since they cannot be expressed in CTL. In CTL, oscillation properties are thus approximated by the necessary but not sufficient formula EG((EF !P ) ∧ (EF P )). We refer to [START_REF] Chabrier | Symbolic model cheking of biochemical networks[END_REF][START_REF] Chabrier-Rivier | Modeling and querying biochemical interaction networks[END_REF] for the expressivity and scalability of this approach in reaction models containing several hundreds of variables and rules.

Model revision from temporal properties

Having the model and the properties defined by a Prolog program, ILP techniques can in principle be used for learning reaction rules from temporal properties, i.e. structure learning of the underlying logic program (see Chapter 1).

Here the positive and negative examples are uniformly given as a list of temporal properties to satisfy (expressed in a language closed by negation), instead of by positive and negative facts. Because of the relative complexity of the model checker in Prolog, this approach is currently limited to reachability properties. For learning from more general temporal properties, the NuSMV model checker is used in BIOCHAM as a black box, within an enumeration algorithm of all possible rule instances of some given rule pattern.

Furthermore, in the general framework of model revision, one wants to discover deletions as well as additions of reaction rules (of some pattern given as a bias) in order to satisfy a set of CTL formulas given as positive and negative examples. CTL properties can be classified into ECTL and ACTL formulas (i.e. formulas containing only E or A path quantifiers respectively) in order to anticipate whether reaction rules need be added or deleted. Indeed if an ECTL (resp. ACTL) formula is false in a Kripke structure, it remains false in a Kripke structure with less (resp. more) transitions. We refer to [START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF] for the details of the model revision algorithm implemented in BIOCHAM along these lines.

We show here our results on the model of example 1. For the structure learning phase, some CTL formulae are entered as a specification, expressing here reachability, oscillation and checkpoint properties:

add_specs({ reachable(Cdc2~{p1}), reachable(Cdc2), reachable(Cyclin), reachable(Cyclin~{p1}), reachable(Cdc2-Cyclin~{p1}), reachable(Cdc2~{p1}-Cyclin~{p1})}). add_specs({ oscil(Cdc2, oscil(Cdc2~{p1})), oscil(Cdc2~{p1}-Cyclin~{p1}), oscil(Cdc2-Cyclin~{p1}), oscil(Cyclin), checkpoint(Cdc2~{p1}-Cyclin~{p1}, Cdc2-Cyclin~{p1})}).
These properties are satisfied by the model and can be automatically checked by the model-checker. The simplest example to illustrate the structural learning method is to delete one rule in the model and let the learning system revise the model in order to satisfy the specification.

biocham: delete_rules(Cyclin+Cdc2~{p1}=>Cdc2~{p1}-Cyclin~{p1}). Cyclin+Cdc2~{p1}=>Cdc2~{p1}-Cyclin~{p1} biocham: check_all. The specification is not satisfied. This formula is the first not verified: Ai(oscil(Cdc2~{p1}-Cyclin~{p1})) biocham: revise_model(more_elementary_interaction_rules). Success Modifications found: Deletion(s):

Addition(s): Cyclin+Cdc2~{p1}=[Cdc2]=>Cdc2~{p1}-Cyclin~{p1}.
The first solution found is correct, even though it does not correspond to the deleted rule. In fact, there are four solutions consisting in adding one rule, the third one corresponds to the original model:

biocham: learn_one_addition(elementary_interaction_rules). (1) Cyclin+Cdc2~{p1}=[Cdc2]=>Cdc2~{p1}-Cyclin~{p1} (2) Cyclin+Cdc2~{p1}=[Cyclin]=>Cdc2~{p1}-Cyclin~{p1} (3) Cyclin+Cdc2~{p1}=>Cdc2~{p1}-Cyclin~{p1} (4) Cyclin+Cdc2~{p1}=[Cdc2~{p1}]=>Cdc2~{p1}-Cyclin~{p1}
It is worth noting that in these algorithms, the use of types [START_REF] Fages | Type inference in systems biology[END_REF] specifying the protein functions for instance, has the effect of reducing the number of possibilities and improving the performances in terms of both adequacy of results and computation time.

Parameter search from quantitative temporal properties

For relatively small networks of less than a hundred of proteins, kinetic models have been proved successful to perform quantitative analyses and predictions. Since the models of most datasets are in SBML, it is quite natural to handle the kinetic expressions provided in those models, especially for relating them to quantitative biological properties. In this section, we recal the two most usual semantics for those expressions, the differential semantics and the stochastic semantics, and relate them to PILP representations. We then show that the Linear Time Logic LTL with numerical constraints provides the expressive power necessary to represent both qualitative and quantitative properties of biological systems. Similarly to what is done in the boolean case, a model-checker is then used as basis for a learning process allowing here to find parameter values fitting a given LTL specification of the biological properties that the model is supposed to reproduce. This is shown on example 1 and is developed in an application in the next section.

Continuous semantics with ODE's

The concentration semantics of BIOCHAM associates to each molecule a real number representing its concentration. Reaction rules are in fact interpreted with their kinetic expressions by a set of nonlinear ordinary differential equations (ODE) 3 . Formally, to a set of BIOCHAM reaction rules E = {e i for S i => S i } i=1,...,n with variables {x 1 , ..., x m }, one associates the system of ODEs :

dx k /dt = n i=1 r i (x k ) * e i - n j=1 l j (x k ) * e j
where r i (x k ) (resp. l i ) is the stoichiometric coefficient of x k in the right (resp. left) member of rule i. Given an initial state, i.e. initial concentrations for each of the objects, the evolution of the system is deterministic, and numerical integration algorithms compute a time series describing the temporal evolution of the system variables. The integration methods actually implemented in BIOCHAM are the adaptive step-size Runge-Kutta method and the Rosenbrock implicit method for stiff systems, which both produce simulation traces with variable time steps and are implemented in Prolog.

Stochastic semantics with SLPs

The stochastic semantics is the most realistic semantics but also the most difficult to compute. This semantics associates to each BIOCHAM object an integer representing the number of molecules in the system. Rules are interpreted as a continuous time Markov chain where transition probabilities are defined by the kinetic expressions of reaction rules.

Stochastic simulation techniques [START_REF] Gillespie | General method for numerically simulating stochastic time evolution of coupled chemical-reactions[END_REF] compute realizations of the process. The results are generally noisy versions of those obtained with the concentration semantics. However, in models with, for instance, very few molecules of some kind, qualitatively different behaviors may appear in the stochastic simulation, and thus justify the recourse to that semantics in such cases. A classical example is the model of the lambda phage virus [START_REF] Gibson | A probabilistic model of a prokaryotic gene and its regulation[END_REF] in which a small number of molecules, promotion factors of two genes, can generate an explosive multiplication (lysis) after a more or less long period of passive wait (lysogeny).

In the stochastic semantics, for a given volume V of the location where a compound is situated, its concentration C is translated into a number of molecules N = C × V × K, where K is Avogadro's number. The kinetic expression e i for the reaction i is converted into a transition rate τ i by replacing all concentrations by the corresponding number of molecules multiplied by volume. After normalization on all possible transitions, this gives the transition probability

p i = τi n j=1 τj .
This semantics is close to SLPs. Two points however render unusable the classical learning techniques, and suggest an extension of the SLP framework :

kinetic expressions, and thus the corresponding transition probabilities τ i , can contain variables representing the molecular concentrations (resp. number) of the reactants in each rule. A faithful translation of those models into SLP would thus involve dynamic probabilities according to variables values, like in the stochastic semantics of BIOCHAM by continuous time Markov chains [START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF]. On the other hand, SLPs as defined in Section 3 of Chapter 2, are restricted to constant probabilities on each rule. in stochastic simulation and Gillespie algorithms [START_REF] Gillespie | General method for numerically simulating stochastic time evolution of coupled chemical-reactions[END_REF], the time is a random variable over reals, which cannot be mixed with SLPs in the current version of the formalism.

Biological properties in LTL with numerical constraints

The Linear Time Logic, LTL is the fragment of CTL * that uses only temporal operators. A first-order version of LTL is used to express temporal properties about the molecular concentrations in the simulation trace. A similar approach is used in the DARPA BioSpice project [START_REF] Antoniotti | Model building and model checking for biochemical processes[END_REF]. The choice of LTL is motivated by the fact that the concentration semantics given by ODEs is deterministic, and there is thus no point in considering path quantifiers. The version of LTL with arithmetic constraints we use, considers first-order atomic formulae with equality, inequality and arithmetic operators ranging over real values of concentrations and of their derivatives.

For instance F([A]>10) expresses that the concentration of A eventually gets above the threshold value 10. G([A]+[B]<[C]

) expresses that the concentration of C is always greater than the sum of the concentrations of A and B. Oscillation properties, abbreviated as oscil(M,K), are defined as a change of sign of the derivative of M at least K times:

F((d[M]/dt>0) & F((d[M]/dt<0) & F((d[M]/dt>0)...))
). The abbreviated formula oscil(M,K,V) adds the constraint that the maximum concentration of M must be above the threshold V in at least K oscillations.

For practical purposes, some limited forms of quantified first-order LTL formulae are also allowed. As an example of this, constraints on the periods of oscillations can be expressed with a formula such as period(A,75), defined as

∃t ∃v F (T ime = t & [A] = v & d([A])/dt > 0 & X(d([A])/dt < 0) & F (T ime = t + 75 & [A] = v & d([A])/dt > 0 & X(d([A])/dt < 0)))
where T ime is the time variable. This very formula is used extensively in the example of the next section.

Note that the notion of next state (operator X) refers to the state of the following time point computed by the (variable step-size) simulation, and thus does not necessarily imply real-time neighborhood. Nevertheless, for computing local maxima as in the formula above for instance, the numerical integration methods do compute the relevant time points with a very good accuracy.

Parameter search from temporal properties

We implemented a dedicated LTL model checker for biochemical properties over simulation traces and proceeded, as in the boolean case, to use it for a learning method. Actually, it is mostly a search method automatically evaluating the fitness of a given parameter set w.r.t. an LTL specification.

The same method could theoretically be used to sample for a probability of satisfaction of an LTL specification for the stochastic semantics, however experimental trials proved to be too computationally expensive. The parameter learning method can be illustrated by changing the value of some parameter like k4 for instance. Figure 1 shows that the model is not oscillating as it should when k4 is set to zero.

The property of oscillation can be added as a temporal logic constraint to the parameter value search system as follows:

learn_parameters([k4],[(0,200)],20,oscil(Cdc2-Cyclin~{p1},3),100). First values found that make oscil(Cdc2-Cyclin~{p1},3) true: parameter(k4,200).

The value 200 found for k4 is close to the original value (180) and satisfies the experimental results formalized in LTL, as depicted in Figure 2.

Note that beacuase of the highly non-linear nature of the kinetics used in most biological models of the literature it is not possible to rely on usual tools of control theory for this kind of parameter estimation. The other available techniques are mostly local optimization based (simulated annealing and derivatives) but require to optimize with respect to a precise quantitative objective function, whereas the presented technique allows to mix qualitative and quantitative data.

Application to modelling the synchronization in period of the cell cycle by the circadian cycle

Cancer treatments based on the administration of medicines at different times of the day have been shown to be more efficient against malign cells and less damaging towards healthy ones. These results might be related to the recent discovery of links between the circadian clock (controlled by the light/dark cycle of a day) and the cell cycle. However, if many models have been developed to describe both of these cycles, to our knowledge none has described a real interaction between them.

In the perspective of the European Union project TEMPO4 on temporal genomics for patient tailored chronotherapies, we developed a coupled model at the molecular level and studied the conditions of synchronization in period of these cycles, by using the parameter learning features of the modeling environment BIOCHAM. More specifically, the learning of parameter values from temporal properties with numerical constraints has been used to search how and where in the parameter space of our model the two cycles get synchronized. The technical report [START_REF] Calzone | Coupling the cell cycle and the circadian cycle[END_REF] describes the conditions of synchronization (i.e. synchronization by forcing the period of the target oscillator to be the same as that of the forcing oscillator) of the cell cycle by the circadian cycle via a common protein kinase WEE1 (see Figure 3).

The proteins chosen to illustrate the cell cycle are MPF, preMPF, the degradation factor of the cyclins, APC, the WEE1 kinase and the CDC25 phosphatase (Figure 3). Early in the cycle, MPF is kept inactive because the cyclin is not synthesized and WEE1 is present. As the cyclin is slowly synthesized, MPF activates and reaches a threshold that both inactivates WEE1 and activates CDC25 which maintains MPF in its active state. The cell enters mitosis. With a short delay, APC is activated and degrades the cyclin component of MPF. The cell exits mitosis and repeats its cycle. The model is composed of two positive feedback loops (CDC25 activates MPF which in turn activates CDC25, and WEE1 inactivates MPF which in turn inactivates WEE1) and a negative feedback loop (MPF activates APC through an intermediary enzyme X and APC degrades the cyclin component of the complex MPF). See Figure 4.

The two models describing the cell and circadian cycles are linked through the transcription of WEE1. In the model of the cell cycle alone, wee1 mRNA was a parameter equal to 1. In the coupled model, the production of Wee1m is a function of the nuclear form of the complex BMAL1/CLOCK (BN) and the unphosphorylated nuclear form of the complex PER/CRY (PCN).

To find values for which synchronization occurs, the parameter space for each parameter has been explored using the BIOCHAM learning features from temporal properties. The values of three parameters appeared to be more significant than others: ksweem, kswee and kimpf. The two parameters ksweem and kswee both control the level of the WEE1 protein and show such similarities that in the following discussion, we will only report on kswee. The parameter values are varied in a given interval and reveal domains of synchronization reported in Fig- For intermediate values of the parameters (region 2), WEE1 starts to play a more significant role in the cell cycle by inhibiting MPF activity, and as a result, disturbing MPF oscillations. It is only when the parameters reach a high value (either kimpf=1.2 or kswee=0.4) that the oscillations of MPF become stable again but with a period similar to that of the circadian cycle (region 3) revealing the synchronization of the cell cycle through WEE1 activity (through kimpf) or protein level (through kswee).

However, the study of kimpf, the parameter controlling the activity of WEE1 on MPF inactivation, shows that the synchronization does not solely depend on the value of the parameter but more particularly on the ratio kimpf/kampf since both CDC25 and WEE1 are involved in the positive feedback loops that activate MPF and therefore responsible for the G2-M transition. To investigate this dual effect, the limit of synchronization is measured as the two parameters kimpf and kampf are varied simultaneously. A linear function of the form: kampf = 2.44832• kimpf + 2.0071 is obtained, the region below the line being the synchronization region (Figure 6).

These studies have been carried out thanks to the machine learning features of BIOCHAM to express the condition of synchronization in period by a first-order LTL formula with time constraints, and to explore the parameter values satisfying that formula. These methods are complementary to and have no counterpart in traditional tools for the study of dynamical systems such as the computation of bifurcation diagrams.

Discussion and conclusions

Temporal logic is a powerful formalism for expressing the biological properties of a living system, such as state reachability, checkpoints, stability, oscillations, etc. This can be done both qualitatively and quantitatively, by considering respectively propositional and first-order temporal logic formulae with numerical constraints.

In the propositional case, we have given a Datalog representation of general biochemical reaction models allowing the use of ILP techniques for discovering reaction rules from given CTL properties. Because of the relative complexity of the CTL model checker however, this approach was limited to reachability properties. For general CTL properties, our approach was to use the symbolic OBDD model checker NuSMV as a black box within a model revision algorithm that searches for rule additions and deletions in order to satisfy a CTL specification. The first results are encouraging but also show some limitations concerning the boolean abstraction that simply forgets the kinetic expressions. Less crude abstractions are however possible and are currently under investigation. Furthermore, when restricting to pure reachability properties between completely defined states, ILP methods have been shown more efficient. There is thus a perspective for combining the best of both methods in this general setting.

Kinetic models of biochemical systems have been considered too with their two most usual interpretations, by ODEs, and by continuous time Markov chains. The second interpretation has been related to PILP representations with a generalized notion of SLPs involving dynamic probabilities according to variable values. However the stochastic interpretation of kinetic expressions is computationally too expensive, while the interpretation by differential equations does scale up to real-size quantitative models. This is the reason why the continuous semantics of BIOCHAM rules based on non-linear ordinary differential equations, instead of the stochastic semantics based on continuous-time Markov chains is used in these applications. We have shown that the inference of parameter values from a temporal logic specification is flexible enough to be easy to use, and provides accurate results with reasonable execution times. These functionalities are completely new and complement the other tools the modeler can use to estimate the range of parameter values. This has been illustrated by an original application to the modelling of the synchronization in period of the cell cycle by the circadian cycle. Fig. 5. Plot of the period as a function of the parameter kimpf from 0.01 to 1.6 and kswee from 0.01 to 1. The system shows synchronization for values superior to 1.2 for kimpf and 0.4 for kswee. For our purposes, constant periods are defined as follows: the last 11 peaks of the simulation over 500 time units show no more than 4% difference in their maxima and the length of the periods. MPF starts with an autonomous period of 22.4h and BN a period of 23.85h. As kimpf and kswee increase, MPF oscillations (accounting for cell cycle) lose stability and are entrained, along with WEE1, for higher values of the parameter with a period of 23.85h.

Fig. 1 .

 1 Fig. 1. Broken cell cycle model with parameter k4 = 0

Fig. 2 .

 2 Fig. 2. Curated cell cycle model showing oscillations with the inferred parameter value k4 = 200.

Fig. 3 .

 3 Fig. 3. Linking the circadian and the cell cycles via WEE1

Fig. 4 .

 4 Fig. 4. Temporal simulation of a generic cell cycle

Fig. 6 .

 6 Fig. 6. BIOCHAM-generated plot of the synchronization in period of the cell cycle by the circadian cycle for different values of kimpf (action of WEE1 on MPF) and kampf (action of CDC25 on MPF). The limit of synchronization computed by BIOCHAM (red crosses) is interpolated by the linear function kampf = 2.44832 • kimpf + 2.0071 (solid line).

  Table 1 recalls the truth value of a formula in a given Kripke structure.

s |= α iff α is a propositional formula true in the state s, s |= Eψ iff there exists a path π starting from s s.t. π |= ψ, s |= Aψ iff for all paths π starting from s

BIOCHAM is available for download at http://contraintes.inria.fr/BIOCHAM

In this context asynchronous refers to the fact that only one transition is fired at a time, even if several are possible. This choice is justified by the fundamental biochemical phenomena of competition and masking between reaction rules.

The kinetic expressions in BIOCHAM can actually contain conditional expressions, in which case the reaction rules are interpreted by a deterministic hybrid automaton.

http://www.chrono-tempo.org

Acknowledgements. The authors would like to thank Stephen Muggleton and Luc de Raedt for extensive discussions on the topics of this paper, and Nathalie Chabrier-Rivier and Laurence Calzone for their contributions. This work has been partly supported by the EC Sixth Framework Project Application of Probabilistic Inductive Logic Programming II (APrIL II) (Grant Ref: FP-508861).