Informations

Équation de mesure du modèle de parcours Tree-Step

Équation de mesure du modèle de parcours Tree-Step


We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

Je ne suis pas spécialiste en biologie mais je souhaite tester une méthode d'estimation des paramètres (basée sur le filtre de Kalman) sur le modèle Tree-Step pathway, mais je ne trouve pas de document qui donne l'équation de mesure du modèle Three-Step Pathway (Moles et all 2003, Parameter Estimation in Biochemical Pathways: A Comparison of Global Optimization Methods"), quelqu'un peut-il me donner quelques références ?


Modèles en biologie : « des descriptions précises de notre pensée pathétique »

Dans cet essai, je vais esquisser quelques idées sur la façon de penser aux modèles en biologie. Je commencerai par essayer de dissiper le mythe selon lequel la modélisation quantitative est en quelque sorte étrangère à la biologie. Je soulignerai ensuite la distinction entre modélisation directe et rétro-modélisation et me concentrerai par la suite sur la première. Au lieu d'entrer dans des détails techniques mathématiques sur différentes variétés de modèles, je me concentrerai sur leur structure logique, en termes d'hypothèses et de conclusions. Un modèle est une machine logique pour déduire le second du premier. Si le modèle est correct, alors, si vous croyez à ses hypothèses, vous devez, en toute logique, croire également à ses conclusions. Cela conduit à considérer les hypothèses qui sous-tendent les modèles. Si celles-ci sont basées sur des lois physiques fondamentales, alors il peut être raisonnable de traiter le modèle comme « prédictif », en ce sens qu'il n'est pas sujet à falsification et que l'on peut se fier à ses conclusions. Cependant, au niveau moléculaire, les modèles sont plus souvent dérivés de la phénoménologie et de conjectures. Dans ce cas, le modèle est un test de ses hypothèses et doit être falsifiable. Je discuterai de trois modèles de cette perspective, chacun d'eux donnant des informations biologiques, et cela conduira à quelques lignes directrices pour les constructeurs de modèles potentiels.


Résumé de l'auteur

L'estimation des paramètres est une question clé en biologie des systèmes, car elle représente l'étape cruciale pour obtenir des prédictions à partir de modèles informatiques de systèmes biologiques. Ce problème est généralement résolu en « adaptant » les simulations du modèle aux données expérimentales observées. Une telle approche ne prend pas pleinement en considération le bruit de mesure. Nous introduisons une nouvelle méthode basée sur la combinaison du filtrage de Kalman, des tests statistiques et des techniques d'optimisation. Le filtre est bien connu dans la théorie du contrôle et de l'estimation et a trouvé des applications dans un large éventail de domaines, tels que les systèmes de guidage inertiel, les prévisions météorologiques et l'économie. Nous montrons comment les statistiques du bruit de mesure peuvent être exploitées de manière optimale et directement intégrées dans la conception de l'algorithme d'estimation afin d'obtenir des résultats plus précis et de valider/invalider les estimations calculées. Nous montrons également qu'un avantage significatif de notre estimateur est qu'il offre un outil puissant pour la sélection de modèles, permettant le rejet ou l'acceptation de modèles concurrents sur la base des mesures bruitées disponibles. Ces résultats ont une application pratique immédiate en biologie computationnelle, et bien que nous démontrions leur utilisation pour deux exemples spécifiques, ils peuvent en fait être utilisés pour étudier une large classe de systèmes biologiques.

Citation: Lillacci G, Khammash M (2010) Estimation des paramètres et sélection de modèles en biologie computationnelle. PLoS Comput Biol 6(3) : e1000696. https://doi.org/10.1371/journal.pcbi.1000696

Éditeur: Anand R. Asthagiri, California Institute of Technology, États-Unis d'Amérique

A reçu: 22 janvier 2009 Accepté: 30 janvier 2010 Publié : 5 mars 2010

Droits d'auteur: © 2010 Lillacci, Khammash. Il s'agit d'un article en libre accès distribué selon les termes de la licence d'attribution Creative Commons, qui permet une utilisation, une distribution et une reproduction sans restriction sur n'importe quel support, à condition que l'auteur et la source d'origine soient crédités.

Le financement: Ce matériel est basé sur des travaux soutenus par la National Science Foundation sous Grant NSF-CDI ECCS-0835847, NSF-ECCS-0802008, NIH R01GM049831, et l'Institute for Collaborative Biotechnologies par Grant DAAD19-03-D-0004 de l'US Army Research Bureau. Les bailleurs de fonds n'ont joué aucun rôle dans la conception de l'étude, la collecte et l'analyse des données, la décision de publier ou la préparation du manuscrit.

Intérêts concurrents : Les auteurs ont déclaré qu'ils n'existaient pas de conflit d'intérêts.


Résultats

Le modèle nominal

Nous supposons qu'un modèle ODE nominal

a été proposé pour décrire la dynamique du système de considération. Le vecteur d'état contient le m variables dynamiques et est la dérivée par rapport au temps t. La valeur initiale du vecteur d'état est . Pour un réseau de réaction biochimique, c'est souvent la concentration ou l'abondance des k-ème espèce. La fonction représente une entrée externe connue du système. La dynamique des variables d'état est déterminée par la fonction et code les hypothèses de modèle faites dans le modèle nominal. Cela peut être représenté par un graphique 20 , où chaque nœud correspond à une variable et un bord dirigé de je à k indique que la dérivée temporelle de dépend de (Fig. 1a). Si est directement influencé par une entrée connue, nous l'illustrons par une flèche verte en zigzag. En règle générale, toutes les variables d'état ne peuvent pas être mesurées directement. Les variables représentent toutes les sorties qui sont expérimentalement accessibles. Dans l'équation (1b), nous supposons que l'application h de l'état X à la sortie oui est connu. Nous utilisons un tilde pour souligner cela et ne sont donc généralement pas parfaitement connus en raison de connaissances limitées ou incertaines sur la véritable dynamique sous-jacente.

Représentation de l'erreur de modèle

La réponse du système naturel réel à un stimulus d'entrée connu est généralement mesurée à des moments discrets et fournit des observations expérimentales pour la sortie. Une partie de ces données est généralement utilisée pour estimer les paramètres du modèle. Nous considérons les estimations initiales des paramètres comme faisant partie de la spécification du modèle nominal dans l'équation (1a).

Le modèle nominal n'est pas satisfaisant lorsque sa sortie n'est pas en accord suffisant avec les données . Une source d'erreur de modèle provient des entrées cachées du système nominal, qui sont causées par des processus dynamiques exogènes au système nominal (Fig. 1b). De plus, il peut y avoir des interactions manquantes ou erronées entre les variables d'état dans le modèle nominal lui-même. Les deux types d'erreur de modèle peuvent être représentés par des entrées cachées agissant sur les nœuds du modèle nominal (Fig. 1c). La « vraie » dynamique du système réel peut être décrite par

Ici, l'état représente les mêmes variables que l'état nominal , mais nous supprimons le tilde pour distinguer les solutions de (2) de celle du modèle nominal. L'erreur du modèle est la différence entre le taux de changement du système réel et le système nominal, évalué le long de la trajectoire de l'état réel. Ainsi, il intègre tout écart entre le système réel et le système nominal. L'entrée connue vous et la fonction de sortie h sont supposés identiques au modèle nominal (1). Cependant, nous discuterons également de l'impact du bruit de mesure ci-dessous.

L'approche typique de l'amélioration du modèle consiste à compenser l'erreur du modèle par des expressions mathématiques explicites, souvent des équations différentielles supplémentaires. Cela augmente le nombre de variables et de paramètres dans le modèle. Ici, on procède différemment en estimant l'erreur de modèle w à partir des données, ce qui permet aussi de corriger le biais de l'estimation d'état induit par le modèle nominal.

Estimation de la dynamique non modélisée

Pour estimer l'erreur du modèle, nous utilisons le système d'observateur

qui est une copie des équations (2a) et (2b). Le chapeau marque les estimations de l'état, de la sortie et de l'erreur de modèle. Cette dernière est obtenue en minimisant la fonctionnelle d'erreur

Le premier terme de l'équation (3c) est l'erreur quadratique moyenne pondérée entre les sorties mesurées et les sorties du système d'observation dans les équations (3a) et (3b). La norme du carré pondéré

contient la matrice de pondération symétrique , qui est souvent choisie pour être diagonale et peut être utilisée pour transformer des sorties de magnitude très différente en une échelle commune ou pour incorporer des estimations de précision des mesures à différents moments . Le terme de régularisation

est nécessaire pour éviter le surajustement des données par des estimations trop complexes . Les paramètres non négatifs et déterminent les contributions relatives de la norme dans l'équation (3f) et de la norme dans (3g). La minimisation de l'équation (3c) sous les contraintes des équations (3a) et (3b) est un problème de contrôle optimal 21,22,23 , qui doit être résolu numériquement (voir Méthodes et texte supplémentaire).

La régularisation combinée dans l'équation (3e) rappelle la pénalité de filet élastique utilisée dans les modèles de régression 18 . Par conséquent, nous avons appelé notre approche le réseau élastique dynamique. Par analogie avec la régression, le terme entraîne la réduction à zéro de certaines composantes de l'erreur de modèle estimée (texte supplémentaire). La quantité de rétrécissement est déterminée par , qui peut être choisie pour supprimer les petits signaux d'erreur ou le bruit répartis sur de nombreux composants de l'estimation. L'estimation clairsemée résultante est utile, car elle fournit des informations sur les états du système qui sont ciblés par des erreurs de modèle systématiques, représentées par des entrées cachées.

Contrairement à la régression, une régularisation pure ou de type Lasso 24 n'est pas utile dans le cadre dynamique, car la solution pour peut entraîner des estimations non bornées de . Même lorsque des contraintes supplémentaires sur sont imposées, la solution résultante n'est pas lisse et nulle ou aux frontières des contraintes 25 . Ces informations sur le problème de contrôle optimal peuvent être obtenues à partir du principe minimum de Pontryagin 21,22 , tel qu'il est détaillé dans le texte supplémentaire avec quelques stratégies pour choisir les paramètres de régularisation appropriés et . En plus des estimations clairsemées mais régulières de l'erreur du modèle, le réseau élastique dynamique fournit automatiquement une estimation de l'état . Il s'agit souvent d'informations très intéressantes, lorsque toutes les variables d'état ne sont pas accessibles expérimentalement.

Le problème de contrôle optimal dans les équations (3a–c) pour nécessite la spécification d'une condition initiale , qui est souvent inconnue ou incertaine. Alternativement, on peut ajouter la contrainte supplémentaire

à (3a–c), où est une tolérance prédéfinie donnée pour l'ajustement de à au temps . De même, une tolérance peut être prescrite pour l'ajustement au dernier point de données en

Les paramètres de tolérance et de ces contraintes optionnelles peuvent souvent être obtenus à partir des barres d'erreur des mesures.

Validation du filet élastique dynamique

Exemple de signalisation JAK-STAT

Pour illustrer l'estimateur dynamique élastique-net pour un petit modèle compréhensible, nous avons utilisé des données expérimentales établies pour la voie de transduction du signal JAK-STAT 4 . Les quatre variables d'état du système représentent la STAT5 cytoplasmique non phosphorylée, la STAT5 monomère phosphorylée, la STAT5 dimère phosphorylée et la STAT5 dimère nucléaire. Le modèle nominal 4

décrit la phosphorylation de STAT5 cytoplasmique lors de l'activation du récepteur de l'érythropoïétine (entrée connue vous), la dimérisation de STAT5 phosphorylée et l'export vers le noyau (Fig. 2). Les données d'évolution temporelle 4 pour la quantité de STAT5 phosphorylée cytoplasmique et de STAT5 cytoplasmique totale ont été utilisées pour calibrer les paramètres. Cependant, la présence d'une erreur de modèle systématique est apparente à partir de l'écart inaltérable entre les données expérimentales et le modèle nominal incorporant des valeurs de paramètres optimisées (Fig. 2b,c).

Estimation de l'erreur de modèle pour la voie JAK-STAT.

(une) L'entrée connue est donnée par des mesures de phosphorylation interpolées linéairement pour le récepteur de l'érythropoïétine 4 . (b,c) Les mesures de sortie 4 (noir) pour STAT5 phosphorylée et STAT5 total dans le cytoplasme par rapport aux sorties du modèle nominal (bleu) et l'ajustement du filet élastique dynamique (rouge). () Graphique du modèle nominal (bleu) et du système observateur (rouge) avec les variables d'état STAT5 cytoplasmique, STAT5 monomère phosphorylé, STAT5 dimère phosphorylé et STAT5 nucléaire. (e) Estimations dynamiques du réseau élastique de l'erreur du modèle et de l'aire sous la courbe (AUC) pour l'amplitude de chaque composante . (F) Les estimations d'état obtenues à partir du modèle nominal (bleu) et de l'observateur dynamique à réseau élastique (voir en rouge).

Pour estimer cette erreur de modèle, nous avons ajusté numériquement le réseau élastique dynamique (3) avec le modèle nominal (4) aux mesures de sortie. Pour quantifier l'amplitude des différentes composantes, nous avons calculé numériquement l'aire sous la courbe (AUC) de chacun, c'est-à-dire . L'AUC et l'évolution temporelle estimée de l'erreur du modèle indiquent (Fig. 2e), que les contributions dominantes et de l'erreur du modèle ciblent les états et , représentant la quantité de STAT5 cytosplasmique non phosphorylée et de STAT5 nucléaire. La seconde composante de l'estimation dynamique du réseau élastique est identiquement nulle pour tout l'intervalle de temps (Fig. 2e). Mis à part le petit signal initié après environ 40 minutes, cela est cohérent avec le modèle de cycle nucléocytoplasmique amélioré rapporté dans 26, qui est basé sur les mêmes données 4 et intègre la relocalisation des molécules nucléaires STAT5 déphosphorylées dans le cytoplasme. Il est important de noter que le réseau élastique dynamique fournit également des estimations modifiées pour les quatre variables d'état STAT5 (Fig. 2f), qui sont également en bon accord avec le modèle de cycle nucléocytoplasmique (texte supplémentaire).

Un problème important avec les approches de régularisation est le choix des paramètres de régularisation et . Nous avons utilisé et dans la Fig. 2, mais nous avons trouvé empiriquement que les valeurs AUC indiquent clairement les points cibles de l'erreur du modèle pour une large gamme de valeurs (Fig. supplémentaire S2). Le paramètre a été choisi pour équilibrer la régularité et la précision de l'ajustement aux mesures de sortie. De plus, le biais induit par la double régularisation 18 peut être compensé par une stratégie de seuillage simple : étant donné une estimation initiale de l'erreur du modèle, nous réajustons le réseau élastique dynamique en contraignant toutes les composantes avec une petite AUC à zéro. Le seuil est connu dans le contexte de la régression 27 et nous avons constaté qu'il améliore les estimations d'état ainsi que les estimations de l'évolution dans le temps des erreurs de modèle restantes (Fig. supplémentaire S3).

L'impact du bruit de mesure et les incertitudes des paramètres

Pour explorer la robustesse du réseau élastique dynamique contre le bruit de mesure, nous avons ajouté des perturbations aléatoires aux données expérimentales 4 . Pour un niveau de bruit donné, nous avons généré 500 ensembles de données perturbées en ajoutant des nombres aléatoires gaussiens avec une moyenne de zéro et un écart type mis à l'échelle par un multiple de l'écart type empirique (voir les barres d'erreur sur les Fig. 2b, c) à chaque point de données expérimentales. Ainsi, le niveau de bruit est défini comme un multiple de l'écart type empirique. Le réseau élastique dynamique a ensuite été ajusté à chaque échantillon de sortie et l'aire correspondante sous la courbe pour chaque composante de l'erreur de modèle estimée a été calculée. Les tracés de ces valeurs AUC en fonction du niveau de bruit sont illustrés à la figure 3a. Les valeurs médianes de l'AUC pour les composants sont largement indépendantes du niveau de bruit, mais la variabilité des estimations de l'AUC augmente avec le bruit de mesure. Néanmoins, les valeurs d'AUC pour et sont toujours bien supérieures à zéro, alors que l'AUC de et est proche ou même égale à zéro pour de nombreux échantillons avec un niveau de bruit plus élevé. Cela augmente la confiance que les nœuds et (Fig. 2d) du modèle nominal JAK-STAT (4) sont les principaux points cibles de l'erreur du modèle.

L'impact du bruit de mesure simulé et de l'incertitude des paramètres sur l'estimation dynamique du réseau élastique dans le modèle JAK-STAT.

(une) Boîtes à moustaches visualisant la variation de l'AUC de pour les estimations de filet élastique dynamique causées par le bruit de mesure (voir le texte principal pour plus de détails). Pour faciliter la visualisation, les box plots à un niveau de bruit donné sont légèrement décalés. (b) La variation de l'AUC causée par l'incertitude des paramètres.

L'impact de l'incertitude des paramètres dans le modèle nominal a été évalué de manière similaire. Les algorithmes d'estimation des paramètres 4,10,26 appliqués au modèle nominal en utilisant les données expérimentales (Fig. 2b,c) fournissent des estimations ponctuelles et des intervalles de confiance pour chaque composante du vecteur de paramètres. Ces intervalles de confiance ont été à nouveau mis à l'échelle par le niveau de bruit, donnant un intervalle pour chaque paramètre à partir duquel des échantillons aléatoires uniformes ont été tirés. Encore une fois, nous avons généré 500 vecteurs de paramètres modifiés par niveau de bruit. Pour chaque échantillon de paramètres, le système (4) a été pris comme modèle nominal et l'AUC des estimations résultantes a été enregistrée (Fig. 3b). Encore une fois, il n'y a pas de tendance systématique pour l'ASC des différentes composantes de l'erreur estimée . Cependant, la variation de l'AUC augmente beaucoup plus rapidement que sur la figure 3a. Outre les différentes distributions d'échantillonnage utilisées, cet effet est lié à la définition de l'erreur de modèle w, qui est toujours défini par rapport au modèle nominal (conférer l'équation 2a). Par conséquent, l'erreur de modèle estimée contient des contributions provenant à la fois de spécifications erronées structurelles et de paramètres dans le modèle nominal. Néanmoins, il est toujours possible d'inférer les composantes dominantes et avec un niveau de confiance élevé. Des résultats similaires ont été trouvés pour la sensibilité par rapport au nombre de points de temps de mesure (texte supplémentaire, Fig. S6).

Exemple de signalisation UV-B photomorphogène

Comme cas de test pour un système plus vaste, nous avons utilisé un modèle récent pour la coordination de la signalisation UV-B photomorphogène chez les plantes 19 . Le modèle est constitué de 11 ODE décrivant la dynamique des concentrations de protéines couplées par 10 réactions chimiques (Fig. 4). Nous avons considéré ce modèle comme le modèle nominal afin de tester la méthode du réseau élastique dynamique pour une situation où la vérité terrain est connue. L'erreur du modèle a été simulée en ajoutant les entrées cachées aux nœuds et . La fonction de sortie est une combinaison linéaire de 7 variables d'état différentes (voir le texte supplémentaire pour toutes les équations). Des données synthétiques ont été échantillonnées à des moments discrets à partir des sorties du modèle réel et des perturbations aléatoires gaussiennes ont été ajoutées pour simuler le bruit de mesure (Fig. 4b–f). Le réseau élastique dynamique avec le modèle nominal a été utilisé pour reconstruire l'erreur du modèle et l'état réel à partir de ces données simulées. L'aire absolue sous la courbe pour chaque composante de l'estimation d'erreur du modèle indique clairement que les états et sont ciblés par des entrées cachées (Fig. 4g), alors que toutes les autres composantes sont soit très petites, soit même nulles. Cela illustre la rareté de l'estimation dynamique du réseau élastique, qui est un net avantage par rapport L2 régularisation. L'écart entre l'erreur du modèle et l'estimation correspondante par rapport à l'amplitude de la vraie erreur du modèle est d'au plus 10 % (Fig. 4h) et est principalement causé par des imprécisions numériques. Plus important encore, l'écart entre la trajectoire d'état réelle et estimée est presque nul (Fig. 4i), indiquant l'excellente performance du réseau élastique dynamique en tant qu'observateur d'état.

L'exemple de la signalisation UV-B photomorphogénique.

(une) Le graphique (sans auto-boucles) des états du modèle 19 . Les points cibles des erreurs du modèle simulé sont indiqués par les flèches rouges. (bF) La sortie simulée avec des barres d'erreur (noir), la sortie du modèle nominal (bleu) et la sortie du réseau élastique dynamique (rouge). (g) L'AUC des erreurs absolues du modèle. (h) Les composantes de par rapport à l'amplitude UNE de la véritable erreur de modèle. (je) L'écart entre l'état réel et l'état nominal du modèle (bleu) par rapport à l'écart du réseau élastique dynamique (rouge).

Tester les limites

Comme pour toute méthode inverse, il existe des limites à la méthode dynamique du réseau élastique. Certaines erreurs de modèle ne sont pas observables, car il existe une fonction d'entrée cachée différente qui génère une sortie identique à la sortie obtenue pour , voir le texte supplémentaire pour un exemple simple. D'autres erreurs de modèle peuvent être pratiquement inobservables, car la sortie d'une autre fonction d'entrée cachée peut ne pas être distinguable dans les erreurs de mesure. Un cas particulier est celui des erreurs de modèle qui n'ont aucun ou presque aucun effet sur la sortie. Ceux-ci ne seront pas remarqués lors de la modélisation et le modèle nominal sera accepté.

Pour tester davantage la capacité du réseau élastique dynamique à déduire les états ciblés par l'erreur de modèle, c'est-à-dire les composantes non nulles de la véritable erreur de modèle, nous avons systématiquement simulé des perturbations sur différents nœuds et paires de nœuds. Tout d'abord, nous avons simulé des erreurs de modèle ciblant un seul nœud k de la même manière qu'avant. Pour les nœuds et il n'y a eu aucun effet sur la sortie (voir à nouveau Fig. 4b–f) et donc ces nœuds ont été omis d'une analyse plus approfondie. De plus, nous avons simulé des entrées cachées pour toutes les combinaisons de nœuds restantes. Pour chacun de ces 36 modèles réels simulés, nous avons testé la capacité du réseau élastique dynamique à récupérer les nœuds cibles corrects à partir de l'AUC de l'estimation. Nous avons considéré qu'un nœud ou une paire de nœuds était correctement récupéré, si son AUC était d'au moins 85 % de l'AUC totale sur tous les nœuds. Par ce critère rigoureux, nous avons constaté que deux erreurs de nœud unique ciblant ou n'étaient pas correctement détectées et qu'un autre nœud unique était prédit comme étant la cible de l'erreur de modèle (Fig. 5a). Cela indique que ces erreurs de modèle ne sont pas observables et que les données de sortie observées peuvent être expliquées par différentes entrées à différents nœuds. A deux exceptions près ((8, 3) et (7, 6)), les erreurs commises par l'algorithme pour les erreurs de modèle par paires simulées concernent ces deux nœuds d'état 1 et 4. Cependant, à l'exception de la combinaison (1, 4), au moins un nœud est correctement prédit.

Détection des nœuds cibles des erreurs de modèle simulées dans le réseau de signalisation UV-B.

(une) Tous les nœuds et toutes les paires de nœuds ont été perturbés par une erreur de modèle simulée. Les nœuds et sont omis, car le signal d'erreur simulé n'a eu aucun effet sur la sortie. Les lignes et les colonnes correspondent aux véritables nœuds cibles de l'erreur de modèle et les nombres dans les cellules sont les nœuds trouvés par le réseau élastique dynamique (NA signifie qu'aucun deuxième nœud n'a été affecté). Les cellules grises indiquent les erreurs commises par le réseau élastique dynamique pour les erreurs de modèle non observables. (b) Un exemple d'erreur de modèle non observable. Les véritables nœuds cibles de l'erreur de modèle sont , mais le réseau élastique dynamique prédit les nœuds cibles . (c) Le remontage du filet élastique dynamique sous contrainte offre une solution alternative. Les deux autres combinaisons et des nœuds ne correspondaient pas aux données de sortie.

Ces résultats démontrent les limites inhérentes à toute tentative de récupérer l'erreur du modèle à partir des sorties observées. Pour une erreur de modèle non observable, la vraie erreur de modèle pourrait correspondre à une valeur légèrement plus grande de la fonctionnelle d'erreur (3c) que le minimum obtenu par le réseau élastique dynamique. Une approche heuristique pour explorer certaines de ces solutions légèrement sous-optimales consiste à réexécuter le réseau élastique dynamique avec certains des nœuds cibles estimés (à partir de la première exécution) exclus et à vérifier si les données de sortie peuvent être adaptées de manière satisfaisante avec le même niveau de parcimonie. Ceci est illustré sur la figure 5b pour la paire de nœuds (9, 1), qui a été prédit être (9, 3) par notre critère. Le remontage de l'élastique dynamique sous la contrainte identifie les nœuds corrects (9, 1), voir Fig. 5c. Les deux autres combinaisons et ne fournissent pas un ajustement satisfaisant aux données (Fig. supplémentaire S9). Pour le réseau de signalisation UVB, nous constatons que les solutions légèrement sous-optimales identifiées par cette heuristique contiennent toujours la configuration de nœud cible correcte. L'explosion combinatoire de cette stratégie ne devrait généralement pas être un problème, grâce à la rareté des prédictions dynamiques du réseau élastique. La décision, laquelle des ensembles de nœuds cibles prédits, ou , est la bonne ne peut en pratique être prise que lorsque des états supplémentaires sont mesurés. Cependant, cet exemple montre comment le réseau élastique dynamique fournit des informations utiles pour sélectionner d'autres états pour l'observation expérimentale 20,28.


Identificabilité

La résolution du problème inverse dépend (a) du modèle mathématique (b) de la signification des données et (c) des erreurs expérimentales. Dans ce qui suit, nous supposons que le modèle est correctement mis à l'échelle de sorte que les valeurs des paramètres et les variables d'état soient du même ordre de grandeur. Sinon, une mise à l'échelle appropriée doit être appliquée au modèle.

Définitions

Un paramètre est globalement identifiable s'il peut être déterminé de manière unique compte tenu du profil d'entrée vous(t) et en supposant des données continues et sans erreur pour les observables du modèle. S'il existe un nombre dénombrable de solutions, le paramètre est identifiable localement, il est non identifiable s'il existe un nombre indénombrable de solutions. Un modèle est structurellement globalement/localement identifiable si tous ses paramètres sont globalement/localement identifiables 2 2 Notons que ces définitions ne sont pas toujours les mêmes. D'autres définitions sont : Un modèle est structurellement identifiable si sa matrice de sensibilité satisfait à deux conditions : chaque colonne a au moins une grande entrée et la matrice a un rang complet [ [4] ]. Un modèle est localement identifiable s'il est globalement identifiable dans un voisinage du paramètre [ [5] ].
.

Pratique ou a postériori l'analyse d'identifiabilité étudie si les paramètres peuvent être déterminés globalement ou localement avec les données expérimentales bruitées disponibles. Dans ce cas, localement signifie au voisinage du paramètre obtenu.

A priori identifiabilité

Il existe plusieurs techniques pour déterminer a priori identifiabilité globale du modèle, mais pour des situations réalistes (c'est-à-dire des modèles non linéaires d'une certaine taille), il est très difficile d'obtenir des résultats. Néanmoins, il est conseillé de toujours effectuer une a priori analyse parce que les méthodes d'estimation des paramètres peuvent avoir des problèmes avec des systèmes localement identifiables ou non identifiables. Les packages d'algèbre symbolique comme maple [ [15] ] et mathematica [ [16] ] peuvent être d'une grande aide.

Pour les modèles linéaires, l'approche de la transformée de Laplace ou de la fonction de transfert peut être appliquée. Pour les modèles non linéaires, la méthode la plus ancienne et la plus simple à comprendre est le développement de Taylor ou en séries entières [ [17] ]. La fonction observable est développée dans une série de Taylor à un moment donné. Les dérivées temporelles sont évaluées en termes de paramètres, résultant en un système d'équations non linéaires pour les paramètres. Si ce système a une solution unique, le modèle est structurellement identifiable. Pour des exemples simples utilisant la transformée de Laplace (modèle linéaire) et les séries de Taylor (cinétique de Michaelis–Menten), nous renvoyons à Godfrey et Fitch [ [18] ]. Une autre méthode classique est l'approche de transformation de similarité [ [19-21] ]. Ces deux méthodes ont été comparées sans préférence décisive [ [22] ]. Récemment, des méthodes ont été développées qui utilisent des techniques d'algèbre différentielle [ [23] ]. En outre, un outil logiciel accessible au public, daisy [ [24] ], est disponible et peut vérifier l'identifiabilité d'un système non linéaire. daisy est implémenté dans le langage symbolique reduce [ [25] ].

A postériori identifiabilité

La difficulté d'estimer les paramètres dans un modèle mathématique quantitatif n'est pas tant de savoir comment les calculer, mais plus d'évaluer la qualité des paramètres obtenus car cela ne dépend pas seulement de la façon dont le modèle décrit le phénomène étudié et de l'existence de un ensemble unique de paramètres, mais aussi de savoir si les données expérimentales sont en nombre suffisant, suffisamment significatives et suffisamment précises. En ce qui concerne les deux premières exigences, une quantité de données suffisante et significative, il est clair que, quelle que soit la méthode utilisée pour ajuster un modèle avec des données expérimentales : estimer m paramètres inconnus, il faut au moins m valeurs expérimentales. D'autre part, il n'est pas nécessaire d'avoir des données expérimentales pour toutes les variables d'état impliquées dans le modèle à tous les moments possibles, souvent seules quelques mesures pour le bon observable à des moments significatifs sont nécessaires. La dernière question, des données suffisamment précises, est liée au fait que les erreurs de mesure impliquent que nous n'avons pas de points de données précis pour ajuster notre modèle, mais que chaque point représente tout un nuage de valeurs de données possibles, impliquant également que les paramètres inférés ne sont pas des valeurs ponctuelles mais sont contenues dans un nuage. Selon le modèle, le nuage de valeurs de paramètres possibles varie en taille et en forme et peut être beaucoup plus grand que l'incertitude d'origine dans les données.

La méthode la plus appliquée [ [12, 26] ] pour étudier cette incertitude dans les paramètres est de calculer la matrice de sensibilité J de l'équation (6) évaluée pour les points de données donnés et le vecteur de paramètre obtenu par l'ajustement des données. Cela peut être fait soit par différenciation finie, soit en résolvant les équations variationnelles 3 3 Les équations variationnelles sont obtenues en prenant la dérivée du système DAE (1) par rapport aux paramètres. Cela se traduit par m Systèmes DAE dans les variables ∂X(t,p)/∂pje, je = 1,…,m.
. Notez qu'il s'agit d'une analyse linéaire et locale à la fois par rapport et par rapport aux points de données donnés.

Une représentation graphique de l'ellipsoïde et de la sphère est donnée sur la figure 1 pour le cas 2D. Supposons que seul le premier k les plus grandes valeurs singulières satisfont (12), alors seulement la première k entrées de z sont estimés avec la précision requise. Si un axe principal de l'ellipsoïde fait un angle significatif avec l'axe dans l'espace des paramètres (c'est-à-dire qu'il existe plus d'une entrée significative dans le vecteur propre), cela correspond à la présence d'une corrélation entre les paramètres dans . Dans ce cas, seule une combinaison de paramètres peut être déterminée.

Exemple d'une région de confiance ellipsoïdale et d'une sphère de précision dans les paramètres du cas 2D p1 et p2 sont corrélées, la combinaison linéaire z1 est bien déterminé, alors que z2 n'est pas. L'intervalle de confiance dépendant, pje, pour un paramètre est donné par l'intersection de l'ellipsoïde avec l'axe des paramètres l'intervalle de confiance indépendant, Δ je pje, par la projection sur l'axe.

Pour résumer, le niveau de bruit dans les données, en combinaison avec l'exigence de précision pour les estimations des paramètres, définit le seuil pour les valeurs singulières significatives dans la matrice . Le nombre de valeurs singulières dépassant ce seuil détermine le nombre de relations paramétriques qui peuvent être dérivées de l'expérience. La relation entre ces relations et les paramètres individuels est décrite par les colonnes correspondantes de la matrice V.

Enfin, Eqn (10) indique qu'avoir, par exemple, des données deux fois plus précises de sorte que l'écart type ?? est réduit de moitié diminuera les rayons le long des axes principaux de l'ellipsoïde d'un facteur 2. Par conséquent, en cas de très petites valeurs singulières ??je (c'est-à-dire des ellipsoïdes fortement allongés), des données plus précises obtenues par l'expérimentateur n'amélioreront pas beaucoup la qualité des estimations des paramètres correspondants. Dans un tel cas, on a certainement besoin de mesures supplémentaires d'un type différent (par exemple, différentes composantes, différents points temporels, ou dans le cas d'équations aux dérivées partielles, différents points spatiaux).

Autres approches

Hengl et al. [ [27] ] proposent une analyse non linéaire : ajustement répété pour différentes estimations initiales du vecteur de paramètre. La matrice vectorielle de paramètres résultante est analysée avec Alternating Conditional Expectation [ [28] ], ce qui entraîne des transformations optimales pour que les paramètres aboutissent à un modèle identifiable. Cette approche est implémentée dans matlab [ [29] ]/PottersWheel [ [30] ].

Enfin, nous souhaitons mentionner une idée intéressante décrite [ [31, 32] ] concernant l'estimation des paramètres basée sur les clusters. This approach uses the sensitivity matrix to define subsets of state variables that depend on a subset of the parameters. The parameters are then split into two classes: global if state variables from more than one cluster depend on it and local otherwise. A hierarchical parameter estimation is performed reducing the dimensionality. On the high level, the global parameters are fitted by optimization of the clusters and, recursively, parameters in each cluster are estimated.

Example insignificant data

On the basis of a very simple artificial example [ [33, 34] ], we show the influence of the experimental data on the parameter determinability.

Suppose the initial concentration of the state variables, [S0], [E0] and [C0] is known, and the concentration of [C] is measured rather precisely at regular time points t = 1,…,20. For this example, the ‘measurements’ are generated artificially by adding an independent, normally distributed perturbance with zero expectation and a fixed variance to the model results (red +-marks in Fig. 2). The initial parameter values are p0 = (6,0.8,1.2). With these parameter values, the model results are given by the solid lines in the left plot in Fig. 2. Fitting the model to these measurements with the Levenberg–Marquardt method (see below) results in the parameter vector (for the model results, see Fig. 2, right).

Model results for initial (left) and final (right) parameter vector, black: [S], red: [C], green: [E] and measurements of [C]: red +.

Confidence region Δk (cf. Eqn 7) in parameter space around computed parameter vector (origin in the plots) and its projection on the parameter-planes. The region contains the true parameter vector with a 95% probability. Left: 20 measurements at t = 1,…,20 right: 20 measurements at time points distributed uniformly over [ [6, 20] ].

Suppose now that it is not possible to measure before time t = 6 but that we take 20 samples of the complex at regular times from t = 6,…,20. Suppose also that the same parameter vector results from minimizing the least squares error e T e. In this case, the confidence region gives much more reason for distrusting the result. As can be seen in Fig. 3 (right), the true parameter vector now lies in a long elongated ‘cigar’ and especially for k1 et k2 we can not even trust the order of magnitude.

Comparing these matrices, corresponding to 20 measurements uniformly distributed in the time interval [ [6, 20] ], with the matrices in Eqn (19), which correspond to 20 measurements at t = 1,2,…,20, it is clear that k3 still can be determined with good accuracy, and even the combination of k1 et k2 can be determined reasonably well, but the third principal axis of the ellipsoidal confidence region has increased almost by a factor of 30! This implies that it is no longer possible to determine k1 et k2 individually.

From the discussion above, it is clear that it is not easy to a priori give an indication whether experimental data are sufficient in number and sufficiently significant. With three ‘lucky’ data points, one can estimate three parameters, but 20 data points in a region where ‘nothing happens’ are not sufficient.

Next, we examine the influence of experimental noise, (i.e. whether the experimental data are sufficiently accurate). Parce que C(??) is proportional to the variance of the measurement error distribution, the principal axes of the ellipsoidal confidence region are proportional to the standard deviation. Roughly speaking: reducing the (standard deviation of the) error by a factor of two, implies that a parameter, or combination of parameters, can be determined more accurately by a factor of two. This means that to shrink the ellipsoidal confidence region for the t > 6 experiment (Fig. 3, right) such that it fits into an ‘accuracy’-sphere that is equal to the experiment with measurements between 1 and 20, one has to reduce the variance of the experimental error beyond reasonable experimental accuracy.

Finally, if we just look at the computable information from the Fisher matrix we get for the confidence intervals:

Exp. ?? (k1) ?? (k2) ?? (k3) ?? je (k1) ?? je (k2) ?? je (k3)
[1,20] 0.033 0.028 0.005 0.076 0.067 0.005
[6,20] 0.074 0.047 0.004 2.217 1.267 0.060

Also, this simple to compute information shows that, for the second case, the parameters are strongly correlated and the model is not identifiable.


Practical applications

Good quantitative models have practical applications. Dieter Oesterhelt (Max Plank Institute of Biochemistry, Martinsried, Germany) demonstrated the down-to-earth potential of models in systems biology using examples from halophilic archaea, including a simple mechanism for regulating the flagellar motor in a way that guarantees an optimal response of the cell to light. He also presented a flux-balance analysis model with about 700 reactions that is capable of predicting growth rates in response to the addition of 16 amino acids to the growth medium. Nicolas Le Novère (European Bioinformatics Institute, Cambridge, UK) presented models of neurological signaling at multiple scales, where lower scales are used to build helpful abstractions at higher scales. Such models are designed to contribute towards understanding dopamine biochemistry and can generate up to 10 terabytes of data for simulating 1 second in the life of a synapse.

At a simpler level, Carolyn Talcott (SRI International, Menlo Park, USA) presented an abstract model of the central pattern generator located in the buccal ganglia of the marine mollusc Aplysie, a long-time model organism for neurobiology. Her simple model considered only a small number of discrete levels of excitation for each neuron but was able to faithfully reproduce the predicted behavior of a more detailed continuous, ordinary differential equation model. Thus, instead of the continuous-state-space approach of differential equations, a logic-based discrete-state-space description of the system is used. This is conceptually much simpler.

The future is likely to see many more mechanistically understood models of actual biological systems that rely on computational methods in systems biology for their construction. The conference showed that there is the potential to greatly enhance our understanding of life if theoretical computer scientists, practical programmers, theoretical biologists and experimental biologists work closely together in order to develop the tools that are needed for constructing and analyzing systems biological models.


Measurement equation of Tree-Step pathway model - Biology

a School of Chemistry, University of Manchester, Sackville St., Manchester, UK
E-mail: [email protected], [email protected]
Fax: +44 (0)1613064556
Tél : +44 (0)1613064492

b The Manchester Interdisciplinary Biocentre, The University of Manchester, 131 Princess St., Manchester, UK

c School of Electrical and Electronic Engineering, University of Manchester, Manchester, UK

d School of Mathematics, University of Manchester, Manchester, UK

Résumé

Mathematical modelling offers a variety of useful techniques to help in understanding the intrinsic behaviour of complex signal transduction networks. From the system engineering point of view, the dynamics of metabolic and signal transduction models can always be described by nonlinear ordinary differential equations (ODEs) following mass balance principles. Based on the state-space formulation, many methods from the area of automatic control can conveniently be applied to the modelling, analysis and design of cell networks. In the present study, dynamic sensitivity analysis is performed on a model of the IκB–NF-κB signal pathway system. Univariate analysis of the Euclidean-form overall sensitivities shows that only 8 out of the 64 parameters in the model have major influence on the nuclear NF-κB oscillations. The sensitivity matrix is then used to address correlation analysis, identifiability assessment and measurement set selection within the framework of least squares estimation and multivariate analysis. It is shown that certain pairs of parameters are exactly or highly correlated to each other in terms of their effects on the measured variables. The experimental design strategy provides guidance on which proteins should best be considered for measurement such that the unknown parameters can be estimated with the best statistical precision. The whole analysis scheme we describe provides efficient parameter estimation techniques for complex cell networks.


INTRODUCTION

Systems biology may be broadly defined as the integration of diverse data into useful biological models that allow scientists to easily observe complex cellular behaviors and to predict the outcomes of metabolic and genetic perturbations. In this report we present kMech, a Cellerator ( Shapiro et al., 2003) language extension that describes a suite of enzyme mechanisms suitable for the mathematical modeling of enzyme-related pathways. Cellerator is a tool for generating reaction network models of cellular processes. It contains reaction mechanisms based on mass action kinetics for elementary catalytic and non-catalytic reactions, along with other regulatory relationships ( Shapiro et al., 2003). kMech supplements Cellerator with a suite of catalytic reaction mechanisms involving multiple substrates and products and regulatory mechanisms. In kMech, each mechanism has been codified to generate a set of elementary reactions that can be translated by Cellerator into ordinary differential equations (ODEs) solvable by Mathematica™, a widely used commercial computer algebra system that integrates numeric and symbolic computational engines with a graphical output and a programming language. Alternatively, Cellerator can generate ODEs in systems biology markup language (SBML) for other simulators ( Hucka et al., 2003).

Traditional enzyme modeling approaches, for example, GEPASI ( Mendes, 1997) or JARNAC ( Sauro et al., 2003), use the Michaelis–Menten kinetic equation for one substrate/ one product reactions and the King–Altman method to derive equations for more complex multiple reactant reactions. These types of equations are called steady-state velocity equations since the derivatives of the concentration of each reactant in the model over time are set to zero to simplify a set of non-linear differential equations to linear algebric equations ( Falk et al., 1998). Therefore, the kinetic model based on this approach has embedded in it the steady-state hypothesis. In contrast, the model generated by kMech/Cellerator consists of non-simplified, non-linear, differential equations that describe the rates of change of each reactant in the model over time. For the solution of these differential equations, users can input the rate constants (kF, kr) if they have been experimentally measured. Alternatively, these rate constants can be approximated from easily measured kinetic constants (Km, kchat) using the Lambda and Omega approximation methods presented in this paper. Moreover, if the assumption of equilibrium between substrate and enzyme is not favored, values of Lambda and Omega can be adjusted to fit the experimental data.

kMech is built in a modular manner such that complex enzyme mechanisms may be constructed from simpler ones. This design is sufficiently versatile to accommodate enzyme mechanisms such as single and multiple substrate (e.g. Ping Pong and Bi Bi) enzyme kinetic reactions, and feedback inhibition (e.g. allosteric, competitive and non-competitive) mechanisms.


The IUPS Physiome Project: a framework for computational physiology

The IUPS Physiome Project 1 is an internationally collaborative open-source project to provide a public domain framework for computational physiology, including the development of modelling standards, computational tools and web-accessible databases of models of structure and function at all spatial scales. A number of papers in this volume deal with the development of specific mathematical models of physiological processes. This paper stands back from the detail of individual models and reviews the current state of the IUPS Physiome Project including organ and organ system continuum models, the interpretation of constitutive law parameters in terms of micro-structural models, and markup languages for standardizing cellular processes. Some current practical applications of the physiome models are given and some of the challenges for the next 5 years of the Physiome Project at the level of organs, cells and proteins are proposed.


Information additionnelle

Intérêts concurrents

Les auteurs déclarent une absence d'intérêts financiers en compétition.

Authors’ contributions

WL, JSL, and XSL designed the statistical model. WL and JK developed the algorithm, designed and performed the analysis. XH and CHC performed sgRNA efficiency prediction analysis. WL, JK, and XSL wrote the manuscript with help from all other authors. XSL and MB supervised the whole project. All authors read and approved the final manuscript.


Voir la vidéo: Modèle de Drude de la conduction électrique (Janvier 2023).