WO2010089363A1

WO2010089363A1 - Procede d'estimation des parametres de la distribution des temps de reponse de particules d'un systeme, applique notamment aux mesures de fluorescence

Info

Publication number: WO2010089363A1
Application number: PCT/EP2010/051404
Authority: WO
Inventors: Tabea Rebafka; François ROUEFF; Antoine Souloumiac
Original assignee: Commissariat a lEnergie Atomique CEA; Telecom ParisTech; Commissariat a lEnergie Atomique et aux Energies Alternatives CEA
Current assignee: Commissariat a lEnergie Atomique et aux Energies Alternatives CEA; Telecom Paris
Priority date: 2009-02-06
Filing date: 2010-02-05
Publication date: 2010-08-12
Anticipated expiration: 2011-08-06
Also published as: FR2942040B1; FR2942040A1

Abstract

Procédé d'estimation des paramètres de la distribution des temps de réponse de particules d'un système, appliqué notamment aux mesures de fluorescence La présente invention concerne un procédé d'estimation de la distribution des temps de réponse de particules d'un système, dans le but notamment de le caractériser. Elle s'applique par exemple aux mesures de fluorescence résolues en temps. Le procédé comporte au moins: -une étape (31) de détermination d'une fonction approchant la log- vraisemblance de la distribution cible ?o correspondant à une distribution sans effet d'empilement à partir des mesures de temps de réponse de particules détectées avec effet d'empilement, c'est-à-dire des mesures de la distribution g?o; -une étape (32) de maximisation de la log-vraisemblance approchée sur l'ensemble des paramètres (?) afin de déterminer l'estimateur comme le vecteur ? qui maximise la log-vraisemblance approchée.

Description

PROCEDE D'ESTIMATION DES PARAMETRES DE LA DISTRIBUTION DES TEMPS DE REPONSE DE PARTICULES

D'UN SYSTEME, APPLIQUE NOTAMMENT AUX MESURES DE FLUORESCENCE

La présente invention concerne un procédé d'estimation de la distribution des temps de réponse de particules d'un système, dans le but notamment de le caractériser. Elle s'applique par exemple aux mesures de fluorescence résolues en temps.

Des mesures de fluorescence permettent de caractériser notamment des molécules. Plus particulièrement, on éclaire des molécules par des lasers, puis les molécules ainsi excitées libèrent des photons. Les temps entre l'instant d'éclairement par le laser et les instants d'arrivée des photons sur un détecteur permettent de caractériser ces molécules de façon connue. Le nombre de photons émis lors de chaque éclairement est indéterminé, à titre indicatif plus l'impulsion émise par le laser est forte, plus il y a de photons. Il y a par ailleurs un aléa supplémentaire sur les temps de réponse. Un problème technique général est la mesure de la distribution des temps de réponse d'un système donné dans le but de le caractériser à l'image de l'exemple précédent.

La distribution de ces temps de réponse est généralement paramétrée par un ensemble de paramètres numériques et ce sont ces paramètres que l'on veut mesurer. On notera par la suite cette dépendance envers un jeu de paramètre en appelant la densité de probabilité de la distribution des

temps de réponse, θo désignant un vecteur de paramètres appartenant à un ensemble Θ de vecteurs θ. A tout élément θ dans Θ on associe une densité de probabilité et ainsi on définit une famille de distribution

D'un point de vue statistique, cette approche est un problème d'estimation paramétrique. Par exemple, dans le contexte de la fluorescence résolue en temps, où le système étudié est alors un ensemble de molécules fluorescentes, on suppose que la fonction de distribution, paramétrée par un jeu de paramètres, est un mélange d'exponentielles et

est l'ensemble des mélanges de densités exponentielles. La modélisation est plus complexe lorsqu'il s'agit de la distribution des temps de réponse d'un système optique ou d'un milieu transparent comme l'atmosphère. Cependant l'extension de cette approche paramétrique à d'autres domaines que la fluorescence résolue en temps est possible.

Les mesures des temps de réponse sont obtenues au moyen d'une excitation puisée du système à caractériser. On émet une excitation très brève, on mesure le temps mis par le système à répondre à cette excitation et on enregistre ce temps. Ces trois opérations sont répétées un très grand nombre de fois jusqu'à ce qu'on dispose de milliers ou de millions de temps de réponse. Enfin, c'est à partir de cette collection de temps de réponse collectés indépendamment que l'on estime les paramètres de leur distribution commune.

Le domaine de l'invention concerne les contextes où le système peut émettre un nombre aléatoire de réponses à l'excitation et où seule la première réponse est détectée et enregistrée. Dans le contexte de la fluorescence résolue en temps, l'excitation est une impulsion laser envoyée sur l'échantillon fluorescent et la réponse de cet échantillon est constituée en général de plusieurs photons, car il y a un très grand nombre de molécules fluorescentes présentes. On enregistre seulement le temps de réponse correspondant au premier photon car ce premier photon aveugle l'appareillage électronique de détection et d'enregistrement. Ce phénomène d'oubli des deuxième, troisième photons et des suivants est appelé empilement ou « pile-up » dans la littérature anglo-saxonne. L'invention concerne donc les contextes qui comportent notamment les quatre caractéristiques suivantes : - Mesure de la distribution de temps de réponse d'un système,

- Distribution paramétrée,

- Mise à disposition d'une collection généralement grande de temps de réponse mesurés indépendamment,

- Le système étudié émet un nombre aléatoire de réponses à chaque excitation mais seule la première réponse est vue par l'appareillage de mesure.

Ces quatre conditions sont rencontrées dans la méthode de mesure dite « Time Correlated Single Photon Counting » (TCSPC) en fluorescence résolue en temps. Des données du type « pile up » sont rencontrées en fluorescence résolue en temps. Il s'agit d'estimer la distribution des durées de vie en fluorescence. La constante de durée de vie d'une molécule est le temps moyen entre l'instant d'excitation de la molécule et l'instant d'émission d'un photon. Elle contient des précieuses informations sur les molécules et leur milieu, par exemple le pH, la polarisation, ou la viscosité. Elle est donc de grand intérêt en biologie, en médecine ou en physique notamment. Typiquement, la distribution de la durée de vie d'une molécule est modélisée par une loi exponentielle. Or, la distribution des temps d'arrivée d'un mélange des molécules correspond à un mélange des lois exponentielles. En outre, le nombre de photons émis par impulsion de laser d'excitation est modélisé par une loi de Poisson.

Les observations disponibles pour l'estimation des paramètres exponentiels peuvent être obtenues par la méthode TCSPC, aussi appelée la méthode SPT (Single Photon Timing), comme le montrent notamment les documents de B.Valeur « Molecular Fluorescence », WILEY-VCH, Weinheim, 2002, de J.R.Lakowicz « Principles of Fluorescence Spectroscopy », Academic/Plenium, New- York, 1999, ou de D.O'Connors et al « Time- Correlated Single Photon Counting », Académie Press, London, 1984. Dans la méthode TCSPC un coup de laser excite un nombre aléatoire de photons, mais pour des raisons techniques seulement le temps de réponse du photon qui arrive le premier au détecteur peut être enregistré. Les autres photons fluorescents ne sont pas observables. En TCSPC le réglage de l'intensité de laser, ou l'atténuation de la chaîne optique, détermine le nombre moyen de photons détectés par excitation.

Un problème technique est notamment le problème inverse. Plus particulièrement, on veut estimer les paramètres d'une densité de probabilité M) , mais on ne dispose pas d'observations directes de la densité M> Pourtant on observe le minimum d'un nombre aléatoire des variables aléatoires indépendantes de la densité hh . Il est évident que la densité des observations, notée ->'^, n'est pas la même que la densité de probabilité M« . En effet, $&$ résulte de h'h par une distorsion non linéaire, effet de l'empilement « pile up ». On ne peut pas appliquer les mêmes estimateurs aux données sans « pile up », issues de la densité M« , qu'aux données avec pile up, issues de la densité ^. Le problème d'estimation vient notamment de la forme compliquée de la densité

En effet, les estimateurs classiques, comme l'estimateur du maximum de vraisemblance ou les méthodes des moments, ne sont pas faisable.

Il existe au moins deux façons d'affronter le problème inverse où une observation représente le minimum d'un nombre aléatoire de variables aléatoires. En TCSPC il est courant de travailler avec une intensité de laser tellement faible qu'il y a au plus un photon émis par excitation avec une très forte probabilité. Donc, le problème d'une perte de photons, c'est-à-dire de ne pas observer tous les photons qui arrivent sur le détecteur, est évité. Dans ce cas, les temps d'arrivée observés peuvent être considérés comme des réalisations indépendantes d'un mélange exponentiel. Des estimateurs de paramètres du mélange exponentiel sont obtenues par des méthodes statistiques classiques comme la méthode des moindres carrées décrite dans le document de B.Valeur et al précité, ou la méthode du maximum d'entropie décrite dans le document de J.C.Brochon « Maximum entropy method of data analysis in time-resolved spectroscopy » Methods Enzymology, 240 : 262-311 , 1994.

Or, une très basse intensité de laser entraîne un long temps d'acquisition, car sur 100 excitations il n'y a qu'un à trois photons détectés en général. Les autres excitations ne sont pas suivies par la détection d'un photon, donc elles ne contiennent aucune information sur les paramètres à estimer. Il faut alors un grand nombre d'excitations pour obtenir une quantité suffisante de temps d'arrivée pour une estimation fiable des paramètres. Une étude de la borne de Cramér-Rao montre qu'une minuscule intensité de laser entraîne une forte augmentation de la vahance des estimateurs et qu'il est préférable de travailler avec un laser plus fort que ceux utilisés habituellement afin d'augmenter l'information sur les paramètres qui est contenue dans les données. Cependant, avec une plus forte intensité, les données subissent l'effet de « pile up » et il faut des méthodes statistiques adaptées au modèle avec « pile up ».

La littérature sur le modèle statistique de « pile up », où uniquement le minimum d'un nombre aléatoire de variables aléatoires est observé, est très limitée et, en TCSPC, ces méthodes sont rarement utilisées dans la pratique. Une première méthode, décrite dans le document de P.B.Coates « The correction of photon 'pile-up' in the measurement of relative lifetimes » Journal of Physics E : Scientific Instruments, 2(1 ) : 878-879, 1968, est une correction de l'histogramme pour des données en TCSPC. Il s'agit de repondérer les groupes ou classes de l'histogramme, encore appelés « bins », afin d'obtenir un histogramme de la distribution cible

Une généralisation de l'approche de P.B.Coates à une grande classe de distributions est présentée dans le document de A.Souloumiac « Apparatus and method for characterizing a System for counting elementary particles », demande de brevet WO2007088173(A1 ). En plus, le document de A.Souloumiac étend la méthode de la repondération à la fonction de répartition afin de corriger l'effet de « pile up ». Dans le document de J.G.Walker « Itérative correction for 'pile-up' in single-photon lifetime measurement », Optics Communications, 201 : 271 -277, 2002, la méthode de P.B.Coates est adaptée au cas d'un laser non stable. Un inconvénient de ces approches de repondération est qu'elles ne fournissent pas d'estimateurs de paramètres. Après la repondération de l'histogramme ou de la fonction de répartition, il faut appliquer une autre méthode statistique pour arriver à des valeurs estimées des paramètres. La repondération ne représente qu'une première étape dans l'estimation. Des méthodes d'estimation de paramètres directe pour le modèle « pile up » ont été développées. Le document de L.Okano et al « Abstracting reliable parameters from time-correlated single photon counting experiments », Journal of Photochemistry and Photobiology A : Chemistry, 175, 221 -225, 2005 propose une méthode des moindres carrées pour ajuster directement la fonction de densité aux données avec empilement. Les méthodes des

moindres carrées sont des « large sample methods » ce qui veut dire que l'on a besoin d'un très grand nombre de temps de réponse afin d'obtenir des estimateurs fiables. Comme nous cherchons à diminuer le temps d'acquisition, ce procédé n'est pas adapté. Une autre méthode directe est présentée dans le document de T. Rebafka « An MCMC approach for estimating a fluorescence lifetime with pile-up distorsion », Actes 21 ^ιeme Colloque GRETSI sur le Traitement du Signal et des Images (GRETSI 2007), Troyes, France, où un échantillonneur de Gibbs est développé qui est bien adaptée au modèle de mélange exponentiel et à un nombre Poissonien de photons par excitation. Cette méthode présente cependant plusieurs inconvénients, en particulier elle nécessite un long temps de calcul et la difficulté d'inclure la fonction d'appareil dans le modèle. Elle ne peut donc pas être utilisée en TCSPC où les mesures contiennent généralement du bruit additif.

Enfin les deux méthodes précédentes ne sont pas très souples dans le sens qu'elles sont restreintes au cas multi-exponentiel. Même si la loi exponentielle est la plus utilisée en fluorescence, il y a des applications avec des distribution de Kohlrausch, comme décrit dans le document de M.N.Berberan-Santos et al « Mathematical functions for the analysis of luminescence decays with underlying distribution : 1. Kohlrausch decay function (stretched exponential), Chemical Physics, 315 : 171-182, 2005, ou des distributions de Becquerel comme décrit dans le document de M.N.Berberan-Santos et al « Mathematical functions for the analysis of luminescence decays with underlying distribution : 2. Becquerel (compressed hyperbola) and related decay functions, Chemical Physics, 317 : 57-62, 2005.

Un but de l'invention est notamment de pallier les inconvénients précités en permettant d'obtenir une estimation directe des paramètres de la distribution des temps, adapté à l'effet de « pile-up », ne nécessitant donc pas l'atténuation de la puissance du laser utilisé. A cet effet, l'invention a pour objet un procédé d'estimation des paramètres θo de la distribution des temps de réponse de particules réémises par un système, en réponse à une série d'émissions de particules vers ledit système, un temps de réponse d'une particule correspondant étant mesuré entre un instant d'émission et l'instant de détection de ladite particule, le procédé comportant au moins :

- une étape de détermination d'une fonction approchant la log- vraisemblance de la distribution cible correspondant à une

distribution sans effet d'empilement à partir des mesures de temps de réponse de particules détectées avec effet d'empilement, c'est-à-dire des mesures de la distribution

- une étape de maximisation de la log-vraisemblance approchée sur l'ensemble des paramètres

afin de déterminer l'estimateur comme le vecteur θ qui maximise la log-vraisemblance approchée.

La fonction approchée est une estimation de la log-vraisemblance de la fonction de distribution cible basée sur des mesures observées

distribuées selon la fonction de distribution ^ = avec effet d'empilement. Supposons que Zi,...,Z_n représentent les n mesures de temps de réponse. La fonction approchée est obtenue en pondérant chaque terme log(f_θ(Zi)) de la log-vraisemblance par un poids qui compense le fait que Z, suit une loi ^- > et non à cause du phénomène du « pile-up ». On peut prendre un poids

de la forme 1

où m note la fonction génératrice de N,

et m' la dérivée de m et &u est un estimateur de la fonction de répartition. Si l'on choisit l'estimateur empirique de la fonction de répartition pour ce dernier et si l'on réordonne les observations dans un ordre croissant sous la forme

alors le poids

correspondant au terme log(f_θ(Z_(l))) s'écrit

La fonction de distribution cible

est par exemple considérée comme étant un mélange de densités exponentielles.

Avantageusement, dans un mode de mise en œuvre particulier, la fonction de distribution cible prend en compte la fonction d'appareil de mesures, la fonction cible étant le produit de convolution de la fonction cible sans fonction d'appareil (h\-> ) avec la densité de probabilité de la fonction d'appareil η.

Les particules sont par exemple des photons émis par un laser, le système réémettant des photons.

Le système peut être un ensemble de molécules à caractériser par la fonction de distribution cible des temps de réponse.

D'autres caractéristiques et avantages de l'invention apparaîtront à l'aide de la description qui suit faite en regard de dessins annexés qui représentent : - la figure 1 , une illustration de l'effet de « pile-up » sur les distributions de mesures de temps de réponse ; - les figures 2a, 2b et 2c, des illustrations d'histogrammes de mesures avec effet de « pile-up » pour différentes intensités d'émissions de photons ;

- la figure 3, une présentation des étapes possibles pour la mise en œuvre d'un procédé selon l'invention.

La figure 1 illustre par deux courbes 1 , 2 l'effet de « pile-up » évoqué précédemment dans un système d'axes, où les abscisses représentent les temps de réponse mesurés et l'axe des ordonnées le nombre de photons détectés correspondant à ces temps de réponse. Les valeurs indiquées sont normalisées. La première courbe 1 et la deuxième courbe 2 illustrent respectivement les densités de probabilité

et évoquées

précédemment.

Dans la méthode TCSPC on éclaire donc un échantillon de molécules avec un rayon laser par exemple. Les molécules excitées par des photons issus du laser émettent à leur tour des photons qui sont détectés par des moyens appropriés, utilisant notamment des filtres adaptés. On mesure et on enregistre le temps d'arrivée des photons. On obtient ainsi une distribution des temps de réponse caractérisant les molécules. Les temps étant de l'ordre de plusieurs nanosecondes pour une précision exigée de l'ordre d'une picoseconde.

En réalité, il n'y a pas un seul photon. Il y a en fait deux aléas. Le premier aléa porte sur le nombre de photons réémis par les molécules et le deuxième aléa porte sur les temps de réponse, car chaque photon a son temps de réponse associé qui lui est propre et donc indéterminé. Plus le laser émet des impulsions fortes, plus il y a de photons réémis par les molécules. Les temps étant très courts, le chronomètre utilisé ne peut mesurer qu'une seule durée, la durée correspondant au temps du premier photon reçu. Le temps que le chronomètre se réinitialise pour engager une deuxième est alors trop long en regard du phénomène physique observé. La fluorescence est en effet terminée et tous les photons déjà reçus. A chaque impulsion laser, on a donc un certain nombre de photons qui arrivent mais seul le temps de réponse du premier est mesuré. La distribution obtenue n'est donc construite que pour le premier photon reçu, et donc pour le temps de réponse le plus court, conduisant ainsi à une image compressée et finalement à une distorsion de mesure. C'est l'effet d'empilement ou « pile- up », où seul le premier photon est détecté et mesuré. La fonction de densité associée aux mesures est illustrée par la première courbe 1. Cette courbe 1 présente une forte concentration de mesures correspondant aux temps de réponse très courts.

Une solution connue consiste à baisser la puissance du faisceau laser émis de façon à ce qu'il n'y ait qu'un seul photon réémis à chaque impulsion laser. La distribution obtenue est satisfaisante car dénuée de distorsion. La deuxième courbe 2, de forme exponentielle, illustre la densité associée à cette solution. La différence de la première courbe 1 par rapport à cette deuxième courbe 2 illustre la distorsion provoquée par l'effet de « pile-up ». L'avantage des mesures correspondant à la courbe 2 est que l'on sait comment en extraire des paramètres θo recherchés avec des méthodes statistiques, puisqu'il s'agit d'une densité exponentielle ou dans des cas plus complexes d'un mélange exponentiel. Cependant il s'ensuit que dans la grande majorité des cas, pratiquement dans environ 95% impulsions, il n'y a pas de photons réémis. Seules 2% à 5% des impulsions laser émises donnent lieu à la réémission d'un photon. L'obtention d'une distribution selon la deuxième courbe 2 se fait donc avec une perte de temps important. Il existe des applications où la distribution doit être obtenue rapidement et où les pertes de temps de la méthode illustrée par cette deuxième courbe 2 n'est pas admissible. En particulier, dans les cas où le système observé est instable et où plus la mesure est longue, plus il y a de possibilités que l'état du système évolue rendant l'ensemble des mesures inaptes à caractériser le système. Il y a aussi des applications qui nécessitent un haut débit de mesures, notamment dans le domaine médical.

Les figures 2a, 2b et 2c illustrent les histogrammes de mesures 21 , 22, 23 avec l'effet de « pile-up » pour des intensités différentes. Comme indiqué précédemment, en TCSPC le réglage de l'intensité de laser, ou l'atténuation de la chaîne optique, détermine le nombre moyen de photons détectés par excitation. L'effet des choix différents de l'intensité sur la distribution des observations est illustré par les trois histogrammes des figures 2a, 2b et 2c correspondant à 1000 observations dans chacun des cas. La densité cible J ^fA> est une exponentielle de moyenne 1 et chaque histogramme est composé de 1000 observations, c'est-à-dire 1000 impulsions de laser, obtenues par simulation Monte-Carlo. Le dernier bin, noté +°°, représenté sur l'échelle des abscisses, contient le nombre d'observations où aucun photon n'est détecté. L'intensité du laser est exprimée par le nombre moyen λ de photons détectés par excitation. Plus précisément, λ est le paramètre de la loi de Poisson qui est la distribution du nombre de photons détectés par excitation. Quand l'intensité est faible, λ = 0,02, il n'y a que 13 photons détectés sur 1000 excitations, cas illustré par la figure 2a. Quand l'intensité est très élevée, λ = 20, la distribution est presque dégénérée car presque toutes les observations se trouvent dans le premier bin, cas illustré par la figure 2b. Dans ces deux cas extrêmes, il est clair que les histogrammes contiennent très peu d'information sur le paramètre exponentiel, de la distribution des temps de réponse, à estimer. Il semble favorable d'utiliser une intensité moyenne. La figure 2c montre que λ = 2 mène à un histogramme « habituel » car les observations où un photon a été détecté sont très nombreuses et les temps d'arrivée sont assez dispersés.

Un but de l'invention est d'obtenir une distribution des temps de réponse conforme à la première courbe 1 , sans perte de temps, c'est-à-dire obtenue beaucoup plus rapidement que la distribution selon la deuxième courbe 2. En particulier, la puissance des faisceaux laser émis n'est pas diminuée autant que pour la deuxième courbe 2. C'est notamment la résolution du problème inverse évoqué précédemment, conduisant à l'obtention des paramètres de la distribution selon la deuxième courbe 2. Dans le cas d'une simple exponentielle, par exemple, le paramètre obtenu est le paramètre τ qui caractérise la densité exponentielle

A cet effet, l'invention comporte au moins deux étapes telles qu'illustrées par la figure 3. Dans une première étape 31 on détermine une fonction simple approchant la log-vraisemblance de la distribution cible, qui est la fonction recherchée, illustrée par la deuxième courbe 2 de la figure 1. La log- vraisemblance, dépendant des mesures de temps de réponse, est une fonction en

le vecteur des paramètres qui caractérise la famille des distributions à laquelle appartient la densité cible

. Puis dans une deuxième étape 32, on maximise cette log-vraisemblance approchée en θ. Rechercher les paramètres θo de la fonction î

revient en fait à rechercher le vecteur θ des paramètres où la log-vraisemblance atteint son maximum. La maximisation est faite sur tout l'ensemble Θ des paramètres θ qui détermine la famille des distributions à laquelle appartient la densité cible

Ainsi dans la première étape 31 , on construit un estimateur de la fonction de log-vraisemblance E[log(fβ(Y))] basée sur des données selon la densité ^, E étant la fonction d'espérance mathématique et Y étant une variable aléatoire de densité h.^ c'est-à-dire E[log(fθ(Y))] = J log(fθ(y)) fβ(y)dy. Soient Zi, ...,Z_n des mesures des temps de réponse des photons. On note les mêmes mesures ordonnées par ordre croissant par Z₍₁₎, ... ,Z_(n). C'est-à-dire que Zf₁) désigne le plus petit temps de réponse mesuré, Z₍₂) le deuxième plus petit et ainsi de suite. L'estimation de la log-vraisemblance est obtenue par une repondération de la log-vraisemblance

de la densité fβ afin de corriger la distorsion, W^_n représentant le poids associé à la donnée Z_(i). Dans cette phase de repondération, on donne un poids aux données Z_(i) qui dépend de leurs positions temporelles. En d'autres termes, les données Z^ sont ordonnées en fonction de leur valeur et leur poids dépend de leur ordre, c'est-à-dire que le poids des photons détectés Z, dépend de leur ordre d'arrivée. En particulier, le temps d'arrivée Z_(n) le plus élevé est associé au plus grand poids dans la fonction

Dans la deuxième étape 32, on maximise la log-vraisemblance ainsi corrigée en utilisant un algorithme adapté, comme l'algorithme EM, abbréviation pour « Expectation-Maximization », présenté dans le document de A.P.Dempster et al, « Maximum likelihood from incomplète data via the EM algorithm », Journal of the Royal Statistical Society, Séries B, 39(1 ), 1 -38, 1977.

On présente par la suite un exemple de méthode d'estimation pour une application avec effet de « pile-up » dans le cas simple d'un mélange de densités exponentielles. Puis on présente une méthode d'estimation pour une application avec une fonction d'appareil.

Dans une approche statistique avec « pile up » pour des mélanges exponentiels, on suppose que la distribution des durées de vie en fluorescence représente un mélange de lois exponentielles. On considère la famille de densités {f_θ, θeΘ} donnée par la relation (1 ) suivante :

où > κ>0 sont les paramètres exponentiels et les poids satisfont

et

\_e vecteur de paramètres et Θ l'ensemble de paramètres θ. Les constantes de durées de vie en fluorescence sont données par les paramètres exponentiels

En plus, K désigne l'ordre du mélange exponentiel qui correspond au nombre de molécules différentes. On suppose que K est connue et fixe. On note

la fonction de répartition associée à la densité donnée par

On note Y1 Y2, ... une suite de variables aléatoires indépendantes de densité où

est le paramètre inconnu à estimer. Les

représentent les temps d'arrivée des photons de fluorescence qui frappent le détecteur.

Le nombre N de photons qui frappent le détecteur par impulsion laser est modélisé par une loi de Poisson de paramètre λ (restreinte à {

. Les probabilités sont données par la relation suivante :

On suppose que N est indépendant de la suite Y1 Y2, .... On définit alors l'observation avec « pile up » par

la mesure retenue Z correspond, comme indiqué précédemment, au temps de réponse le plus court. La fonction de répartition avec « pile-up » '"^:! de la variable aléatoire Z est donnée par la relation suivante :

et la densité avec « pile up » -^ par la relation suivante

Pour effectuer l'estimation, la méthode selon l'invention approche la log- vraisemblance de

par un estimateur basé sur des observations avec « pile up ». Ensuite, on obtient un estimateur du paramètre θo en maximisant cette log-vraisemblance approchée. La forme simple de cette log-vraisemblance permet l'application de l'algorithme EM, défini par la suite, pour déterminer son maximum.

On procède ensuite à l'approximation de la log-vraisemblance de fe. La log-vraisemblance de la fonction f_θ associée à un échantillon Yi, ...,Y_m de variables aléatoires indépendantes et identiquement distribuées (i.i.d) de densité ^

est donnée par la relation suivante :

et l'estimateur du maximum de vraisemblance classique Θ_ML est défini par la relation suivante :

II découle de la loi forte des grands nombres que L₀(θ) est un estimateur consistent de E[log(fe(Y))]. Si Z est une variable aléatoire de densité -^ * , on établit la propriété suivante en utilisant relation (3),

On en déduit que la quantité

est un estimateur consistant de E[log(fo(Y))] basé sur un échantillon Z₁, ... ,Z_n i.i.d. de densité

Pourtant, la quantité L₁(O; θ₀) ne peut pas être utilisée pour estimer le paramètre inconnu θo, car elle-même dépend de θo. Or, une autre modification est indispensable. On déduit de la relation (2) que :

En remplaçant

Λ; i, par la fonction de répartition empirique

on établit un nouvel estimateur de

E[log(f_θ(Y))] similaire à L i(θ; θo) par

où Z(ι) désigne la /-ème statistique d'ordre de Z₁, ... ,Z_n. On introduit les poids

On définit un nouvel estimateur de θo par la relation suivante

Donc, on cherche le paramètre où la log-vraisemblance approchée atteint son maximum. On peut considérer comme un estimateur de maximum de vraisemblance approché.

On remarque que l'on peut également utiliser d'autres estimateurs de la fonction de répartition que

Cela changera les poids qui sont associés

aux données mais le principe restera le même.

On effectue par la suite une maximisation par l'algorithme EM.

On constate que le problème de maximisation de la relation (5) diffère du problème de la relation (4) seulement par des constantes positives w

Or, la maximisation de la relation (4) peut se faire par l'algorithme EM. Grâce à la structure semblable des problèmes, l'algorithme EM s'applique également à la relation (5).

On introduit la variable manquante S d'une observation Y du mélange exponentiel fe étant l'étiquette de la composante du mélange qui a généré l'observation Y. Plus précisément, S est une variable aléatoire à valeurs dans La loi jointe de (Y₁S) est donnée par

Comme dans l'algorithme EM classique, on définit la fonction

et on calcule itérativement la suite (θ^(t))_t en maximisant Q(θ,(θ^(t)) , conduisant à la relation suivante :

Ce problème de maximisation est facile car il se décompose en plusieurs problèmes de maximisation indépendants:

où la dernière maximisation est faite sous la contrainte que tous les ctk sont à valeurs dans ]0,1[ et ∑kCik = 1. On peut montrer qu'avec chaque itération la log-vraisemblance approchée L(θ) augmente. Plus précisément on a pour toute valeur initiale θ^<1) e θ et pour tout t ≥1,

Par ailleurs, on peut montrer que la suite converge vers un point

critique de L(θ).

Il y a plusieurs critères d'arrêt pour déterminer le moment quand la suite

a convergé et où on peut arrêter le calcul itératif. Un critère repose sur l'évolution de la log-vraisemblance approchée. On calcule la valeur de L(θ^(t)) à toute itération et on s'arrête si l'accroissement par rapport à l'étape précédente est inférieur à un certain seuil. Un autre critère étudie la variation des valeurs du paramètre θ^m. A toute itération on évalue la différence

désigne la norme infinie. Si la différence est inférieure à un certain seuil, on considère la dernière valeur de θ^m comme limite de la suite

L'algorithme EM pour maximiser la fonction L(O). dans le cas d'un mélange exponentiel, peut être défini par les phases suivantes :

On fixe un seuil σ > 0 pour le critère d'arrêt.

Soit λ la valeur connue du paramètre Poissonien lors de l'expérience. Soit _n les observations avec « pile-up », c'est-à-dire un échantillon

i.i.d. de densité

On note Z

les observations ordonnées, n étant le nombre total de mesures avec photons détectés, et / étant l'ordre d'arrivée du photon

1. On pose t = 1.

On choisit une valeur initiale

On calcule, pour / = 1, ...,n

2. On réitère jusqu'à la convergence les étapes suivantes :

(a) Calculer les probabilités conditionnelles, pour i = 1,...,n et s = 1, ...,K ,

(b) Déterminer les nouveaux paramètres, pour s

On pose

(c) D'après le critère d'arrêt choisi, on évalue

(d) Si σ_t+i < σ, on passe à la phase 3 ci-dessous qui donne l'estimateur,

(e) dans le cas contraire, on pose t = t+1 et on retourne à l'étape 2(a). 3. L'estimateur du paramètre θo est donné par

Le modèle précédent n'est que proche de la réalité pour des appareils de mesure avec un laser qui émet de la lumière instantanément, c'est-à-dire si l'impulsion laser peut être modélisée par une loi de Dirac avec toute la masse en zéro. Cependant, la source de lumière de la plupart des appareils de mesure émet de la lumière pendant une courte durée. On appelle fonction d'appareil la distribution de rémission de lumière qui est différente pour tout appareil, mais qui peut être enregistrée facilement avant d'effectuer l'expérience. Donc une observation TCSPC est le temps que la molécule reste en état excité, la durée de vie, à laquelle s'ajoute une durée aléatoire qui provient de l'appareil. La méthode EM prend en compte les aléas dus à l'appareil.

On note η la densité de probabilité de la fonction d'appareil. On supposera dans la suite que η est connue. Soient (E₁, i≥1} et (Y₁, i≥1} deux suites indépendantes de variables aléatoires indépendantes de loi η respectivement dduu mmééllaannggee eexxppoonneennttiieell

donnée par la relation (1 ).

On pose

La densité de est la convoi ution de

On note la fonction de répartition de

par ' Or, une observation en

TCSPC a la même distribution que : .

La fonction de répartition avec « pile-up » àe est donnée par :

et la densité avec « pile-up » 5Kest donnée par

De la même manière que précédemment, on déduit l'estimateur comme le vecteur θ qui maximise la log-vraisemblance approchée,

où

désigne des observations avec « pile-up » de densité et

la log-vraisemblance approchée est donnée par

Pour l'observation

on considère les deux variables manquantes (S, E) où S est l'étiquette de la composante du mélange qui a généré Y, et E désigne la durée ajouté à la durée de vie par l'appareil. La loi jointe de

est donnée par :

L'algorithme EM pour maximiser

dans le cas d'une fonction d'appareil est défini par les phases suivantes :

On fixe un seuil σ>0 pour le critère d'arrêt.

Soit λ la valeur connue du paramètre Poissonien lors de l'expérience. Soit η une densité connue constante par morceaux donnée par

Soient les observations avec pile up, un échantillon i.i.d. de

densité .9^ . On note ^ les observations ordonnées

1 . On pose t = 1 On choisit une valeur initiale ^v

On calcule, pour i=1, ...,n ,

2. On réitère jusqu'à convergence les étapes suivantes : (a) On calcule, pour i=1,...,n et s=1, ...,K,

où aΛb dénote le minimum de a et b, et on calcule

et

(b) On détermine les nouveaux paramètres, pour s=1, ...,K,

et

On pose ^t

(c) D'après le critère d'arrêt choisi, on évalue

(d) Si

passer à la phase (3) ci-dessous.

(e) Dans le cas contraire on pose t = t+1 et on retourne à l'étape 2(a).

3. L'estimateur du paramètre θ₀ est donné par

L'invention utilise l'algorithme précédent pour le problème inverse général évoqué précédemment, c'est-à-dire pour des familles de densité autres que la famille de mélanges exponentiels et pour d'autres lois discrètes que la loi de Poisson. Ce nouvel algorithme EM, utilisé pour résoudre le problème inverse général, n'est donc pas restreint au mélange exponentiel et la loi de Poisson. On considère

une famille de densités quelconque telle que l'algorithme EM s'applique à un échantillon i.i.d. de densité

Soient

une suite de variables aléatoires indépendantes de densité On note

la fonction de répartition associée à la densité

On considère N^* une variable aléatoire à valeurs dans {1,2, ...} indépendante de la suite i

On note m la fonction génératrice de A/^*,

et m' la dérivée de m. On définit l'observation avec « pile-up » par :

On établit la log-vraisemblance approchée basée sur des observations

_{avec (<} pjie-up » :

OÙ

On définit

où S dénote les variables manquante associées à la variable

et désigne la densité jointe de ^ - ^ k

L'algorithme EM général pour maximiser comporte les phases

suivantes :

On fixe un seuil σ>0 pour le critère d'arrêt. 1. On choisit une valeur d'initialisation On pose t = 1.

2. On réitère jusqu'à la convergence les étapes suivantes : (a) On poser t = t+1. (b) On calcule

(c) D'après le critère d'arrêt choisi, on évalue

(d) Si σt < σ on passe à la phase (3) ci-dessous, sinon on retourne à l'étape 2(a).

3. L'estimateur du paramètre θo est donné par

Une méthode d'estimation selon l'invention peut fournir des mesures très précises quand elle est appliquée à des données réelles. Elle demande peu de temps de calcul, prend en compte la fonction d'appareil, est facile à mettre en œuvre et peut traiter le cas d'un mélange de plusieurs exponentielles.

En comparaison avec des méthodes existantes, elle est notamment adaptée à la distorsion non linéaire du problème général. En particulier, on peut appliquer les résultats indiqués dans Rebafka et al, 2008, concernant le réglage optimal de l'intensité du laser. En d'autres termes, l'expérience peut être effectuée sous des conditions optimales et ainsi une réduction du temps d'acquisition, d'un facteur 10 par exemple, peut être atteinte par rapport aux techniques actuellement utilisées.

Le système à caractériser peut être composé de molécules comme décrit précédemment. L'invention peut néanmoins s'appliquer à d'autres systèmes. Un système à caractériser peut être par exemple une fibre optique dans laquelle on injecte des photons et dont mesure les durées de vie. D'autres types de particules que les photons peuvent aussi être utilisés.

Claims

REVENDICATIONS

1. Procédé d'estimation des paramètres θo de la distribution des temps de réponse de particules réémises par un système, en réponse à une série d'émissions de particules vers ledit système, un temps de réponse d'une particule correspondant étant mesuré entre un instant d'émission et l'instant de détection de ladite particule, caractérisé en ce que le procédé comporte au moins :

- une étape (31 ) de détermination d'une fonction approchant la log- vraisemblance de la distribution cible <M> correspondant à une distribution sans effet d'empilement à partir des mesures de temps de réponse de particules détectées avec effet d'empilement, c'est-à-dire des mesures de la distribution ^o ;

- une étape (32) de maximisation de la log-vraisemblance approchée sur l'ensemble des paramètres Θ afin de déterminer l'estimateur comme le vecteur θ qui maximise la log-vraisemblance approchée.

2. Procédé selon la revendication 1 , caractérisé en ce que la fonction approchée est une estimation de la log-vraisemblance de la fonction de distribution cible <M> basée sur des mesures observées distribuées selon la fonction de distribution ^ avec effet d'empilement.

3. Procédé selon la revendication 2, caractérisé en ce que la fonction approchée est obtenue par une pondération de la log-vraisemblance de la fonction de distribution cible J** .

4. Procédé selon la revendication 3, caractérisé en ce que la pondération de la log-vraisemblance se fait de façon à ordonner les mesures de temps de réponse observées Zi, ...,Z_n en fonction de l'ordre d'arrivée des particules correspondante, le poids affecté à chaque mesure dépendant de son ordre.

5. Procédé selon la revendication 4, caractérisé en ce que le poids affecté à une mesure croît avec son ordre d'arrivée associé.

6. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que n étant le nombre total de particules détectées et / étant l'ordre d'arrivée d'une particule, le poids w^* _hn associé à la mesure du temps de réponse Z_(l) d'une particule d'ordre / est donné par la relation suivante :

où m note la fonction génératrice de N,

et m' la dérivée de m et est un estimateur de la fonction de répartition.

7. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que la fonction de distribution cible est considérée

comme étant un mélange de densités exponentielles.

8. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que la fonction de distribution cible prend en compte la fonction d'appareil de mesures, la fonction cible étant le produit de convolution de la fonction cible sans fonction d'appareil

avec la densité de probabilité de la fonction d'appareil η.

9. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que les particules sont des photons émis par un laser, le système réémettant des photons.

10. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que le système est un ensemble de molécules à caractériser par la fonction de distribution cible des temps de réponse.