Brevet FR3160045A1 — publié le 20250912

L’invention porte sur un procédé de calibration temporelle d’un arbre phylogénétique par traitement numérique, et comprend : -la réception sous format numérique d’un ensemble d’évènements datés d’isolements ou de transmission d’un agent pathogène et d’un ensemble d’évènements à dater d’isolements ou de transmission de l’agent pathogène, organisés sous forme d’arbre phylogénétique selon un ordre de traversée de la racine vers les feuilles. Plusieurs évènements ont une date à déterminer ; -formation d’une matrice Hessienne d’une fonction F d’objectif logarithmique d’un modèle d’horloge moléculaire de l’agent pathogène, définissant la relation entre les distances génétiques, les dates et les durées entre les évènements, la matrice Hessienne H comprenant un bloc T dit local définissant les branches entre lesdits évènements de l’arbre phylogénétique. Le bloc T a une taille proportionnelle au nombre d’évènements de l’arbre phylogénétique dont la date est à déterminer.; -identifier les dates à déterminer par une méthode d’optimisation de Newton. Figure à publier avec l’abrégé : Fig. 1 Procédé de calibration temporelle d’un arbre phylogénétique L’invention concerne le traitement des données épidémiques, et en particulier le traitement de séquences génétiques pour extrapoler des informations concernant les transmissions épidémiques. Une épidémie infectieuse implique des événements de transmission de pathogènes. Les dates de ces événements de transmission sont un outil précieux pour aider la prise de décision pour le contrôle épidémique. L’utilisation de ces dates permet notamment de déterminer l’entrée d’un pathogène dans une unité géographique (service hospitalier, ville, pays) ou d’identifier les pics de fréquence des événements de transmission. Les dates des événements de transmission permettent également d’estimer rétrospectivement le taux de reproduction de l’épidémie appelé R0. Les dates des évènements de transmission ne sont pas connues exactement. Ces dates peuvent être estimées grâce à l’analyse phylodynamique d’échantillons de pathogènes. Cette analyse phylodynamique comporte la récupération d’échantillons à partir d’hôtes infectés ou colonisés, l’obtention de caractéristiques héritables des pathogènes, typiquement des séquences ADN ou ARN, la reconstruction de la phylogénie (généalogie) des pathogènes à partir de ces séquences, et l’estimation des dates de transmission à partir de la phylogénie. La phylodynamique repose sur la création d’un arbre phylogénétique dont les feuilles sont des échantillonnages de pathogènes et dont les nœuds coïncident avec des évènements de transmission du pathogène entre des hôtes lorsque ces nœuds relient des feuilles issues d’hôtes distincts. Chaque feuille représente un pathogène directement observé et chaque nœud représente un évènement de divergence, non directement observé. Les épidémies pouvant porter sur un grand nombre de cas, la détermination des dates des évènements doit être réalisable par des méthodes numériques impliquant potentiellement des ressources de calcul importantes. Les méthodes de datation connues sont basées sur des modèles d’évolution des séquences ARN ou ADN, en exploitant les dates connues des feuilles de l’arbre phylogénétique. On considère généralement que la date d’échantillonnage correspond à chaque feuille. La datation des nœuds de l’arbre est désignée par le terme de calibration temporelle. La calibration temporelle est une étape cruciale de toute analyse phylodynamique. Ainsi, sa précision et sa rapidité d’exécution sont des paramètres prépondérants d’un système de surveillance épidémique basé sur l’analyse de séquences génétiques du pathogène. Le but est donc de fixer le plus précisément possible des évènements de transmission, par exemple en quantifiant les durées entre deux nœuds parents. De nombreuses méthodes de calibration temporelle ont été développées. Pour aboutir à des résultats précis et robustes, les méthodes de calibration temporelle connues présentent un problème commun de temps de calcul qui augmente plus vite que le nombre d’échantillons de l’arbre phylogénétique à analyser (typiquement exponentiellement). Pour une épidémie importante, il n’est actuellement pas possible de réaliser une calibration temporelle en un temps raisonnable. A contrario, en un temps réduit, il n’est pas actuellement possible d’obtenir une précision et une robustesse suffisantes pour le traitement des données. L’invention vise à résoudre un ou plusieurs de ces inconvénients. L’invention porte ainsi sur un procédé de calibration temporelle d'un arbre phylogénétique par traitement numérique, comprenant : -la réception sous format numérique d'un ensemble d'évènements datés d'isolements ou de transmission d'un agent pathogène et d'un ensemble d'évènements à dater d'isolements ou de transmission de l'agent pathogène, les évènements étant organisés sous forme d'arbre phylogénétique comprenant une racine, des feuilles, des nœuds ordonnés selon un ordre de traversée de la racine vers les feuilles dans l'arbre phylogénétique, la racine, les feuilles et les nœuds étant chacun associé à un desdits évènements, la racine, les feuilles et les nœuds étant reliés par des branches ayant une longueur correspondant à une distance génétique entre des évènements, plusieurs desdits évènements ayant une date à déterminer, chacune desdites branches étant associée à un évènement ascendant et à un évènement descendant de l'arbre phylogénétique selon ledit ordre; -formation d'une matrice Hessienne d'une fonction F d'objectif logarithmique d'un modèle d'horloge moléculaire de l'agent pathogène, définissant la relation entre les distances génétiques, les dates et les durées entre les évènements, la matrice Hessienne H comprenant un bloc T dit local définissant les branches entre lesdits évènements de l'arbre phylogénétique, le bloc T ayant une taille proportionnelle au nombre d'évènements de l'arbre phylogénétique dont la date est à déterminer; -identifier les dates à déterminer dans un vecteur solution w, par une méthode d'optimisation de Newton de résolution de l'équation H(F(v)) s = g(F(v)), avec v un vecteur courant, avec g une fonction gradient, et s une inconnue ou direction de Newton, en fixant le vecteur solution w à la valeur du vecteur courant v lorsque la variation du vecteur courant v entre deux itérations de la méthode de Newton est inférieure à un seuil ou lorsque la variation de la fonction d'objectif F entre deux itérations de la méthode de Newton est inférieure à un seuil. L’invention porte également sur les variantes suivantes. L’homme du métier comprendra que chacune des caractéristiques des variantes suivantes peut être combinée indépendamment aux caractéristiques ci-dessus, sans pour autant constituer une généralisation intermédiaire. Selon une variante, la matrice Hessienne H est de la forme : avec P un bloc dit global correspondant à des paramètres à déterminer du modèle d'horloge moléculaire de l'agent pathogène ; avec C et C' des blocs correspondants aux dérivées partielles impliquant chacune un évènement et un paramètre global du modèle d'horloge moléculaire de l'agent pathogène. Selon encore une variante, au moins un paramètre global du modèle d’horloge moléculaire de l’agent pathogène est à déterminer. Selon encore une variante, la résolution du système linéaire dans l'optimisation de Newton comprend une décomposition LU creuse ou une décomposition QR creuse mise en œuvre sur le bloc T. Selon une variante, les évènements associés aux feuilles sont datés et les évènements associés aux nœuds sont à dater. Selon une autre variante, la longueur de chaque branche est définie en unités de distance génétique ou en nombre de substitutions génétiques. Selon encore une variante, le modèle d'horloge moléculaire de l'agent pathogène est un modèle d'horloge relâchée à variance additive, comportant des paramètres définissant un taux d'évolution en unités de distance génétique par unité de temps et une dispersion du taux d'évolution représentative de la variabilité de la vitesse d'évolution entre différentes branches de l'arbre phylogénétique. Selon encore une variante, ladite fonction F est de la forme α * log β – log Γ(α) +( α-1) * log d -β* d avec Γ une distribution Gamma ayant les paramètres α=µ *l / (1+ ω) qui correspond à la forme de la distribution et β= 1/ (1+ ω) qui correspond au taux de la distribution, avec ω un paramètre de dispersion du modèle d’horloge relâchée à variance additive. Selon encore une autre variante, dans ledit bloc T, un indice de ligne i correspondant à un unique nœud et la même valeur d'indice de colonne j correspondant à cet unique nœud, les indices de ligne i et de colonne j étant attribués selon l'ordre inverse de la traversée de l'arbre phylogénétique. Selon une variante, chaque nœud est associé à un unique évènement ascendant et au plus à deux évènements descendants. Selon encore une variante, l'amplitude de déplacement dans la direction du vecteur s lors d'une itération de la résolution de Newton est définie pour respecter la règle qu'un évènement descendant est toujours postérieur à son évènement ascendant dans l'arbre phylogénétique. Selon une autre variante, le procédé comprend un calcul d'un niveau d'incertitude de ladite solution par une méthode de Wald, incluant un calcul de la diagonale de l'inverse de la matrice hessienne H de la fonction de vraisemblance appliquée audit vecteur w. Selon encore une autre variante, on utilise une méthode du complément de Schurr pour obtenir la diagonale de l'inverse de la matrice hessienne H lors du calcul du niveau d'incertitude de la solution. L’invention porte également sur un système numérique, comprenant des moyens de traitement configurés pour mettre en œuvre un procédé de calibration temporelle tel que défini ci-dessus. D'autres caractéristiques et avantages de l'invention ressortiront clairement de la description qui en est faite ci-après, à titre indicatif et nullement limitatif, en référence aux dessins annexés, dans lesquels : FIG. 1 est une illustration schématique d’un exemple d’arbre phylogénétique ; FIG. 2 illustre la structure d’un bloc d’une matrice Hessienne obtenue pour l’arbre phylogénétique de la FIG. 1 ; FIG. 3 illustre un exemple d’infrastructure de mise en œuvre d’une calibration selon l’invention ; FIG. 4 illustre le taux d’évolution estimé pour un procédé selon l’invention et une référence de l’état de la technique, avec un premier exemple d’échantillon ; FIG. 5 illustre le pourcentage d’erreur tMRCA pour le procédé selon l’invention et la référence de l’état de la technique, avec le premier exemple d’échantillon ; FIG. 6 illustre comparativement le pourcentage d’erreur sur les dates pour la référence de l’état de la technique et pour l’invention, avec le premier exemple d’échantillon ; FIG. 7 illustre comparativement le pourcentage d’erreur sur les dates pour le procédé selon l’invention et la référence de l’état de la technique, avec le premier exemple d’échantillon ; FIG. 8 illustre le taux d’évolution estimé pour un procédé selon l’invention et une référence de l’état de la technique, avec un deuxième exemple d’échantillon ; FIG. 9 illustre le pourcentage d’erreur tMRCA pour le procédé selon l’invention et la référence de l’état de la technique, avec le deuxième exemple d’échantillon ; FIG. 10 illustre comparativement le pourcentage d’erreur sur les dates pour le procédé selon l’invention et la référence de l’état de la technique, avec le deuxième exemple d’échantillon ; FIG. 11 illustre comparativement le pourcentage d’erreur sur les dates pour le procédé selon l’invention et la référence de l’état de la technique, avec le deuxième exemple d’échantillon. FIG. 12 illustre le temps d’exécution comparatif entre le procédé de l’invention et le procédé de référence de l’état de la technique, en fonction de la taille de l’échantillon, pour un système de traitement numérique donné ; FIG. 13 illustre l’évolution du temps d’exécution du procédé de l’invention pour des échantillons de grande taille. L’invention porte sur un procédé de calibration temporelle d’un arbre phylogénétique par traitement numérique. Le procédé comprend la réception de données numériques d’entrée. Les données numériques d’entrée correspondent à un arbre phylogénétique à calibrer. Les données numériques d’entrée comprennent ainsi : -un ensemble d’évènements datés d’isolements ou de transmission d’un agent pathogène. Les caractéristiques héritables des pathogènes isolés peuvent être obtenus par l’intermédiaire de séquençages d’ADN ou d’ARN de l’agent pathogène, ou par d’autres sources de données telles que la spectroscopie Raman, l’étude de profils de méthylation. Les analyses permettant d’obtenir les caractéristiques héritables des pathogènes isolés peuvent être réalisées en décalé par rapport aux isolements eux-mêmes. Et -un ensemble d’évènements à dater d’isolements ou de transmission de l’agent pathogène. Les données numériques d’entrée des évènements sont organisées sous forme d’arbre phylogénétique comprenant : une racine, des feuilles, des nœuds ordonnés selon un ordre de traversée de la racine vers les feuilles dans l’arbre phylogénétique. Le but de la calibration temporelle de l’arbre phylogénétique est de déterminer des dates d’isolements ou de transmission pour des évènements, ces dates étant initialement indéterminées. La racine, les feuilles et les nœuds sont chacun associés à un des évènements. La racine, les feuilles et les nœuds sont reliés par des liaisons ou branches. Les branches ont une longueur correspondant à une distance génétique entre les évènements qu’elles relient. Chacune des branches est associée à un évènement ascendant et à un évènement descendant de l’arbre phylogénétique. La durée des branches peut se déduire des dates des évènements auxquels ces branches sont associées. Dans un souci de simplification, les dates des évènements datés correspondent aux feuilles, c’est-à-dire des évènements d’isolements, les dates des évènements à dater correspondent à des nœuds de transmission. On peut également envisager que des feuilles aient des dates inconnues ou au contraire que des nœuds soient associés à des dates d’évènements connus. Les branches de l’arbre phylogénétique peuvent être représentées par un ensemble de triplets (identification d’un évènement ascendant, identification d’un évènement descendant, longueur de la branche), chaque triplet représentant une branche entre un évènement ascendant et un évènement descendant. La longueur connue de cette branche est par exemple définie en unités de distance génétique, par exemple en nombre de substitutions ou de mutations. Les feuilles de l'arbre représentent typiquement les séquences d’ADN ou d’ARN observées. Une feuille n'a pas de descendant et n'apparaît donc pas comme un ascendant dans la liste des triplets. L'arbre comporte également une racine (nœud le plus ancien). La racine n'a pas d'ascendant et n'apparaît donc pas comme un descendant dans la liste des triplets. Lors d’une étape antérieure de création de l’arbre phylogénétique, connue en soi de l’art antérieur, la racine est typiquement déterminée par une méthode de groupe extérieur ou par une analyse systématique de toutes les racines possibles pour sélectionner la racine qui maximise la corrélation entre, d'une part, la distance génétique entre la racine et les feuilles et, d'autre part, la date des feuilles. Cette méthode est généralement désignée par une régression de la racine vers les feuilles. La relation entre les distances génétiques et les durées entre les évènements ou intervalles temporels (aussi appelées durées d’évolution) entre ascendants et descendants est décrite par un modèle d’horloge moléculaire. De nombreux modèles d’horloges moléculaires sont connus en soi. Un modèle d’horloge moléculaire connu et robuste est le modèle d’horloge relâchée à variance additive. Un tel modèle d’horloge moléculaire sera avantageusement mis en œuvre dans le cadre de l’invention. Un modèle d’horloge moléculaire comprend typiquement un ou plusieurs paramètres qui décrivent la relation entre les distances génétiques et les intervalles temporels à l’échelle de l’arbre phylogénétique dans son ensemble. Ces paramètres sont désignés par le terme paramètres globaux. Les paramètres locaux du modèle d’horloge moléculaire sont les dates des évènements (en particulier les nœuds) à estimer. Chaque date est locale pour un évènement ou nœud particulier. Dans le modèle d’horloge relâchée à variance additive, deux paramètres globaux sont définis : le taux d’évolution et la dispersion du taux d’évolution. Le taux d’évolution est exprimé en unités de distance génétique par unité de temps et décrit la vitesse globale d’évolution dans l’arbre phylogénétique. La dispersion du taux d’évolution est une grandeur sans signification biologique et décrit la variabilité de la vitesse d’évolution entre les branches de l’arbre. Un paramètre de dispersion égal à zéro définit un modèle d’horloge dit strict. Dans lequel la vitesse d’évolution est constante dans l’ensemble de l’arbre phylogénétique. Le procédé de calibration temporelle vise à estimer les valeurs les plus probables de l'ensemble des dates inconnues et, conjointement (le cas échéant), des paramètres globaux inconnus du modèle d’horloge moléculaire. Les valeurs les plus probables des dates des évènements sont obtenues en maximisant la valeur d’une fonction d’objectif logarithmique associée au modèle d’horloge moléculaire de l’agent pathogène dépendant des dates inconnues des évènements. Pour la calibration de l’arbre phylogénétique, on forme une matrice Hessienne de cette fonction d’objectif logarithmique. Pour rappel, une matrice Hessienne H d’une fonction F à valeurs réelles, appliquées à des paramètres (x1,…,xn), se définit comme : H ij (f)=δ²f/δx i δx j pour i et j compris entre 1 et n. Cette matrice Hessienne comprend un bloc T dit local. Ce bloc T comprend les secondes dérivées partielles de la fonction d’objectif par-rapport aux dates des événements à déterminer. Le bloc T a ainsi une taille proportionnelle au nombre d’évènements de l’arbre phylogénétique dont la date est à déterminer. Dans le contexte de l’invention, le bloc T a la propriété d’être creux, c’est-à dire qu’il comporte des valeurs forcément nulles. En effet, la seconde dérivée partielle de la fonction d’objectif par rapport aux dates de deux nœuds est nulle, dès lors que ces nœuds sont séparés dans l’arbre par au moins un autre nœud. En pratique, chaque ligne ou colonne du bloc T comporte au plus 4 valeurs non nulles car un nœud a au maximum trois voisins connectés dans un arbre binaire. Le nombre de valeurs non nulles dans le bloc T croît donc linéairement avec le nombre de dates à déterminer. Il existe par ailleurs une permutation du bloc T telle que les indices de ses lignes et colonnes (correspondant chacun à un nœud) respectent l’ordre de traversée de l’arbre depuis les feuilles vers la racine. Une telle permutation permet de réduire encore le temps de calcul pour mettre en œuvre l’optimisation de Newton détaillée par la suite. Les méthodes d’optimisation de Newton impliquent l’estimation et l’inversion de la matrice Hessienne. L’inventeur a par ailleurs constaté que les méthodes de résolution de Newton présentaient en théorie un temps de calcul cubique avec la dimension de la matrice pour des matrices usuelles mais présentait un temps de calcul linéaire avec la dimension de cette matrice si cette matrice était creuse. L’inventeur ayant constaté que la matrice Hessienne était creuse dans son application à un arbre phylogénétique pour un modèle d’horloge relâchée à variance additive, il apparaît que la durée de calibration temporelle est ici linéaire avec le nombre de dates à déterminer. Ainsi, les méthodes d’optimisation de Newton, qui semblent de premier abord incompatibles avec la recherche d’un temps de calcul réduit, se sont avérées particulièrement appropriées pour la calibration temporelle. L’invention permet ainsi de réaliser une calibration temporelle d’un arbre phylogénétique en un temps qui évolue linéairement avec le nombre d’échantillons traités. Ainsi, la réalisation d’une calibration temporelle devient possible en un temps raisonnable avec une grande précision, même pour une épidémie d’ampleur. Le procédé de l’invention est probabiliste et repose sur une estimation au maximum de vraisemblance des paramètres du modèle d’évolution. L’optimisation du procédé selon l’invention est réalisée sans approximation puisqu’elle est basée sur l’expression analytique de l’ensemble des dérivées partielles du modèle d’horloge moléculaire relâchée à variance additive. Les figures 1 et 2 illustrent respectivement un exemple d’arbre phylogénétique simplifié et une permutation feuilles vers racine correspondante d’un exemple de matrice Hessienne pour cet arbre phylogénétique. L’arbre phylogénétique comprend ainsi des nœuds 111 à 119 et des feuilles 201 à 210. Dans le bloc T illustré à la FIG. 2 , les valeurs possiblement non nulles sont représentées par une croix et les valeurs toujours nulles par un zéro. Les numérotations des lignes et des colonnes du bloc T reflètent la position des nœuds intérieurs de l’arbre phylogénétique. L’optimisation de Newton-Raphson implique de rechercher de façon itérative la valeur d’un vecteur de paramètres qui maximise une fonction d’objectif ou, plus formellement, qui annule la dérivée de cette fonction d’objectif. A chaque itération, le vecteur de paramètres est déplacé dans une direction donnée par le produit de l’inverse de la matrice hessienne et du gradient de la fonction d’objectif. Ce déplacement correspond en pratique à la solution exacte du maximum sous l’hypothèse simplificatrice que la fonction d’objectif est une forme quadratique. Lorsque cette hypothèse est vraie, la méthode de Newton-Raphson maximise directement l’objectif en une étape unique. Dans le cas général où cette hypothèse simplificatrice est une approximation (par exemple avec une calibration temporelle), un nombre fini d’itérations s’avère nécessaire pour atteindre l’objectif. A chaque itération, la valeur de l’objectif est comparée à la valeur précédente et l’algorithme stoppe lorsque la valeur d’objectif est virtuellement égale à la précédente, ce qui signale que le procédé a atteint un plateau. La fonction d’objectif est ici la log-vraisemblance du modèle d’horloge moléculaire et le vecteur de paramètres est composé des dates à estimer et des paramètres globaux (le cas échéant). Lorsque des paramètres doivent obéir à des contraintes (positivité, inégalité par exemple), ils sont typiquement remplacés par un paramètre projeté sur les nombres réels (de moins l’infini à plus l’infini) à l’aide d’une fonction de liaison comme la fonction logarithme, qui projette une valeur non-négative sur les nombres réels. Les dates sont ici exprimées en ancienneté par-rapport au présent, le présent étant typiquement défini comme la date d’échantillonnage de la dernière feuille. Les dates sont non-négatives, la date d’un ascendant ne pouvant être postérieure à la date de son descendant le plus ancien (ou feuille). De même, les paramètres globaux du modèle d’horloge moléculaire relâchée à variance additive sont tous deux non-négatifs, une vitesse d’évolution ou une dispersion négative n’ayant pas de sens biologique. Dans ce contexte, tous les paramètres à estimer sont projetés sur les nombres réels par une fonction de liaison logarithmique. Les notations suivantes peuvent être utilisées : -d i est la distance génétique d’une branche se terminant par un évènement d’indice i ; -l i est la durée de la branche se terminant par l’évènement d’indice i ; -t i est l’âge d’une feuille correspondant à l’évènement d’indice i ; en notant p(i) le parent de i et c(i) l’enfant de i, on a donc l i = t p(i) -t i ou l=t p -t pour simplifier la notation ; -le taux d’évolution µ du modèle d’horloge relâchée à variance additive et un paramètre ω un paramètre de dispersion (égal à 1 pour un modèle strict dans la mise en œuvre qui suit, ou égal à 0 pour un modèle strict dans la publication Didelot) selon le principe enseigné par la publication Didelot (Didelot, X., Siveroni, I. & Volz, E. M. ‘Additive Uncorrelated Relaxed Clock Models for the Dating of Genomic Epidemiology Phylogenies. Mol. Biol. Evol. 38, 307–317 (2021)’). La distance génétique d, interprétée comme le nombre de modifications dans une branche dans le modèle d’horloge relâchée à variance additive, suit une distribution Gamma Γ avec les paramètres α=µ *l / (1+ ω) qui correspond à la forme de la distribution ; β= 1/ (1+ ω) qui correspond au taux de la distribution. Cette distribution a une moyenne de µ *l, une variance de µ *l / (1+ ω) et un coefficient de variation de √((1+ ω)/ µ *l La densité correspondante (simplifiée en retirant l’indice i) peut être définie par : P(d|l, µ , ω)= (β α /Γ(α)) * d α-1 * e - βd = d ( µ *l / (1+ ω) )-1 * e -d / (1+ ω) * (1+ ω) -( µ *l /(1+ ω)) / Γ(µ *l /(1+ ω)) On peut alors définir une fonction d’objectif ou de probabilité logarithmique de branche, associée au modèle d’horloge moléculaire : L = α * log β – log Γ(α) +( α-1) * log d -β* d =α * log (β*d) – log Γ(α) - log d -β* d = (µ *l / (1+ ω)) * log (d/ (1+ ω))-log Γ(µ *l / (1+ ω))- log d -d/(1+ ω) On peut alors définir une fonction d’objectif ou de probabilité logarithmique de l’arbre par : Par un procédé d’optimisation de Newton-Raphson, on maximise conjointement les âges des nœuds t 1 , …, t n-1 ainsi que les paramètres du modèle d’horloge. Le procédé d’optimisation doit garantir la contrainte d’ordre de l’arbre, à savoir t p(i) > t i >t c ( i) pour l’ensemble des évènements i. Tous les paramètres incluant des dates de nœuds sont strictement positifs et leur erreur est proportionnelle, on optimise donc leur logarithme. L’optimisation de Newton-Raphson s’avère avantageuse sur des arbres phylogénétiques de grande taille, du fait de sa convergence rapide en nombre d’itérations, même si chaque itération est en théorie plus consommatrice de ressources de calcul. L’optimisation de Newton-Raphson est contrainte. L’algorithme d’optimisation doit garantir le respect des contraintes durant les itérations, la fonction d’optimisation logarithmique étant indéfinie hors de l’intervalle de faisabilité. On désigne par x s le vecteur de paramètres à l’étape s (incluant les paramètres de modèle, ici après les dates des n-1 nœuds). Soit Δx s = ∇Λ(x s )=d Λ/dx s le gradient de la fonction d’objectif logarithmique et Δ²x s = ∇² Λ (x s )=d² Λ /dx s ² la Hessienne de la fonction d’objectif logarithmique pour x s . Chaque nœud a une date t i = e xi et la contrainte d’ordre de l’arbre phylogénétique impose que t p ( i ) >t i . On peut définir l’amplitude d’une itération γ s et la règle d’itération x s+1 = x s − γ s ( Δ²x s ) −1 Δx s , avec ( Δ²x s ) −1 Δx s la direction de Newton. L’itération est faisable (dont le résultat a une fonction d’objectif finie) si et seulement si, pour toutes les valeurs de i, t p(i),s+1 >t i,s+1 , c’est-à dire que e x p(i),s+1 > e x i,s+1 en notant que la fonction exponentielle augmente de façon monotone, Pour tout i>1 : 0 p(i),s+1 - x i,s+1 = x p(i),s+1 + γ s ( Δx p(i),s )- x i,s - γ s Δx i,s Ainsi, une amplitude optimale d’itération γ s respecte la règle qui suit : 0 s * p(i),s - x i,s )/( Δx i,s - Δx p(i),s )= γ s + pour i>1. Dans cet intervalle, l’amplitude optimale d’itération est celle qui maximise la fonction d’objectif. On utilise avantageusement la méthode de Brent pour déterminer l’amplitude optimale qui maximise la fonction d’objectif dans cet intervalle. Lorsque la distance entre deux valeurs successives du vecteur de paramètres x s est inférieure à un certain seuil (par exemple 1/1 000 000), on peut stopper les itérations de l’optimisation de Newton. L’optimisation de Newton suppose l’inversion de la matrice Hessienne. L’inversion d’une Hessienne peut s’avérer instable. Afin de contourner cette problématique, on peut mettre en œuvre une régularisation de Tikhonov pour l’étape de Newton, avec p= (Ф+1) * (Δ²x+Ф* tr(Δ²x). I ) -1 Δx I étant une matrice identité. P tend vers l’itération de Newton avant régularisation lorsque Ф tend vers 0 et vers le gradient de cette itération lorsque Ф tend vers l’infini. En cas d’arrêt prématuré de l’exécution de l’optimisation de Newton par le fait d’une instabilité numérique, on augmente la valeur de Ф. La matrice Hessienne peut contribuer à accélérer la convergence de l’optimisation et fournit une information d’approximation concernant l’indice de fiabilité. L’optimisation de Newton mise en œuvre a une complexité linéaire avec la variable n du fait du caractère creux de la matrice Hessienne, ce qui permet d’obtenir une convergence rapide pour un coût de traitement assez réduit à chaque itération. On peut partitionner le gradient des paramètres en un vecteur t des logarithmes des gradients de date des nœuds, et en un vecteur du gradient des paramètres liés au modèle, soit Δx=[t,θ] T Ce qui peut se traduire par la forme suivante de la Hessienne : Le bloc P apparait lorsque des paramètres du modèle d’horloge moléculaire sont inconnus préalablement à la calibration temporelle. Le bloc P est dit global car il implique des paramètres communs à tous les nœuds de l’arbre. La taille du bloc P est définie par le modèle d’horloge moléculaire choisi, et est indépendante de la taille de l’arbre. Les blocs C et C’ incluent les dérivées partielles impliquant un nœud et un paramètre global. Plutôt que de réaliser l’inversion complète de la matrice Hessienne, on peut aussi réaliser une inversion par blocs de matrice dite méthode de complément de Schur, par décomposition LU creuse, lorsque la matrice creuse H respecte la permutation racine vers les feuilles. Cette inversion par bloc peut s’exprimer comme suit dans le cas plus générique illustré ci-dessus : La matrice A est une matrice de k par k dont l’inversion est peu gourmande en temps de calcul. La direction de Newton-Raphson peut être exprimée comme suit : Par ailleurs, la vitesse de convergence de la calibration temporelle est quadratique : chaque itération réduit l’erreur au carré. Une fois les valeurs optimales (et le cas échéant des paramètres globaux) ont été obtenues, une estimation de leur incertitude (par exemple sous la forme d’un intervalle de confiance à 95%) peut être obtenue par la méthode de Wald. Cette méthode repose sur le caractère asymptotiquement gaussien de la distribution des paramètres estimés au maximum de vraisemblance, la variance de la distribution gaussienne de chaque paramètre étant l’entrée diagonale de l’inverse du négatif de la matrice hessienne. Comme pour la solution de la direction de Newton, la méthode de Wald repose sur l’inversion de la matrice Hessienne. Une telle inversion est en général infaisable en pratique pour les problèmes de grande taille. Cependant, cette matrice est typiquement définie négative au maximum de vraisemblance. La combinaison du caractère défini-négatif, du caractère creux et de la structure en arbre permet d’obtenir les variances de Wald, soit la diagonale de l’inverse du négatif de la matrice Hessienne, de façon très efficace en utilisant la décomposition de Cholesky creuse. En effet, l’inverse d’une matrice creuse est généralement dense. Mais avec une permutation de la racine vers les feuilles, on peut bénéficier du fait que l’inverse de la décomposition de Cholesky du bloc T est creuse. Les entrées diagonales (et uniquement celles-ci) de l’inverse de la matrice Hessienne sont obtenues sans perdre de temps sur les entrées non-diagonales, en multipliant uniquement les vecteurs nécessaires de l’inverse de la décomposition de Cholesky. L’ensemble de cette procédure permet d’obtenir une estimation d’incertitude pour tous les paramètres (dates et paramètres globaux) avec une complexité computationnelle linéaire, là où une implémentation standard de complexité cubique serait infaisable. On peut mettre en œuvre la méthode du complément de Schurr pour obtenir la diagonale de l’inverse de la matrice Hessienne H lors du calcul du niveau d’incertitude de la solution. On peut comparer les performances d’un procédé de calibration temporelle aux procédés TREETIME (approche fréquentiste, approximative et basée sur les séquences, décrite dans Sagulenko, P., Puller, V. & Neher, R. A. TreeTime: Maximum-likelihood phylodynamic analysis. Virus Evol. 4, (2018)), et TREEDATER (approche fréquentiste, approximative et basée sur les distances, décrite dans Volz, E. M. & Frost, S. D. W. Scalable relaxed clock phylogenetic dating. Virus Evol. 3, (2017). Comme illustré sur les diagrammes des figures 4 à 7, on a illustré l’erreur d’estimation d’un procédé selon l’invention avec le procédé dit Treedater, supposé être le procédé connu le plus précis. Les simulations ont été réalisées sur 1024 phylogénies simulées de 50 isolats présentant un signal temporel de bonne qualité (taux d’évolution μ=10, dispersion ω=1, échantillonnage des 1% d’isolats les plus récents). L’erreur d’estimation du procédé selon l’invention est globalement plus faible que celle de Treedater, pour tous les paramètres dont le taux d’évolution μ, l’ancêtre commun le plus récent (tMRCA) et les dates des nœuds. L’erreur sur les dates monte à plus de 60% pour Treedater et ne dépasse pas 30% pour le procédé de l’invention. La FIG. 4 illustre le taux d’évolution estimé (à gauche l’invention, à droite Treedater), la FIG. 5 illustre le pourcentage d’erreur tMRCA (à gauche l’invention, à droite Treedater), la FIG. 6 illustre comparativement le pourcentage d’erreur sur les dates pour treedater (en ordonnée) et pour l’invention (en abscisse), et la FIG. 7 illustre selon une autre présentation le pourcentage d’erreur sur les dates pour le procédé de l’invention et pour treedater (à gauche l’invention, à droite Treedater). Comme illustré sur les diagrammes des figures 8 à 11, on a illustré l’erreur d’estimation d’un procédé selon l’invention avec le procédé dit Treedater. Les simulations ont ici été réalisées sur 1024 phylogénies simulées de 50 isolats présentant un signal temporel de bonne médiocre (taux d’évolution μ=10, dispersion ω=10, échantillonnage des 1% d’isolats les plus récents). L’erreur d’estimation du procédé selon l’invention est globalement plus faible que celle de Treedater, pour tous les paramètres dont le taux d’évolution μ, l’ancêtre commun le plus récent (tMRCA) et les dates des nœuds. L’erreur sur les dates maximale monte à 91% pour l’invention et de 7410% pour Treedater. L’erreur de tMRCA maximale est de 179% pour l’invention et de 25 380% pour Treedater. La FIG. 8 illustre le taux d’évolution estimé (à gauche l’invention, à droite Treedater), la FIG. 9 illustre le pourcentage d’erreur tMRCA (à gauche l’invention, à droite Treedater), la FIG. 10 illustre comparativement le pourcentage d’erreur sur les dates pour treedater (en ordonnée) et pour l’invention (en abscisse), et la FIG. 11 illustre également comparativement le pourcentage d’erreur sur les dates (à gauche l’invention, à droite Treedater). La FIG. 12 illustre (avec une échelle logarithmique) le temps d’exécution comparatif entre le procédé de l’invention et treedater, en fonction de la taille de l’échantillon. Les résultats pour Treedater correspondent à l’enveloppe supérieure, les résultats pour l’invention correspondent à l’enveloppe inférieure. On constate que l’écart de temps d’exécution entre le procédé de l’invention et Treedater s’accroît de façon cubique avec la taille de l’échantillon, comme anticipé. La FIG. 13 illustre l’évolution du temps de calcul du procédé de calibration temporelle selon l’invention en fonction de la taille de l’échantillon. Comme anticipé par la théorie, le temps de calcul évolue sensiblement linéairement avec la taille de l’échantillon et permet donc d’envisager une calibration temporelle pour des échantillons de très grande taille. La FIG. 3 illustre schématiquement un exemple d’infrastructure 1 pour la mise en œuvre de l’invention. L’infrastructure repose sur des terminaux 301 à 303 configurés pour récupérer des dates et des informations sur les évènements d’isolation, correspondant typiquement aux feuilles de l’arbre phylogénétique. Les terminaux 301 à 303 récupèrent des informations sur les évènements de d’isolation ou déterminent ces informations à partir d’échantillons. Les terminaux peuvent 301 à 303 recevoir ou mettre en œuvre des séquençages d’ADN ou d’ARN de ces évènements d’isolation. La reconstruction de la phylogénie peut être réalisée par l’intermédiaire de moyens de traitement dédiés ou sur un système de traitement numérique 320 également destiné à la calibration temporelle de l’arbre phylogénétique. A cet effet, les informations provenant des terminaux 301 à 303 peuvent être collectées par le système de traitement approprié par l’intermédiaire d’un réseau de communication 310. Le système de traitement crée de façon connue en soi un arbre phylogénétique dont les feuilles sont les évènements d’isolements du pathogène et dont les nœuds coïncident avec les évènements de transmission du pathogène entre des hôtes lorsque ces nœuds relient des feuilles issues d’hôtes distincts. Chaque feuille représente un pathogène directement observé et chaque nœud représente un évènement de divergence, généralement non directement observé. Les méthodes de datation connues sont basées sur des modèles d’évolution des séquences ARN ou ADN, en exploitant les dates connues des feuilles de l’arbre phylogénétique. On considère généralement que la date d’isolation correspond à chaque feuille. Un ordre de traversée de la racine vers les feuilles est défini dans l’arbre phylogénétique. Les branches de l’arbre phylogénétique ont une longueur correspondant à une distance génétique calculée par le système de traitement approprié. Le système de traitement 320 récupère l’arbre phylogénétique généré. Le système de traitement 320 met en œuvre le procédé de calibration temporelle tel que décrit précédemment, sur la base de cet arbre phylogénétique. La solution déterminée par ce procédé de calibration temporelle est alors stockée dans une base de données 330. Procédé de calibration temporelle d’un arbre phylogénétique par traitement numérique, comprenant : -la réception sous format numérique d’un ensemble d’évènements datés d’isolements ou de transmission d’un agent pathogène et d’un ensemble d’évènements à dater d’isolements ou de transmission de l’agent pathogène, les évènements étant organisés sous forme d’arbre phylogénétique comprenant une racine, des feuilles, des nœuds ordonnés selon un ordre de traversée de la racine vers les feuilles dans l’arbre phylogénétique, la racine, les feuilles et les nœuds étant chacun associé à un desdits évènements, la racine, les feuilles et les nœuds étant reliés par des branches ayant une longueur correspondant à une distance génétique entre des évènements, plusieurs desdits évènements ayant une date à déterminer, chacune desdites branches étant associée à un évènement ascendant et à un évènement descendant de l’arbre phylogénétique selon ledit ordre; -formation d’une matrice Hessienne d’une fonction F d’objectif logarithmique d’un modèle d’horloge moléculaire de l’agent pathogène, définissant la relation entre les distances génétiques, les dates et les durées entre les évènements, la matrice Hessienne H comprenant un bloc T dit local définissant les branches entre lesdits évènements de l’arbre phylogénétique, le bloc T ayant une taille proportionnelle au nombre d’évènements de l’arbre phylogénétique dont la date est à déterminer; -identifier les dates à déterminer dans un vecteur solution w, par une méthode d’optimisation de Newton de résolution de l’équation H(F(v)) s = g(F(v)), avec v un vecteur courant, avec g une fonction gradient, et s une inconnue ou direction de Newton, en fixant le vecteur solution w à la valeur du vecteur courant v lorsque la variation du vecteur courant v entre deux itérations de la méthode de Newton est inférieure à un seuil ou lorsque la variation de la fonction d’objectif F entre deux itérations de la méthode de Newton est inférieure à un seuil. Procédé de calibration temporelle selon la revendication 1, dans lequel la matrice Hessienne H est de la forme : [Math. 7] avec P un bloc dit global correspondant à des paramètres à déterminer du modèle d’horloge moléculaire de l’agent pathogène ; avec C et C’ des blocs correspondants aux dérivées partielles impliquant chacune un évènement et un paramètre global du modèle d’horloge moléculaire de l’agent pathogène. Procédé de calibration temporelle d’un arbre phylogénétique selon la revendication 2, dans lequel au moins un paramètre global du modèle d’horloge moléculaire de l’agent pathogène est à déterminer. Procédé de calibration temporelle d’un arbre phylogénétique selon l'une quelconque des revendications précédentes, dans lequel la résolution du système linéaire dans l’optimisation de Newton comprend une décomposition LU creuse ou une décomposition QR creuse mise en œuvre sur le bloc T. Procédé de calibration temporelle selon l'une quelconque des revendications précédentes, dans lequel les évènements associés aux feuilles sont datés et les évènements associés aux nœuds sont à dater. Procédé de calibration temporelle d’un arbre phylogénétique par traitement numérique selon l'une quelconque des revendications précédentes, dans lequel la longueur de chaque branche est définie en unités de distance génétique ou en nombre de substitutions génétiques. Procédé de calibration temporelle d’un arbre phylogénétique selon l'une quelconque des revendications précédentes, dans lequel le modèle d’horloge moléculaire de l’agent pathogène est un modèle d’horloge relâchée à variance additive, comportant des paramètres définissant un taux d’évolution en unités de distance génétique par unité de temps et une dispersion du taux d’évolution représentative de la variabilité de la vitesse d’évolution entre différentes branches de l’arbre phylogénétique. Procédé de calibration temporelle selon la revendication 5, dans lequel ladite fonction F est de la forme α * log β – log Γ(α) +( α-1) * log d -β* d avec Γ une distribution Gamma ayant les paramètres α=µ *l / (1+ ω) qui correspond à la forme de la distribution et β= 1/ (1+ ω) qui correspond au taux de la distribution, avec ω un paramètre de dispersion du modèle d’horloge relâchée à variance additive. Procédé de calibration temporelle d’un arbre phylogénétique selon l'une quelconque des revendications précédentes, dans lequel dans ledit bloc T, un indice de ligne i correspondant à un unique nœud et la même valeur d’indice de colonne j correspondant à cet unique nœud, les indices de ligne i et de colonne j étant attribués selon l’ordre inverse de la traversée de l’arbre phylogénétique. Procédé de calibration temporelle d’un arbre phylogénétique selon l'une quelconque des revendications précédentes, dans lequel chaque nœud est associé à un unique évènement ascendant et au plus à deux évènements descendants. Procédé de calibration temporelle d’un arbre phylogénétique selon l'une quelconque des revendications précédentes, dans lequel l’amplitude de déplacement dans la direction du vecteur s lors d’une itération de la résolution de Newton est définie pour respecter la règle qu’un évènement descendant est toujours postérieur à son évènement ascendant dans l’arbre phylogénétique. Procédé de calibration temporelle d’un arbre phylogénétique selon l'une quelconque des revendications précédentes, comprenant un calcul d’un niveau d’incertitude de ladite solution par une méthode de Wald, incluant un calcul de la diagonale de l’inverse de la matrice hessienne H de la fonction de vraisemblance appliquée audit vecteur w. Procédé de calibration temporelle d’un arbre phylogénétique selon la revendication 12, dans lequel on utilise une méthode du complément de Schurr pour obtenir la diagonale de l’inverse de la matrice hessienne H lors du calcul du niveau d’incertitude de la solution. Système numérique, caractérisé en ce qu’il comprend des moyens de traitement configurés pour mettre en œuvre un procédé de calibration temporelle selon l'une quelconque des revendications précédentes.