DE19923588C2 - Verfahren zur Erfassung und Auswertung von Messdaten und zur Durchführung des Verfahrens geeigneter Computer sowie Logikbaustein - Google Patents

Description

Die Erfindung betrifft ein Verfahren zur Erfassung und Auswertung von Messdaten, insbesondere von Messdaten aus der Kernresonanztomographie, wobei die Messdaten erfasst und mit einer Modellfunktion verglichen werden.

Die Erfindung betrifft ferner einen zur Durchführung des Verfahrens geeigneten Computer sowie einen entsprechend geeigneten Logikbaustein.

Ein gattungsgemäßes Verfahren ist aus dem Artikel von Peter A. Bandettini et. al.: Processing Strategies for Time-Course Data Sets in Functional MRI of the Human Brain, Magnetic Resonance in Medicine, 30, p. 161-173, 1993, bekannt. Hierbei wird ein Signalverarbeitungsverfahren eingesetzt, durch das ein zeitlicher Verlauf von Signalanregungen mit einem Paradigma, das heißt einem zeitlichen Ablauf einer Aufgabenreihe für einen Probanden, verglichen wird. In diesem Artikel wird ferner eine Korrelationsanalyse dargestellt, die für einzelne Echosignale erfolgt und die für mehrere, zeitlich nacheinander erfolgende, Relaxationsmessungen durchgeführt wird.

Stand der Technik nach § 3 Abs. 2 PatG ist die Deutsche Patentschrift DE 198 17 228 C1, in der ein Computer zur Auswertung von Daten aus der kernmagnetischen Resonanztomographie sowie ein Verfahren zur Auswertung von Daten aus der kernmagnetischen Resonanztomographie dargestellt sind. Hierbei werden Relaxationssignale von wenigstens zwei verschiedenen Zuständen einer Probe so ermittelt, dass ein Unterschiedssignal aus verschiedenen Relaxationssignalen gebildet wird, dass ein zeitlicher Verlauf des Unterschiedssignals eine Funktion U(t) ermittelt wird und dass zwei Zeiten ausgewählt werden, wobei die Auswahl der Zeiten t₁ und t₂ so erfolgt, dass in einem zwischen ihnen liegenden Zeitintervall t_i ein Verhältnis zwischen dem Unterschiedssignal und einem oder mehreren Rauschsignalen größer ist als in anderen Zeitintervallen und dass innerhalb des Zeitintervalls t_i der Wert der Unterschiedsfunktion U(t) ermittelt wird.

Aus der Deutschen Offenlegungsschrift DE 196 35 038 A1 ist ein Verfahren bekannt zur nicht invasiven Messung des zerebralen Blutflusses durch simultane Bestimmung eines dem Blutfluss direkt proportionalen Blutflussindexes durch das Erfassen der zerebralen Anflutungskinetik eines intravenös injizierten Indikatorfarbstoffes mit Absorptionseigenschaften im nah-infraroten Spektrum mittels Nah-Infrarot-Spektroskopie und der arteriellen Anflutungskinetik mittels Pulsdensitometrie und dafür verwendete Anordnung, bestehend aus einem 3-Kanalgerät mit zwei Nah-Infrarot-Spektroskopen, jeweils bestehend aus einem gepulsten Laser mit einer Wellenlänge im nah-infraroten Spektrum, die dem Absorptionsmaximum der eingesetzten Tracersubstanz möglichst entspricht, sowie einem Photomultipler und einem nicht invasiven Messgerät zur pulsdensitrometrischen Bestimmung der Tracerkonzentration im arteriellen Blut des systemischen Kreislaufes.

Aus dem nationalen Teil DE 689 11 583 T2 der Europäischen Patentschrift EP 0 341 783 B1 ist ein Verfahren bekannt zum Bestimmen von Spektralparametern eines auf spektralanalytische Signale in der Zeitebene bezogenen Spektrums in einer Frequenzebene, bei dem die von einer Substanz abgeleiteten spektralanalytischen Signale zum Erhalten von Abtastwerten abgetastet und durch eine komplexe Modellfunktion angenähert werden, die Spektralparameter und exponentiell gedämpfte Sinuskurven enthält, und bei dem außerdem eine anfängliche Schätzung für die nicht-linear in der Modellfunktion auftretenden Spektralparameter auf der Basis der Abtastwerte vorgenommen wird, wobei die Spektralparameter, von der Modellfunktion ausgehend, durch Iteration mit einer Opitimierung nach der Methode der kleinsten Quadrate sehr genau geschätzt werden und Vorkenntnisse in die Modellfunktion eingeführt werden. Das dort dargestellte Verfahren zeichnet sich dadurch aus, dass die Modellfunktion vor der Iteration in Bezug auf die nicht­ linear auftretenden Spektralparameter linearisiert wird und anschließend nach der Methode der kleinsten Quadrate auftretende Skalarprodukt-Terme in eine analytische Form umgewandelt werden, die während der Iteration verwendet wird.

Eine Echtzeit-Korrelationsanalyse von auf Signalanregungen folgenden Messsignalen ist aus Cox, R. W., Jesmanowicz, A., Hyde, J. S. Magnetic Resonance in Medicine, 33, p. 230, 1995 bekannt.

Bei den Meßdaten handelt es sich insbesondere um funktionale Daten wie beispielsweise neuronale Signale, die auch als Feuern von Neuronen bezeichnet werden. Durch neuronale Aktivierung werden Änderungen des cerebralen Blutflusses hervorgerufen. Diese Änderungen können ermittelt werden. Die Transformation des Signals wird durch eine hämodynamische Responsefunktion beschrieben, die eine Zeitverzögerung gegenüber dem Paradigma und eine Verzerrung des Signals berücksichtigen sollte.

Insbesondere in der medizinischen Forschung besteht ein Bedürfnis, Informationen über funktionale Meßdaten, die beispielsweise neuronaler Aktivität entsprechen, zu ermitteln. Funktionale Meßdaten bedeuten hierbei, daß es sich um Daten handelt, die eine durch einen Stimulus induzierte Aktivität oder eine endogene Aktivität wiedergeben. Bekannte Meßmethoden für funktionale, insbesondere neuronale, Aktivität sind Kernresonanztomographie (Nuclear Magnetic Resonance - NMR), funktionale Kernspintomographie (Functional Magnetic Resonance Imaging - fMRI), Elektro-Enzephalographie (EEG) und Magnet-Enzephalographie (MEG). Die Meßsignale werden insbesondere zu einer Bestimmung von funktionaler Aktivität ermittelt, beispielsweise bei Reaktionen von Probanden auf Stimuli.

Durch geeignete Bildgebungsmethoden werden zu einer Detektion der Daten Bildschichten oder Bildvolumina ausgewählt. Eine derartige Bildgebungsmethode kann so verstanden werden, daß sie eine Überführung von Meßdaten in einen Ortsraum beinhaltet. Dies erfolgt beispielsweise in der Kernresonanztomographie durch eine Fourier-Transformation.

Die Bildinformation liegt vorzugsweise in der Form von einem oder mehreren Schichtbildern vor. Ein rekonstruiertes Schichtbild besteht aus Pixeln (= Picture Element = Bildelement), ein Volumendatensatz aus Voxeln (= Volume Element = Volumenelement). Ein Pixel ist ein zweidimensionales Bildelement, beispielsweise ein Quadrat. Aus den einzelnen Pixeln ist das gesamte Bild zusammengesetzt. Ein Voxel ist ein dreidimensionales Volumenelement, beispielsweise ein Quader. Die Abmessungen eines Pixels liegen in der Größenordnung von 1 mm², die eines Voxels in der Größenordnung von 1 mm³. Die Geometrien und Ausdehnungen der Pixel beziehungsweise der Voxel können variabel sein.

Da aus experimentellen Gründen bei Schichtbildern nicht von einer streng zweidimensionalen Ebene ausgegangen werden kann, wird häufig auch hier der Begriff Voxel verwendet, welcher berücksichtigt, daß die Bildebenen eine Eindringtiefe in die dritte Dimension aufweisen.

Durch einen Vergleich der gemessenen Daten, insbesondere eines gemessenen Signalverlaufs in jedem Pixel, mit den Daten einer Modellfunktion, insbesondere deren zeitlichem Verlauf, kann eine für den jeweiligen Stimulus spezifische funktionale Aktivierung detektiert und räumlich lokalisiert werden. Ein Stimulus kann zum Beispiel ein somatosensorischer, akustischer, visueller oder olfaktorischer Reiz sowie eine mentale oder motorische Aufgabe sein. Die Modellfunktion, beziehungsweise die Modell-Zeitreihe, beschreibt eine erwartete Änderung des Signals, insbesondere eines magnetischen Resonanzsignals, infolge funktionaler Aktivierung.

Magnetische Resonanzsignale werden vorzugsweise mittels funktionaler Bildgebung (Functional Magnetic Resonance Imaging - fMRI) aufgenommen. Üblicherweise ist dabei ein kontrasterzeugender Effekt wie der BOLD (Blood Oxygenation Level Dependence) - Effekt der grundlegende biologische Wirkmechanismus. Er beruht darauf, daß nach einer neuronalen Aktivierung nähr- und sauerstoffreicheres Blut (mit mehr Oxyhämoglobin und weniger Deoxyhämoglobin) zu Neuronen transportiert wird. Die neuronale Aktivierung äußert sich in einer Zunahme des Blutflusses in aktivierten Gehirnarealen und somit in einer Abnahme der Deoxyhämoglobinkonzentration im Blut. Da Deoxyhämoglobin (DOH) ein paramagnetischer Stoff ist, welcher die Magnetfeldhomogenität verringert und damit die Signalrelaxation beschleunigt, führt seine Abnahme zu einem Signalanstieg. Anstelle des endogenen Kontrastmittels DOH können auch andere Kontrastmittel treten, die ebenfalls eine Änderung in der Suszeptibilität hervorrufen.

Ein weiteres Auswerteverfahren für fMRI ist in dem Artikel von O. Josephs et. al.: Event-Related MRI, Human Brain Mapping 5, S. 243-248, 1997, dargestellt. Hierbei wird eine Aktivierung mittels vielfacher linearer Regression ermittelt.

Ein weiteres Verfahren zur Ermittlung von neuronaler Aktivität mit fMRI ist in dem Artikel J. C. Rajapakse et. al.: Neuronal and Hemodynamic Responses from Functional MRI Time- Series: A Computational Model, Progress in Connectionist- Based Information Systems (ICONIP'97), N. Kasabov, R. Kozma, K. Ko, R. O'Shea, G. Coghill, T. Gedeon (eds.), Springer- Verlag, Singapur, 1997, S. 30-34, dargestellt.

Die Zeitverzögerung gegenüber dem Paradigma und die Verzerrung des Signals (Verbreiterung des neuronalen Response) rühren von der durch die neuronale Aktivierung hervorgerufenen, langsamen (auf der Zeitskala des Feuerns der Neuronen) Regulation des Blutflusses her. Wenn diese Funktion bekannt ist, besteht die geeignetste Vorgehensweise darin, mit dieser aus dem Zeitverlauf der neuronalen Aktivierung, welche mit dem Zeitverlauf des verwendeten Paradigmas gleichgesetzt werden kann, einen Referenzvektor zu generieren.

In einem Artikel von Geoffrey M. Boynton, Stephen A. Engel, Gary H. Glover und David J. Heeger Linear Systems Analysis of Functional magnetic Resonance Imaging in Human V1:, Journal of Neuroscience, 16 (13): 4207-4221, 1996 wurde offenbart, daß dieser Prozeß sich mathematisch durch eine Faltung beschreiben läßt, da die Blutflußregulation in guter Näherung als lineares System betrachtet werden kann. Die Abweichungen vom linearen Verhalten wurden u. a. in Nonlinear Aspects of the BOLD Response in Functional MRI: Alberto L. Vazquez and Douglas C. Noll; Neuroimage 7: 108-118, 1998 untersucht.

Der Erfindung liegt die Aufgabe zugrunde, ein Verfahren, einen Computer und einen Logikbaustein der jeweils angegebenen Art zu schaffen, die einen schnellen und zuverlässigen Vergleich von Meßdaten mit einer Modellfunktion ermöglichen.

Erfindungsgemäß wird diese Aufgabe dadurch gelöst, daß ein Verfahren zur Erfassung und Auswertung von Meßdaten, wobei die Meßdaten erfaßt und mit einer Modellfunktion verglichen werden, so durchgeführt wird, daß eine Unterschiedsfunktion U zwischen den Meßdaten und der Modellfunktion an einer ersten Stelle (a, b) und an einer zweiten Stelle (a + sda, b + s db) ermittelt wird, daß eine erste partielle Ableitung h₁ der Unterschiedsfunktion U an der ersten Stelle (a, b) ermittelt wird, daß eine zweite partielle Ableitung h₂ der Unterschiedsfunktion U an der ersten Stelle (a, b) ermittelt wird, daß ermittelt wird, ob h₂ < 0 ist und daß bei Bejahung einer Aussage h₂ < 0 die nachfolgenden Schritte durchgeführt werden: ein Quotient h₁/h₂ wird ermittelt, und es wird ermittelt, ob (a + h₁/h₂da) innerhalb eines Intervalls mit Parametergrenzen a_min beziehungsweise a_max liegt und ob (b + h₁/h₂db) innerhalb eines Intervalls mit Parametergrenzen b_min beziehungsweise b_max liegt.

Ferner ist ein Computer zur Erfassung und Auswertung von Meßdaten Gegenstand der Erfindung. Erfindungsgemäß zeichnet sich der Computer zur Durchführung des Verfahrens dadurch aus, daß er mit wenigstens einem Mittel zur Bestimmung einer Unterschiedsfunktion U zwischen den Meßdaten und einer Modellfunktion an einer ersten Stelle (a, b) und an einer zweiten Stelle (a + sda, b + sdb) arbeitet, daß der Computer mit wenigstens einem Mittel zur Bestimmung einer ersten partiellen Ableitung h₁ der Unterschiedsfunktion U an der ersten Stelle (a, b) arbeitet, daß der Computer mit wenigstens einem Mittel zur Bestimmung einer zweiten partiellen Ableitung h₂ der Unterschiedsfunktion U an der ersten Stelle (a, b) arbeitet, daß der Computer mit wenigstens einem Mittel zur Bestimmung eines Quotienten h₁/h₂ arbeitet, daß der Computer mit wenigstens einem Mittel zum Vergleich des Quotienten h₁/h₂ mit einem vorgegebenen Wert arbeitet und daß der Computer in dem Fall, daß h₁/h₂ kleiner ist als der vorgegebene Wert, eine Übereinstimmung zwischen den Meßdaten und der Modellfunktion feststellt.

Der Begriff Computer ist in keiner Weise einschränkend zu verstehen. Es kann sich hierbei um eine beliebige zur Durchführung von Berechnungen geeignete Einheit handeln. Es ist möglich, daß der Computer ein Großrechner, eine Workstation, ein Personalcomputer, ein Mikrocomputer oder eine zur Durchführung von Berechnungen geeignete Schaltung ist. Es ist besonders zweckmäßig, daß die Schaltung in einem Logikbaustein enthalten ist.

Eine bevorzugte Ausführungsform des Computers und des Verfahrens zeichnen sich dadurch aus, daß, falls a + α_jda innerhalb des Intervalls aus den Parametergrenzen a_min beziehungsweise a_max liegt und b + α_jdb innerhalb des Intervalls aus den Parametergrenzen b_min beziehungsweise b_max liegt, überprüft wird, ob der Wert U(a + α_jda, b + α_jdb) an der zweiten Stelle (a + α_jda, b + α_jdb) kleiner ist als der Wert der Unterschiedsfunktion U(a, b) an der ersten Stelle (a, b).

Hierbei ist es zweckmäßig, daß, falls U(a + α_jda, b + α_jdb) an der zweiten Stelle (a + α_jda, b + α_jdb) kleiner ist als der Wert der Unterschiedsfunktion U(a, b) an der ersten Stelle (a, b), der Quotient h₁/h₂ mit einem vorgegebenen Wert ε verglichen wird.

Vorzugsweise wird in dem Fall, daß α_j kleiner ist als der vorgegebene Wert ε, eine Übereinstimmung zwischen den Meßdaten und der Modellfunktion festgestellt.

Außerdem ist es vorteilhaft, daß in dem Fall, daß α_j nicht kleiner ist als der vorgegebene Wert ε, neue Werte a und/oder b aus a + α_jda und/oder b + α_jdb ermittelt werden.

Es ist zweckmäßig, den Computer so zu gestalten beziehungsweise das Verfahren so durchzuführen, daß bei Verneinung der Aussage h₂ < 0 ein Schrittweitenparameter α_j gebildet wird.

Dies geschieht vorzugsweise dadurch, daß der modifizierte Wert von α_j wie folgt gebildet wird: α_j = -λsgn(h₁). Hierbei bezeichnet sgn(h₁) ein Signumfunktion mit einem Wert von +1 für positive Werte, 0 bei 0 und -1 bei negativen Werten.

Es ist vorteilhaft, daß mit dem Schrittweitenparameter α_j ermittelt wird, ob a + α_jda innerhalb eines Intervalls aus Parametergrenzen a_min beziehungsweise a_max liegt und ob b + α_jdb innerhalb eines Intervalls mit Parametergrenzen b_min beziehungsweise b_max liegt.

Falls mit dem Schrittweitenparameter α_j die Summe a + α_jda innerhalb des Intervalls aus den Parametergrenzen a_min beziehungsweise a_max liegt und die Summe b + α_jdb innerhalb des Intervalls mit den Parametergrenzen b_min beziehungsweise b_max liegt, ist es zweckmäßig, daß überprüft wird, ob die Summe der aktuellen Schrittweite α_j und einer früheren Schrittweite kleiner ist als ein vorgegebener Wert δ. Insbesondere wird überprüft, ob α_j + α_j-1 < δ ist. Hierbei bezeichnet j eine Iterationsvariable, die vorzugsweise einer Anzahl von Verfahrensdurchläufen einschließlich des aktuellen Verfahrensdurchlaufes entspricht. Vorzugsweise ist δ klein oder gleich Null.

Vorzugsweise werden, falls die Aussage: α_j + α₁ = < δ mit 1 < j nicht zutrifft, beziehungsweise insbesondere falls die Aussage: α_j + α_j-1 < δ nicht zutrifft, neue Werte a und/oder b aus a + α_jda und/oder b + α_jdb ermittelt, wobei j eine Iterationsvariable bezeichnet, die vorzugsweise höchstens einer tatsächlichen Anzahl von früheren Verfahrensdurchläufen entspricht. Hierdurch wird eine zu häufige Durchführung des Verfahrens vermieden.

Insbesondere sieht die Erfindung vor, daß eine Unterschiedsfunktion U zwischen den Meßdaten und der Modellfunktion an einer ersten Stelle (a, b) und an einer zweiten Stelle (a + sda, b + sdb) ermittelt wird, daß eine erste partielle Ableitung h₁ der Unterschiedsfunktion U an der ersten Stelle (a, b) ermittelt wird, daß eine zweite partielle Ableitung h₂ der Unterschiedsfunktion U an der ersten Stelle (a, b) ermittelt wird, daß ein Quotient α_j = h₁/h₂ ermittelt wird, daß der Quotient α_j mit einem vorgegebenen Wert verglichen wird und daß in dem Fall, daß α_j kleiner ist als der vorgegebene Wert, eine Übereinstimmung zwischen den Meßdaten und der Modellfunktion festgestellt wird.

Eine Anpassung der Meßdaten an die Modellfunktion ist insbesondere dann zweckmäßig, wenn die Werte des Datensatzes ein Rauschen oder eine systematische Veränderung, beispielsweise einen zeitlichen Trend (insbesondere eine Variabilität aufgrund einer Nichtreproduzierbarkeit neuronaler Schwankungen), aufweisen. Ein Beispiel hierfür ist eine Korrelationsanalyse, wie sie unter anderem bei der Detektion neuronaler Aktivität mittels kernmagnetischer Resonanztomographie eingesetzt wird. Sie wird durch die Berechnung der Kreuzkorrelationskoeffizienten zwischen einem Referenzvektor, also einer Modell-Zeitreihe und den Zeitreihen von betrachteten Voxeln (volume element = Volumenelement) realisiert.

Ferner ist es zweckmäßig, eine Unterdrückung von Störsignalen mittels eines Detrending-Verfahrens zu unterstützen. Bei einem Detrending-Verfahren wird der Versuch unternommen, den Effekt von nicht durch einen Stimulus induzierten Signaländerungen in gemessenen Zeitreihen zu reduzieren. Mathematisch heißt das, daß der Meßvektor, also der Vektor, der die gemessene Zeitreihe eines Voxels beinhaltet, in die Summe zweier orthogonaler Vektoren aufgespalten wird. Der Anteil, der durch eine Linear-Kombination von Detrending- Vektoren beschrieben wird, wird verworfen. Die Detrending- Vektoren, die mathematisch gesehen eine Basis bilden, beinhalten die Zeitreihen, aus deren gewichteter Summe der Gesamtstöranteil der gemessenen Zeitreihen beschrieben wird. Durch Anwendung von Detrending auf die in der Korrelationsberechnung eingehenden Meß- und Referenzvektoren kann der Effekt von nicht durch einen Stimulus induzierten Signaländerungen auf die Korrelationsbilder reduziert werden.

Um bessere Korrelationsergebnisse zu erzielen, kann entweder die hämodynamische Responsefunktion oder der Referenzvektor bezüglich bestimmter Parameter, wie der Blutfluß-Latenz-, Anstiegs- und Abfallkonstanten, optimiert werden. Da der Referenzvektor durch Faltung des Paradigmas mit der Responsefunktion hervorgeht, sind beide Methoden mathematisch äquivalent.

Für funktionale Bildgebung in Echtzeit ist diese Optimierung ein wichtiger Bestandteil, da es aufgrund der Geschwindigkeitsanforderungen prinzipiell unmöglich ist, dieselben Daten online mit einem hand-optimierten Vektor erneut auszuwerten. Um regionale Unterschiede in der Gehirnantwort zu berücksichtigen und zu erkunden, sollte der einzige, für alle Bildpunkte identische Referenzvektor zu einem Satz von Referenzvektoren verallgemeinert werden, von denen jeder einem bestimmten Hirnareal zugeordnet wird.

Unabhängig von der verwendeten Methode ist sicherzustellen, daß die Optimierung immer für eine Region mit der "richtigen" Größe erfolgt. Für zu kleine Regionen und eine große Zahl von Freiheitsgraden bei der Optimierung könnte die Referenz- Zeitreihe lediglich die Meß-Zeitreihe replizieren, gleichbedeutend mit einer generellen Erhöhung der Korrelationskoeffizienten. Im entgegengesetzten Fall einer Optimierung für sehr große Regionen ergibt sich ein Referenzvektor, der hauptsächlich jenen globalen Intensitätsschwankungen entspricht, die auf den Herzschlag zurückzuführen sind.

Die einfachste Optimierung kann dadurch erfolgen, daß eine Zeitreihe mit mehreren Referenzvektoren korreliert wird. Diese könnten bis auf unterschiedliche Delayzeiten identisch gewählt werden. Für eine Visualisierung könnte der maximale Korrelationskoeffizient dargestellt werden. Ein wichtiger Vorteil dieses Ansatzes liegt darin, daß er eine iterative Berechnung der Korellationskoeffizienten ermöglicht.

In allen anderen betrachteten Fällen ist dieses vor allen dann vorteilhaft, wenn die Referenzvektoren nur einmal in einer Zeitspanne vollständig ersetzt werden, die lang ist im Vergleich zur Breite des Sliding-Windows. Bei einer Aktualisierung der Referenzvektoren in jedem Zeitschritt werden die Korrelationskoeffizienten am effizientesten direkt gemäß seiner Definitionsgleichung berechnet.

Eine sehr einfache iterative Methode für die Optimierung des Referenzvektors besteht darin - ausgehend von einem initialen Referenzvektor - dadurch iterativ neue Referenzvektoren zu erzeugen, daß der Referenzvektor unter Einbeziehung früherer Meßzeitreihen, insbesondere solcher Meßzeitreihen, deren Korrelationskoeffizienten oberhalb eines bestimmten Schwellwertes liegen, gewählt wird.

Vorzugsweise wird für jede Korrelationsberechnung jeweils der neueste Referenzvektor verwendet. Bei der Anwendung dieser Methode auf reale Meßdaten konnte Aktivierung "zurückgeholt" werden, die ursprünglich durch eine absichtliche Dejustage der Verzögerung der Referenzzeitreihen verloren wurde.

Vorzugsweise wird das Verfahren so durchgeführt, daß eine Berechnung sowohl für aktivierte als auch für nicht aktivierte Regionen erfolgt. Hierdurch wird die Unterscheidung zwischen den aktivierten und den nicht aktivierten Regionen während der Datenauswertung (Fitten) vorgenommen, so daß diese Unterscheidung zuverlässiger ist.

Bevorzugte Ausführungsformen, weitere Vorteile, Besonderheiten und zweckmäßige Weiterbildungen der Erfindung ergeben sich aus der nachfolgenden Darstellung bevorzugter Ausführungsbeispiele der Erfindung anhand von Beispielsrechnungen, Zeichnungen und einer Tabelle.

Von den Zeichnungen zeigt

Fig. 1 eine Meßkurve eines fMRI-Signals im visuellen Kortex mit einer Zeitreihe aus einer 3 × 3 benachbarte Voxel umfassende Region-of-Interest (ROI) im visuellen Kortex, die mit einer wiederholten visuellen Stimulation gemessen wurde,

Fig. 2 ein Blockdiagramm für einen besonders bevorzugten Verfahrensdurchlauf.

Fig. 3 a: Paradigma; b: Referenzvektor, der mit dem Paradigma übereinstimmt; c: wie b. aber mit Verzögerungszeit td; d: Referenzvektor mit exponentiellem Anstieg und Abfall, der noch besser an die Physiologie angepaßt ist; e: aus den Messungen gewonnener Referenzvektor,

Fig. 4 Verlauf eines Korrelationskoeffizienten, der sich mit und ohne Referenzvektor-Optimierung ergibt. Man beachte die durch die Optimierung erzielte Zunahme in der Amplitude und die verbesserte Stabilität,

Fig. 5 Vergleich von zwei verschiedenen hämodynamischen Response-Funktionen mit ähnlicher Dispersion. Die linear-exponentielle Funktion (gepunktete Kurve), die hauptsächlich aus mathematischen Gründen gewählt wurde, hat typischerweise einen schnelleren Signalanstieg und einen langsameren Signalabfall als die Gauß-Funktion (durchgezogene Linie). Diese Asymmetrie der linear-exponentiellen Modellfunktion ist nützlich für die Modellierung von in-vivo fMRI- Daten. Durch Normierung stimmen in diesem Plot beide Funktionen in der Nähe des Maximums überein (a = 3 s, b = 3 s, σ = 3 s), d. h. die Funktionswerte und die zweiten Ableitungen sind am Maximum für beide Funktionen identisch,

Fig. 6 Zeitverlauf der hämodynamischen Abklingzeit, der sich aus der Referenzvektor-Optimierung für die in Fig. 1 gezeigten Daten ergibt. Der Referenzvektor wurde für jeden Zeitpunkt durch Faltung des Paradigmas für die letzten 50 Zeitpunkte mit der in Fig. 4 gezeigten Response-Funktion bestimmt. Die Akkumulationsphase bleibt dadurch unberücksichtigt, d. h. die Berechnung umfaßte nur 450 der 500 Meßzeitpunkte. Die Punkte mit dem kleineren Durchmesser geben die Standardabweichung der hämodynamischen Parameter an.

Fig. 7 Zeitverlauf der hämodynamischen Verzögerungszeit, der sich aus der Referenzvektor-Optimierung für die in Fig. 1 gezeigten Daten ergibt. Der Referenzvektor wurde wiederum für jeden Zeitpunkt durch Faltung des Paradigmas für die letzten 50 Zeitpunkte mit der in Fig. 4 gezeigten Response- Funktion bestimmt. Die Akkumulationsphase bleibt dadurch unberücksichtigt, d. h. die Berechnung umfaßte nur 450 der 500 Meßzeitpunkte. Die Punkte mit dem kleineren Durchmesser geben die Standardabweichung der hämodynamischen Parameter an.

Fig. 8 Überlappung aller 500 Referenzvektoren, von denen jeder 50 Zeitpunkte lang ist, für die in 5 dargestellten Parameterwerte. Der Plot veranschaulicht die gute Stabilität der Optimierungsroutine und offenbart eine leichte Variabilität der Form der Referenzvektoren, die auf eine von der linear-exponentiellen Modellfunktion abweichende Responsefunktion hinweist.

Fig. 9 eine Darstellung einer trapezoiden Referenzfunktion,

Fig. 10 Die obere Kurve stellt den zeitlichen Verlauf des Abklingparameters und die untere den der Verzögerungszeit dar, die für den Fall einer trapezoiden Referenzfunktion erhalten wurde. Diese Kurvenverläufe sind mit denen in Fig. 6 zu vergleichen. Der Unterschied der beiden Referenzfunktionen (trapezoide Funktion und Faltung aus linear-exponentieller Funktion mit Boxcar- Funktion) manifestiert sich hauptsächlich in einem Offset für beide Parameter und stärker verrauschten Werten für den Abklingparameter.

Bei der in Fig. 1 dargestellten Meßkurve handelt es sich um ein fMRI-Signal im visuellen Kortex bei einer Stimulation mit einem Licht oszillierender Intensität. Die Messung erfolgt zweckmäßigerweise mit einem zeitlichen Abstand der Meßpunkte von etwa 1 Sekunde. Die fMRI-Signale werden durch kernspintomographische Untersuchungen des Gehirns von Versuchspersonen ermittelt. In unmittelbarer Nähe des Gesichts der Versuchspersonen wird eine Lichtquelle, insbesondere eine Matrix von Lumineszenzdioden (Light Emitting Diode LED), positioniert und zu Signalblitzen angeregt. Die Anregungsfrequenz liegt vorzugsweise bei etwa 8 Hz. Ein Einwirken der Signalblitze erfolgt über ein mit einem Trägersignal eines Scanners synchronisiertes Zeitintervall von mehreren Sekunden, beispielsweise 5 Sekunden, an das sich ein etwa gleichlanges Ruheintervall anschließt. Bei dem Scanner handelt es sich um einen Vision 1,5 Tesla Ganzkörperscanner der Siemens Medical Systems, Erlangen, mit EPI-Booster und Magnetfeldgradienten von 25 mT/m. Ein derartiger Scanner ist in der Lage, Gradientenfelder innerhalb von etwa 300 µs umzuschalten.

Bei der Bildgebungsmethode handelt es sich vorzugsweise um eine Echo-Planar-Bildgebungsmethode, beispielsweise eine Single-Shot-Echo-Planar-Methode.

Alternativ ist eine spektroskopische Bildgebungsmethode, insbesondere eine wiederholte zweidimensionale Echo- Bildgebungsmethode, einsetzbar. Diese beinhaltet beispielsweise eine wiederholte Anwendung einer zweidimensionalen Echo-Planar-Bildkodierung. Eine räumliche Kodierung erfolgt in einem möglichst kurzen Zeitraum, welcher während eines Signalabfalls mehrfach wiederholt wird und vorzugsweise 20 bis 100 ms beträgt. Durch die mehrfache Wiederholung der Echo-Planar-Kodierung während eines Signalabfalls wird ein Verlauf des Signalabfalls in der Abfolge von rekonstruierten Einzelbildern dargestellt. Eine derartige, gleichfalls vorteilhafte Implementierung der erfindungsgemäßen Methode erfolgt vorzugsweise mit Turbo- PEPSI, wobei PEPSI für Proton-Echo-Planar-Spectroscopic- Imaging steht.

Bei dem Bildgebungsverfahren wird ein Referenzvektor zum Beispiel gemäß dem Ausdruck

verwendet.

Hierbei bezeichnet t die Meßzeit, t_delay die Verzögerungszeit und 1 die zeitliche Länge eines Aktivierungszyklus.

Die Größe m bezeichnet eine Steigung des Signalanstiegs beziehungsweise des Signalabfalls.

Bei den Meßdaten handelt es sich um eine Zeitreihe, die durch eine räumliche Mittelung einer geeigneten Anzahl von vorzugsweise benachbarten Voxeln, insbesondere 3 × 3 Voxeln, gewonnen wird.

Bei einer relativen Verschiebungszeit handelt es sich bespielsweise um die Zeit, die einer Verschiebung des Erwartungswerts t gegenüber dem Erwartungswert des ersten Meßzyklus entspricht.

Eine Berechnung der Verschiebungszeit erfolgt hierbei mit einem Sliding-Window-Verfahren mit N = 50 Werten.

Fig. 2 zeigt ein Blockdiagramm für einen besonders bevorzugten Verfahrensdurchlauf.

Zunächst werden Richtungsvektoren als partielle Ableitungen einer Unterschiedsfunktion gebildet. Anschließend werden die Richtungsvektoren normiert. Die Größen a und b können beispielsweise Zeit- oder Ortsparameter sein. Anschließend wird die Unterschiedsfunktion extrapoliert, wobei Werte der Unterschiedsfunktion an neuen Stellen a + sda, b + sdb ermittelt werden. Von der Unterschiedsfunktion werden an den ursprünglichen Stellen (a, b) eine erste Ableitung h₁ und eine zweite Ableitung h₂ gebildet.

Anschließend wird ermittelt, ob die zweite Ableitung h₂ an der ersten Stelle (a, b) größer als 0 ist. Falls ja, wird anschließend ein Koeffizient α_j = h₁/h₂ gebildet. Die Untersuchung, ob die zweite Ableitung h₂ an der ersten Stelle (a, b) größer als 0 ist, erfolgt, damit in dem Quotienten h₁/h₂ nur solche Werte berücksichtigt werden, die sicherstellen, daß ein Minimum anstelle eines Maximums oder eines Wendepunktes erreicht wird.

Der Koeffizient α_j wird auch als Schrittweitenparameter bezeichnet. Er gibt eine mögliche Entfernung von neuen Werten zu bereits ermittelten Werten an. Anschließend wird ermittelt, ob a + α_jda innerhalb eines Intervalls aus Parametergrenzen a_min beziehungsweise a_max liegt und ob b + α_j db innerhalb eines Intervalls aus Parametergrenzen b_min beziehungsweise b_max liegt.

Falls diese Aussage wahr ist, wird anschließend eine weitere Abfrage vorgenommen. Bei der weiteren Abfrage wird überprüft, ob der Funktionswert U(a + α_jda, b + α_jdb) kleiner ist als der ursprüngliche Funktionswert der Unterschiedsfunktion U(a, b) an der ursprünglichen Stelle. Falls auch diese Abfrage bejaht wird, wird anschließend ermittelt, ob die aktuelle Schrittweite α_j eine vorgegebene minimale Schrittweite ε unterschreitet. Falls diese Abfrage bejaht fragt, wird die Liniensuche beendet und der ermittelte Wert wird gegebenenfalls für eine Liniensuche in einer anderen Richtung verwendet.

Die Unterschiedsfunktion U(a, b) entspricht einer Differenz zwischen gemessenen Daten und einer Modellfunktion. Daher ist es möglich, bei einer Bejahung des eben dargestellten letzten Schrittes von einer maximalen Übereinstimmung zwischen der Modellfunktion und den gemessenen Daten auszugehen.

Falls die letzte Abfrage: α_j kleiner ε verneint wird, werden neue Werte a aus den bisherigen Werten a + α_jda, beziehungsweise ein neuer Wert b = b + α_jdb ermittelt, das heißt, daß folgende Zuweisungen vorgenommen werden: a ← a + α_jda beziehungsweise b ← b + α_jdb. Anschließend wird auch dem Index j ein neuer Wert j + 1 zugewiesen j ← j + 1. Nach dieser Zuweisung erfolgt eine neue Abfrage, wobei ermittelt wird, ob das neue j kleiner ist als ein maximaler Wert des Index j. Falls diese Frage verneint wird, wird die Bearbeitungsroutine abgebrochen. Falls die Abfrage bejaht wird, erfolgt ein neuer Durchlauf der Verfahrensschritte, beginnend mit einer Ermittlung von neuen Zeitableitungen der Unterschiedsfunktion.

Bei diesem neuen Durchlauf wird nach Ermittlung der ersten und der zweiten Zeitableitung überprüft, ob die zweite Zeitableitung h₂ größer als 0 ist. Falls diese Abfrage verneint wird, wird anschließend ein modifizierter Wert von α_j gewählt, α_j = -λsgn(h₁), wobei h₁ die erste Zeitableitung der Unterschiedsfunktion bezeichnet. Beispielsweise ist λ = 0,5.

Hiernach erfolgt eine Abfrage, ob a + α_jda innerhalb eines Intervalls aus einem minimalen Wert von a_min und a_max liegt und ob b + α_jdb innerhalb eines durch b_min und b_max begrenzten Intervalls liegt. Diese Abfrage erfolgt, damit nur solche Werte von a + α_jda beziehungsweise b + α_jdb in das weitere Verfahren eingehen, die innerhalb eines geeigneten Intervalls liegen. Hierdurch wird vermieden, daß physikalisch unrealistische Parameter im weiteren Verfahren berücksichtigt werden. Somit ist es möglich, physikalische Randbedingungen zu berücksichtigen. Falls die Abfrage verneint wird, wird das Verfahren abgebrochen.

Falls die Abfrage bejaht wird, erfolgt anschließend eine neue Abfrage, bei der ermittelt wird, ob α_j + α_j-1 ≅ 0 ist. Falls diese Abfrage bejaht wird, wird das Verfahren abgebrochen. Falls diese Abfrage verneint wird, werden anschließend auf die weiter oben beschriebene Weise neue Werte für a und b eingesetzt und es erfolgt ein Übergang von j zu j + 1.

Der Übergang zu neuen Schrittweiten ist lediglich beispielhaft zu verstehen. Insbesondere ist es möglich, die Schrittweiten auf eine andere Weise zu modifizieren, beispielsweise dadurch, daß die Schrittweite zwischen verschiedenen Verfahrensdurchläufen jeweils verringert wird.

Hierdurch wird eine erzwungene Konvergenz erzielt und es werden Lösungen gewonnen, die jeweils gegenüber einer Lösung in einem vorhergehenden Verfahrensdurchlauf verbessert sind. Hierdurch ist das Verfahren besonders robust und für Praxiseinsätze, beispielsweise für eine Echtzeitverarbeitung der Meßdaten, besonders geeignet.

Konventionelle Gradientenabstiegsmethoden eignen sich vorzugsweise für Daten, die wenig verrauscht sind. Bei tatsächlichen Meßdaten, die bedingt durch äußere Einflüsse, eine Streuung aufweisen, eignen sich die bekannten Verfahren nur begrenzt.

Durch das erfindungsgemäße Verfahren werden einer Optimierung zugrundeliegende Schrittweiten an einen tatsächlichen Verlauf der Daten angepaßt. Daher kann das Verfahren weitgehend von Ausgangsparametern unabhängig durchgeführt werden.

Ein weiterer Vorteil der Erfindung ist, daß geeignete Parameter bereits nach wenigen Verfahrensdurchläufen ermittelt werden. Hierdurch ist das Verfahren besonders schnell.

Es hat sich gezeigt, daß durch das erfindungsgemäße Verfahren mittels einer Parameterbegrenzung unphysikalisch große Sprungbreiten vermieden werden.

Durch die Wahl von anfänglich großen Schrittweiten wird eine Wahrscheinlichkeit für ein Festsetzen in lokalen Minima verringert.

Interpolierte Zeitreihen sind in Fig. 3 in Teilbildern (a. bis e.) symbolisch dargestellt. Hierbei bezeichnet a: Werte für das Paradigma; b: Werte für einen Referenzvektor, der mit dem Paradigma übereinstimmt; c: wie b. aber mit Verzögerungszeit td; d: Werte für einen Referenzvektor mit exponentiellem Anstieg und Abfall, der noch besser an die Physiologie angepaßt ist. Die Zeitkonstante für den Anstieg (1 - 1/e des Maximalwertes) ist tri, die für den exponentiellen Abfall (1/e des Maximalwertes) ist tf; e: Werte für einen aus der Messungen gewonnenen Referenzvektor.

Anhand der in Fig. 3 dargestellten Verschiebung des Erwartungswerts t wird festgestellt, daß Veränderungen des Erwartungswerts t mit einer Genauigkeit detektiert werden können, die höher ist als es der physikalischen Zeitauflösung entspricht. Dieses überraschende Ergebnis, das sich bei einem möglichst hohen Signal-zu-Rausch-Verhältnis der gemessenen Daten erzielen läßt, kann wie folgt erklärt werden: Bereits nach einem kurzen Zeitraum tritt eine signifikante, auswertbare Änderung der Signalamplitude auf, aus der sich auf die Zeitverzögerung rückschließen läßt.

In Fig. 4 ist eine Abhängigkeit des Korrelationskoeffizienten von der gemessenen Zeit mit einer konventionellen Optimierung und einer erfindungsgemäßen Sliding-Window-Optimierung in Abhängigkeit von der Meßzeit dargestellt.

Bei der eingesetzten konventionellen Optimierung handelt es sich um den Einsatz einer Modellfunktion mit konstanten, jedoch an die gesamte Menge der Daten angepaßten Parametern. Die Anpassung erfolgte so, daß im Mittel die Korrelationskoeffizienten einen möglichst hohen Wert erhielten. Somit stellt die hier dargestellte Kurve ein Optimum für die bekannten Auswerteverfahren dar.

Erfindungsgemäß wird vorzugsweise ein Sliding-Window- Verfahren eingesetzt, bei dem im Anschluß an eine anfängliche Auffüllphase von Daten eine Phase eines "stationären Gleichgewichts" eintritt, das heißt die Anzahl der Datenwerte bleibt konstant, wobei mit fortschreitender Messung der jeweils älteste Datenwert entfällt und der jüngste Datenwert hinzugefügt wird. Das Sliding-Window-Verfahren wurde so durchgeführt, daß eine Überschreibung mit dem neuen Wert so erfolgt, daß er jeweils einen Wert mit gleicher Phasenlage überschreibt.

Nachfolgend wird eine bevorzugte Variante des Verfahrens mit einer Optimierung parameterisierter Referenzvektoren dargestellt:

Der Referenzvektor hängt von einem Satz von Parametern {ω_j} ab, wie z. B. dem Zeitpunkt für den Beginn und das Ende einer Stimulus-Ein-Phase oder den Zeiten für die allmähliche Zu- oder Abnahme des Blutflusses. Dann wird die folgende quadratische Form bezüglich der Parameter ω_j optimiert, wofür die Gradientenabstiegsmethode oder die (in der Nähe eines Minimums) schneller konvergierende Methode der konjugierten Gradienten in Zusammenhang mit dem anhand von Fig. 2 erläuterten Verfahren verwendet werden kann

E = |_S - α_S(ω_j)|², (1)

wobei sich das Minimum von E für den sogenannten Überlappwert

α = (_{S S})/_S ² = ρ|_S|/|_S| (2)

ergibt. Diese Beziehung erlaubt die Berechnung von α. Für linear-unabhängige Detrending-Vektoren ist die quadratische Form Eq. (1) äquivalent zu

E = |P_S - αP_S(ω_j)|², (3)

mit P_S = I - S(S^TS)^-1S^T, wobei I eine Einheitsmatrix ist.

Daraus folgt

Mit Gleichung

für den Überlappwert α kann nun eine Minimierung zum Beispiel durch Änderung der Parameter um

erfolgen, wofür die Gradientenabstiegsmethode oder die (in der Nähe eines Minimums) schneller konvergierende Methode der konjugierten Gradienten in Zusammenhang mit dem anhand von Fig. 2 erläuterten Verfahren verwendet werden kann. Hierbei sollte der Schrittweiten-Parameter ε, der auch als Lernrate bezeichnet wird, hinreichend klein gewählt werden. Wie bereits erwähnt wurde, ist es vorteilhaft, den Referenzvektor für Regionen, nicht aber für einzelne Voxel zu optimieren. Das läßt sich dadurch erreichen, daß Gleichung (1) um eine Summe über alle Meßvektoren einer Region-of-Interest erweitert wird.

Es soll nun das Optimierungsverfahren an einem Beispiel genauer beschrieben werden. Das Detrending erfolgt nur mit einem Detrending-Vektor, bei dem alle Einträge auf 1 gesetzt werden; dieses entspricht einer Korrektur in Form einer Basislinienverschiebung.

Damit ein Optimierungs-Algorithmus vollkommen automatisch, d. h. ohne Eingriff durch den Benutzer arbeiten kann, muß er besonderen Anforderungen an die numerische Stabilität genügen. Für eine Berechnung in Echtzeit ist ferner eine schnelle Konvergenz erforderlich. Daher sollte die hämodynamische Responsefunktion zusätzliche Bedingungen erfüllen. Für eine experimentelle Untersuchung zur hämodynamischen Response siehe: "The variability of human BOLD hemodynamic responses: Aguirre; NeuroImage, 1998, Vol. 8(4), p. 360-369, den in der Literatur "J. Rajapakse, F. Kruggel, D. Y. von Cramon, Neuronal and hemodynamic responses from functional MRI time-series: A comutational model, in "Progress in Connectionist-Based Information Systems (ICONIP'97)" (N. Kasabov, R. Kozma, K. Ko, R. O'Shea, G. Coghill, T. Gedeon, Eds.), p. 30-34, Springer, Singapur, 1997" und "Modeling Hemodynamic Response for Analysis of Functional MRI Time-Series: Jagath c. Rajapakse, Frithjof Kruggel, Jose M. Maisog und D. Yves von Cramon; Human Brain Mapping 6: 283-300, 1998" mit vorgeschlagenen Gauß- und Poisson-Funktionen Die Antwortfunktion soll eine Referenzfunktion ergeben, deren zweite Ableitungen existieren und durch analytische Ausdrücke gegeben sind. Eine geeignete Funktion ist

(t) = (t - b)e^-a(t-b)θ(t - b) (7)

(wobei θ(x), die Heaviside θ-Funktion, definiert ist als θ(x < 0) = 0, θ(x < 0) = 1 und θ(x = 0) = 1/2), die im folgenden als linear-exponentiell bezeichnet werden soll. Diese ist in Fig. 5 zusammen mit einer Gauß-Funktion dargestellt.

Die Betrachtung gemessener Zeitverläufe ergibt, daß diese asymmetrische Funktion eine geeignete Wahl für die hämodynamische Responsefunktion ist. In Fig. 5 sind beide Funktionen so normiert, daß sie denselben Maximalwert und am Maximum dieselben zweiten Ableitungen haben. Nach dem Maximum fällt die Gauß-Funktion wesentlich schneller ab als die linear-exponentielle Funktion.

Die Funktion h(t) hat nur zwei freie Parameter. Die Zahl der Parameter muß hinreichend klein sein, damit auch bei einer kurzen Fensterlänge der Referenzvektor nicht zu stark vom Rauschen beeinflußt wird. Damit wird insbesondere die Anpassung an Zeitreihen vermieden, die keine stimuluskorrelierte Aktivierung zeigen. Der Parameter a gibt eine inverse Abfallzeit und b eine zusätzliche Verzögerungszeit an. Beide Parameter zusammen geben die Lage des Maximums von h(t) an, welches bei b + 1/a liegt. Der Verlauf der neuronalen Aktivierung, von dem angenommen wird, daß er sich zum experimentellen Paradigma synchron verhält, wird durch eine Box-car-Funktion beschrieben, deren steigende Flanke bei t_c und deren fallende Flanke bei t_d liegt. Die Referenzfunktion ist gegeben durch eine Faltung der Box-car- Funktion mit h(t):

Hierbei wurde die θ-Funktion in den Integralgrenzen verarbeitet, wodurch die Funktion

(t) = (t - b)e^-a(t-b) (9)

anstelle von h(t) auftritt. Die analytische Lösung für r(t):

wird in Gleichung (1) eingesetzt. Im Prinzip liegt ein 2D- Optimierungsproblem in den Parametern a und b vor. Die analytische Berechnung der partiellen Ableitungen kann aber wesentlich vereinfacht werden, wenn das 2D- künstlich zu einem 3D-Problem erweitert wird, indem α, das von a und b abhängt, als zusätzlicher Optimierungsparameter behandelt wird. Genau am Minimum der durch Gleichung (1) gegebenen Funktion entspricht α dem mit Gleichung (2) gefundenen Wert.

Die Referenzfunktion r(t) und ihre ersten und zweiten Ableitungen bezüglich a, b und α wurden mit dem Computer- Algebra-System MAPLE (Waterloo Maple Inc., Kanada) ermittelt, mit dem auch die analytischen Ausdrücke in optimierten C- Sourcecode übersetzt wurden. Die Ableitungen können für eine besonders vorteilhafte Implementierung des erfindungsgemäßen Verfahrens verwendet werden.

Die beschriebene Referenzvektor-Optimierung wurde zur Auswertung einer Einzel-Zeitreihe verwendet. Die Einzel- Zeitreihe beinhaltet einen Datensatz, der wie folgt gewonnen wurde: Es erfolgte eine visuelle Stimulation mit 8 Hz Flickerlicht (effektives T_R: 1 s, T_E: 67 ms, Anzahl der Messungen: 512, Matrixgröße: 32 × 32, Voxelgröße 6 × 6 × 3 mm³, 12 Baseline-Messungen zu Beginn, 10 Stimulations-Zyklen, bei der die Stimulations-Bedingung 10 Messungen und die Kontrollbedingung 40 Messungen umfaßte).

Das Timing der visuellen Stimulation war mit dem Pulssequenz- Trigger des Scanners synchronisiert. Für den Worstcase-Fall, nämlich daß das Parameter-Feedback von einer Sliding-Window- Position zur nächsten innerhalb der Optimierungsroutine deaktiviert wird, beträgt die Rechenzeit für ein einzelnes Voxel und für eine Position des Sliding-Windows 0,2 s auf einer SUN Ultra Sparc 1 (143 MHz). Wenn das Feedback aktiviert ist, erhöht sich die Rechengeschwindigkeit um mindestens eine Größenordnung.

Vor der Implementierung der 3D- wurde zunächst eine 4D- Optimierung implementiert und getestet. Diese benötigte 50% mehr Rechenleistung und ergab einen wesentlich verrauschteren Zeitverlauf für den Abkling-Parameter. Die Eliminierung des verbliebenen "künstlichen" Optimierungsparameters α in der 3D-Optimierung durch Einsetzen von Gleichung (2) in die Kostenfunktion könnte zu einer weiteren Geschwindigkeitssteigerung der numerischen Optimierung führen. Das setzt voraus, daß die Anzahl der erforderlichen Iterationsschritte stärker ab - als die Komplexität der analytischen Ausdrücke für die ersten und zweiten Ableitungen zunimmt.

Wenn vollständige Datensätze für eine Schicht oder ein Volumen verarbeitet werden, sind die angegebenen Berechnungszeiten mit der Anzahl der Kompartimente, für die die Optimierung erfolgen soll, zu multiplizieren. Fig. 4 zeigt eine Zeitreihe, die durch Mittelung über 3 × 3 Voxel im visuellen Kortex erzeugt wurde. In Fig. 6 sind die aus ihr mittels iterativer Optimierung berechneten Zeitverläufe für die Abklingzeit 1/a dargestellt. In Fig. 7 sind die aus ihr gleichfalls mittels iterativer Optimierung berechneten Zeitverläufe für die Verzögerungszeit b + 1/a dargestellt. Die statistische Signifikanz der ermittelten Parameteränderungen wurde mit einer Monte-Carlo-Simulation untersucht, indem das Fitting für die in Fig. 4 gezeigten Daten 100mal wiederholt wurde, wobei in jedem Durchlauf gaußverteilte Zufallszahlen (σ = 25) hinzuaddiert wurden. Die Standardabweichung der Verzögerungszeit, 0,5 s, und die der Abklingzeit, 0,3 s, zeigen, daß die berechneten Änderungen der Abklingparameter statistisch signifikanter sind als die der Verzögerungszeit.

Diese bei einer exakten Wiederholung einer visuellen Stimulation auftretenden Schwankungen der BOLD-Response sind überraschend, da angenommen werden könnte, daß die auf der niedrigsten Ebene der visuellen Verarbeitung V1 ablaufenden Prozesse reproduzierbar sind.

In Fig. 4 sind die Korrelationskoeffizienten, die sich mit und ohne Referenzvektor-Optimierung ergeben, in Abhängigkeit von der Sliding-Window-Position dargestellt. Der Mittelwert und die Standardabweichung des Korrelationskoeffizienten sind 0,93 ± 0,01 beziehungsweise 0,75 ± 0,03. In dem Fall, wo keine Optimierung erfolgt, wird eine Boxcar-Funktion mit einem Delay von 4 s verwendet. Von allen ganzzahligen Verzögerungszeiten (T_R = 1 s) kleiner als 8 s ergab diese Verzögerungszeit den größten Korrelationskoeffizienten. Die Zunahme des Korrelationskoeffizienten ist kein globaler Effekt, sondern erfolgt aktivierungsspezifisch, wie die Auswertung einer Zeitreihe ohne stimulus-induzierte Signalanteile gezeigt hat. Die verwendete fMRI-Zeitreihe, die aus einer Basislinien-Messung, also einer Messung ohne Stimulus-Präsentation, stammt, enthält lediglich "physiologisches Rauschen". Es wurde vom Anmelder gezeigt, daß der Korrelationskoeffizient für die Rausch-Zeitreihe und einem fiktiven Paradigma unverändert bleibt: Der Mittelwert und die Standardabweichung für die Korrelationskoeffizienten, die sich beide nach iterativer Optimierung des Referenzvektors ergaben, betrug bei der gewählten, als repräsentativ anzusehenden Basislinien-Zeitreihe 0,064 ± 0,125.

Bei fehlender Aktivierung verlangsamt sich das Fitting sehr drastisch - um einen Faktor 6 bis 7 für die untersuchte Zeitreihe - weil dann kein ausgeprägtes globales Minimum existiert. Es ist daher zweckmäßig, der Verarbeitung von 2D- oder 3D-Datensätzen einen einfachen Test vorzuschalten, um eine Referenzvektor-Optimierung für Voxel zu vermeiden, die in nicht-aktivierten Arealen liegen. Mit einem unabhängigen Test könnte in vielen Fällen auch das Auftreten von lokalen Minima detektiert werden. Bereits ein oder zwei Gradienten- Schritte genügen einer Routine für eine Berechnung von konjugierten Gradienten, um den Korrelationskoeffizienten auf einen Wert zu erhöhen, der nahe am Endergebnis liegt. Weitere Iterationen werden jedoch benötigt, damit sich die errechneten Parameterwerte stabilisieren.

Fig. 8 gibt eine Überlagerung von allen optimierten Referenzvektoren an, die für die verschiedenen Zeitpunkte berechnet wurden. Diese Überlagerung, welche als eine Art "gefilterte" Gehirn-Antwortfunktion aufgefaßt werden kann, zeigt eine gute Kohärenz der Referenzvektoren für aufeinanderfolgende Sliding-Window-Positionen. Dieser Plot demonstriert die Stabilität des implementierten Optimierungs- Algorithmus. Um die Modellabhängigkeit der berechneten Parameteränderungen zu bestimmen, wurden dieselben Berechnungen mit einer in Fig. 9 dargestellten, aus zwei tanh-Funktionen bestehenden, trapezoiden Referenzfunktion durchgeführt. Der Zeitverlauf der Parameteränderungen konnte - wie in Fig. 10 dargestellt - im wesentlichen reproduziert werden.

Für verrauschte Aktivierungsdaten (Signal-zu-Rausch- Verhältnis = 1 : 1) ist ein Fitting immer noch möglich; allerdings wiesen die berechneten Parameter dann größere Fluktuationen auf.

Zu Zeitpunkten, zu denen ein zweites Aktivierungssignal im Sliding-Window erscheint, sind kleinere Instabilitäten des Verzögerns- und des Abkling-Parameters zu beobachten, deren Ursache im wesentlichen in der unterschiedlichen Amplitude aufeinanderfolgender Aktivierungssignale liegt (Fig. 6 und Fig. 7).

Das erfindungsgemäße Verfahren ist nicht auf zwei Parameter beschränkt, sondern kann auch kann mit einer größeren Anzahl von Parametern durchgeführt werden.

Claims (20)

1. Verfahren zur Erfassung und Auswertung von Messdaten, wobei die Messdaten erfasst und mit einer Modellfunktion verglichen werden, dadurch gekennzeichnet,

• - dass eine Unterschiedsfunktion U zwischen den Messdaten und der Modellfunktion an einer ersten Stelle (a, b) und an einer zweiten Stelle (a + sda, b + sdb) ermittelt wird,
• - dass eine erste partielle Ableitung h₁ der Unterschiedsfunktion U an der ersten Stelle (a, b) ermittelt wird,
• - dass eine zweite partielle Ableitung h₂ der Unterschiedsfunktion U an der ersten Stelle (a, b) ermittelt wird,
• - dass ermittelt wird, ob h₂ < 0 ist und dass bei Bejahung einer Aussage h₂ < 0 die nachfolgenden Schritte durchgeführt werden:
• - ein Quotient h₁/h₂ wird ermittelt und
• - es wird ermittelt, ob (a + h₁/h₂da) innerhalb eines Intervalls mit Parametergrenzen a_min beziehungsweise a_max liegt und ob (b + h₁/h₂db) innerhalb eines Intervalls mit Parametergrenzen b_min beziehungsweise b_max liegt.

2. Verfahren nach Anspruch 1, dadurch gekennzeichnet,

• - dass, falls (a + h₁/h₂da) innerhalb des Intervalls mit den Parametergrenzen a_min beziehungsweise a_max liegt und (b + h₁/h₂db) innerhalb des Intervalls mit den Parametergrenzen b_min beziehungsweise b_max liegt, überprüft wird, ob der Wert der Unterschiedsfunktion U(a + h₁/h₂da, b + h₁/h₂db) an der zweiten Stelle (a + h₁/h₂da, b + h₁/h₂db) kleiner ist als der Wert der Unterschiedsfunktion U(a, b) an der ersten Stelle (a, b).

3. Verfahren nach Anspruch 2, dadurch gekennzeichnet,

• - dass, falls der Wert der Unterschiedsfunktion U(a + h₁/h₂da, b + h₁/h₂db) an der zweiten Stelle (a + h₁/h₂ da, b + h₁/h₂db) kleiner ist als der Wert der Unterschiedsfunktion U(a, b) an der ersten Stelle (a, b), der Quotient h₁/h₂ mit einem vorgegebenen Wert ε verglichen wird.

4. Verfahren nach Anspruch 3, dadurch gekennzeichnet,

• - dass in dem Fall, daß h₁/h₂ kleiner ist als der vorgegebene Wert ε, eine Übereinstimmung zwischen den Messdaten und der Modellfunktion festgestellt wird.

5. Verfahren nach Anspruch 4, dadurch gekennzeichnet,

• - dass in dem Fall, dass h₁/h₂ nicht kleiner ist als der vorgegebene Wert ε, neue Werte a und/oder b aus (a + h₁/h₂da) und/oder (b + h₁/h₂db) ermittelt werden.

6. Verfahren nach einem oder mehreren der Ansprüche 1 bis 5, dadurch gekenn­ zeichnet,

• - dass bei Verneinung der Aussage h₂ < 0 ein modifizierter Wert α_j gebildet wird.

7. Verfahren nach Anspruch 6, dadurch gekennzeichnet,

• - dass der modifizierte Wert α_j wie folgt gebildet wird:
  α_j = -λsgn(h₁).

8. Verfahren nach einem der Ansprüche 6 oder 7, dadurch gekennzeich­ net,

• - dass mit dem modifizierten Wert α_j ermittelt wird, ob (a + α_jda) innerhalb eines Intervalls mit den Parametergrenzen a_min beziehungsweise a_max liegt und ob, (b + α_jdb) innerhalb eines Intervalls mit den Parametergrenzen b_min beziehungsweise b_max liegt.

9. Verfahren nach Anspruch 8, dadurch gekennzeichnet,

• - dass, falls mit dem modifizierten Wert α_j die Summe (a + α_jda) innerhalb des Intervalls mit den Parametergrenzen a_min beziehungsweise a_max liegt und die Summe (b + α_jdb) innerhalb des Intervalls mit den Parametergrenzen b_min beziehungsweise b_max liegt,
• - überprüft wird, ob eine Summe aus der aktuellen Schrittweite α_j und der früheren Schrittweite kleiner ist als ein vorgegebener Wert 6.

10. Verfahren nach Anspruch 9, dadurch gekennzeichnet,

• - dass, falls die Aussage, dass die Summe aus der aktuellen Schrittweite α_j und der früheren Schrittweite kleiner ist als ein vorgegebener Wert δ, nicht zutrifft, neue Werte a und/oder b aus a + α_jda und/oder b + α_jdb ermittelt werden.

11. Verfahren nach Anspruch 10, dadurch gekennzeichnet, dass mit fortschreitender Aufnahme von Messdaten ältere Datenwerte entfallen und durch neuere Datenwerte ersetzt werden.

12. Verfahren nach Anspruch 11, dadurch gekennzeichnet, dass eine Überschreibung mit den neuen Datenwerten so erfolgt, dass die neuen Datenwerte jeweils einen alten Datenwert mit gleicher Phasenlage überschreiben.

13. Verfahren nach einem oder mehreren der Ansprüche 1 bis 12, dadurch gekenn­ zeichnet, dass ein Referenzvektor gebildet wird.

14. Verfahren nach Anspruch 13, dadurch gekennzeichnet, dass der Referenzvektor parameterisiert wird.

15. Verfahren nach einem der Ansprüche 13 oder 14, dadurch gekenn­ zeichnet, dass der Referenzvektor durch eine Referenzfunktion r(t) wiedergegeben wird.

16. Verfahren nach Anspruch 15, dadurch gekennzeichnet, dass die Referenzfunktion r(t) lautet:

17. Verfahren nach einem oder mehreren der Ansprüche 1 bis 16, dadurch gekenn­ zeichnet, dass eine Berechnung sowohl für aktivierte als auch für nicht aktivierte Regionen erfolgt.

18. Computer zur Durchführung des Verfahrens nach einem der Ansprüche 1 bis 17, dadurch gekennzeichnet,

• - dass der Computer mit wenigstens einem Mittel zur Bestimmung einer Unterschiedsfunktion U zwischen den Messdaten und einer Modellfunktion an einer ersten Stelle (a, b) und an einer zweiten Stelle (a + sda, b + sdb) arbeitet,
• - dass der Computer mit wenigstens einem Mittel zur Bestimmung einer ersten partiellen Ableitung h₁ der Unterschiedsfunktion U an der ersten Stelle (a, b) arbeitet,
• - dass der Computer mit wenigstens einem Mittel zur Bestimmung einer zweiten partiellen Ableitung h₂ der Unterschiedsfunktion U an der ersten Stelle (a, b) arbeitet,
• - dass der Computer mit wenigstens einem Mittel zur Bestimmung eines Quotienten h₁/h₂ arbeitet,
• - dass der Computer mit wenigstens einem Mittel zum Vergleich des Quotienten h₁/h₂ mit einem vorgegebenen Wert arbeitet.

19. Computer nach Anspruch 18, dadurch gekennzeichnet,

• - dass der Computer in dem Fall, daß h₁/h₂ kleiner ist als der vorgegebene Wert, eine Übereinstimmung zwischen den Messdaten und der Modellfunktion feststellt.

20. Logikbaustein zur Durchführung des Verfahrens nach einem der Ansprüche 1 bis 17, dadurch gekennzeichnet,

• - dass der Logikbaustein mit wenigstens einem Mittel zur Bestimmung einer Unterschiedsfunktion U zwischen den Messdaten und einer Modellfunktion an einer ersten Stelle (a, b) und an einer zweiten Stelle (a + sda, b + sdb) arbeitet,
• - dass der Logikbaustein mit wenigstens einem Mittel zur Bestimmung einer ersten partiellen Ableitung h₁ der Unterschiedsfunktion U an der ersten Stelle (a, b) arbeitet,
• - dass der Logikbaustein mit wenigstens einem Mittel zur Bestimmung einer zweiten partiellen Ableitung h₂ der Unterschiedsfunktion U an der ersten Stelle (a, b) arbeitet,
• - dass der Logikbaustein mit wenigstens einem Mittel zur Bestimmung eines Quotienten h₁/h₂ arbeitet und
• - dass der Logikbaustein mit wenigstens einem Mittel zum Vergleich des Quotienten h₁/h₂ mit einem vorgegebenen Wert arbeitet.