Modellierung tumorinduzierter Gewebedeformation als Optimierungsproblem mit weicher Nebenbedingung

Abstract

Kurzfassung. Ein ungelo¨stes Problem in der nicht-rigiden Bildregistrierung ist die Behandlung von pathologie-bedingten, morphologischen Unterschieden, wie sie beispielsweise bei der ra¨umlichen Normalisierung neuroradiologischer Datensa¨tze, die Tumorpathologie abbilden, auftreten. Mit der vorliegenden Arbeit liefern wir einen Baustein fu¨r einen modellbasierten Lo¨sungsansatz. Wir schlagen vor, die entstehenden Irregularita¨ten durch eine explizite Modellierung der Pathologie zu umgehen. Im Detail stellen wir die Erweiterung eines auf der Formulierung eines Optimierungsproblems basierenden Ansatzes zur Modellierung tumorinduzierter Gewebedeformation vor. Dieser bietet potentiell die Mo¨glichkeit einer direkten Integration in Verfahren der nicht-rigiden Bildregistrierung. Neben einer Darstellung des theoretischen Zusammenhangs mit existierenden Verfahren, zeigen wir experimentell, dass die Hinzunahme einer nicht-linearen Wichtung der Terme des Zielfunktionals eine adaptivere Steuerung der resultierenden Deformationsmuster erlaubt.

-

Die nicht-rigide Bildregistrierung [ 1, 2 ] ist ein ma¨chtiges und etabliertes Werkzeug zur Analyse medizinischer Bilddaten. Es gilt, eine ra¨umliche Abbildung y zwischen zwei oder mehreren Bildern zu finden, die korrespondierende Strukturen zur u¨berlagerung bringt. Aus mathematischer Sicht handelt es sich hierbei um ein schlecht gestelltes, inverses Problem. Um den Lo¨sungsraum auf physikalisch plausible Abbildungen zu beschra¨nken und damit einhergehend, Forderungen an die Stetigkeit der Funktion y zu stellen, wird das Zielfunktional typischerweise um einen Regularisierungsterm (Gla¨tter) erweitert [ 1 ].

Betrachten wir das Problem der nicht-rigiden Registrierung eines neuroradiologischen Datensatzes, der Tumorpathologie abbildet, mit einem standardisierten Referenzatlas (Anwendungsbeispiele sind die atlasbasierte Segmentierung oder die Generierung statistischer Pathologieatlanten), so sind die genannten Grundvoraussetzungen bzgl. Korrespondenz und Stetigkeit verletzt: Es existiert weder eine der Pathologie entsprechende, morphologische Struktur im Atlas noch ist die durch den Tumor induzierte Gewebedeformation zwangsweise ”glatt“. Aufgrund dieser Irregularita¨ten sind existierende Registrierverfahren [ 1, 2 ] per Definition nicht in der Lage, eine geeignete Abbildung y zwischen den gegebenen Bilddaten zu bestimmen. Es ergibt sich die Forderung nach neuen Lo¨sungsstrategien. Eine naheliegende Vorgehensweise ist, den durch die Pathologie betroffenen Bereich aus der Optimierung des Zielfunktionals auszublenden [ 3 ]. Eine ku¨rzlich vorgeschlagene [ 4 ], elaboriertere Strategie – die nicht nur eine Lo¨sung fu¨r das Korrespondenzproblem bietet, sondern auch die Mo¨glichkeit der Scha¨tzung patientenindividueller Modellparameter erlaubt – ist die explizite Integration eines bio-physikalischen Pathologiemodells [ 5, 6 ] in die Registrierung. Wir folgen diesem Ansatz und stellen eine Erweiterung des in [ 6 ] vorgestellten Modells vor.

Der wesentliche Beitrag der vorliegenden Arbeit ist (i) die Einfu¨hrung einer nicht-linearen Wichtung der Terme des Zielfunktionals in Abha¨ngigkeit der vorliegenden Zelldichte, (ii) eine explizite Bestrafung der Deformation rigider Strukturen, (iii) eine adaptive Auswertung des Zielfunktionals durch eine numerische Integration u¨ber das Teilgebiet des Einflussbereiches der einzelnen Parameter des Deformationsmodells und (iv) die Verwendung analytischer Ableitung fu¨r die Optimierung des Zielfunktionals. Weiter zeigen wir mittels Variationsrechnung den theoretischen Zusammenhang mit existierenden Arbeiten [ 4 ]. 2

Material und Methoden

In den folgenden Abschnitten rekapitulieren wir den Ansatz zur Modellierung der raumzeitlichen Dynamik kanzero¨ser Zellen und berichten u¨ber die theoretischen und numerischen Hintergru¨nde des vorgeschlagenen Modells. 2.1

Mathematisches Modell

Wir beginnen mit der Beschreibung des kontinuierlichen Modells fu¨r die raumzeitliche Dynamik der Tumorzelldichte. Gesucht wird eine Funktion u˜ : ΩB × [0, τ ] → R 0, ΩB ⊂ Ω := (ω11, ω12) × · · · × (ωd1, ωd2) ⊂ Rd, d ∈ {2, 3}, welche die partielle Differentialgleichung ∂tu˜(x, t) = ∇ · (D(x)∇u˜(x, t)) + ρu˜(x, t)(1 − u˜(x, t)), auf ΩB × (0, τ ] (1) mit den Rand- und Anfangsbedingungen ∂n^u˜(x, t) = 0 auf ∂ΩB × (0, τ ] und u˜(x, 0) = u˜S exp(−∥x − xS ∥22/(2σ2)) auf Ω¯B (unter der Annahme u˜ ∈ C2(ΩB × (0, τ ), R 0) ∩ C(Ω¯B × [0, τ ])) erfu¨llt. Hierbei repra¨sentiert ΩB mit Abschluss ¯ ΩB := ΩB∪∂ΩB das durch das zerebrale Gewebe eingenommene Gebiet, ρ ∈ R 0 die Wachstumsrate, D : Ω → Rd d einen linear gewichteten Diffusionstensor, nˆ ∈ Rd den Normaleneinheitsvektor auf dem Rand ∂ΩB, xS ∈ Rd den Saatpunkt und u˜S ∈ R 0 die initiale Zelldichte mit Standardabweichung σ ∈ R 0. Fu¨r eine detailliertere Beschreibung verweisen wir auf [ 6 ].

Die Modellierung der tumorinduzierten Gewebedeformation erfolgt basierend auf der zu einem diskreten Simulationszeitpunkt tj = qj, j = 0, . . . , m, q = τ /m, vorliegenden Zelldichte u˜. Fu¨r eine anschauliche Beschreibung fu¨hren wir die zeitdiskrete Repra¨sentation u˜j (x) := u˜(x, tj ) unter der Annahme u˜j ∈ U := L2(Ω, R 0) ein. Entsprechend markieren wir im Folgenden alle zu einem diskreten Zeitpunkt tj ausgewerteten Funktionen mit einem oberen Index j. Die Deformation wird u¨ber einen parametrischen Ansatz beschrieben. Die abbildende Funktion yj (x; µ) := y(x; µ, tj ), yj ∈ Y := C2(Ω, Rd) ist gegeben durch

υ yj (x) = x + ∑ · · · υ ∑ d ∏ bki (µi)ψlj1+k1,...,ld+kd kd=0 i=1 k1=0 Hierbei ist bki eine B-Spline-Basisfunktion des Grades υ ∈ N. Weiter ist µi = xi/hiG − ⌊xi/hiG⌋, µi ∈ R 0 die relative Position bezu¨glich des dem Parameter ψl ∈ Rd zugeo¨rigen Gitterpunktes mit Index l = (l1, . . . , ld) ∈ Zd, li = ⌊xi/hiG⌋− 1 innerhalb des Gitters G = h1GGm1G × · · · × hdGGmdG ∈ Rd(m1G+3) (mdG+3), GmiG = {li ∈ Z : −1 ≤ li ≤ miG + 1} mit Zellgro¨ße hG = (h1G, . . . , hdG) ∈ Rd>0.

Die Modellierung der Deformation des zerebralen Gewebes erfolgt u¨ber die Formulierung eines Optimierungsproblems. Gesucht wird diejenige Abbildung yj , welche die Bedingung J : Y × U × U → R, Ju~,v~[yj ] := J [yj , u˜j , v˜j ] yj Ju~,v~[yj ] = Pu~,v~[yj ] + α Cu~[yj ] −→ min erfu¨llt. Hierbei ist P : Y × U × U → R, Pu~,v~[yj ] := P[yj , u˜j , v˜j ], ein Potential, welches die Deformation antreibt, C : Y × U → R, Cu~[yj ] := C[yj , u˜j ] eine weiche Nebenbedigung und α ∈ R 0 ein Regularisierungsparameter. Das Potential in (3) ist gegeben durch (2) (3) (4) (5) Ω wobei das nicht-lineare Funktional w(u˜j (x)) = p1 exp(−p2(u˜j (x)/κ) 2 − p2(2 − (u˜j (x)/κ)2) 1) eine Wichtung (gesteuert durch p1, p2 ∈ R 0, κ ∈ R>0) in Abha¨ngigkeit der vorliegenden Zellkonzentration u˜j erlaubt (Abb. 1 (links)); v˜j repra¨sentiert eine ‘Referenzzelldichte’. Basierend auf der Annahme, dass Gewebe nahezu inkompressibel ist [ 5 ], wa¨hlen wir als weiche Nebenbedingung [ 7 ] ∫ Cu~[yj ] = β(x, u˜(x))(log(det((∂xi2 yij1 (x))i1,i2 )))2 dx.

1 ∫ 2 Ω Pu~,v~[yj ] =

(v˜j (x) − w(u˜j (x))u˜j (yj (x)))2dx Abb. 1. Wichtung des Potentials P (links; unterschiedliche Werte fu¨r p1; p2; und der Nebenbedingung C (rechts; Variation von ). = 1:0) Fu¨r eine bio-physikalisch plausible Anna¨herung der Gewebeeigenschaften fu¨hren wir das nicht-lineare Wichtungsfunktional β(x, u˜(x)) = 10ζ(x)ξ(u˜j (x)) = 10ζ(x) (2 − 2(1 + exp(−γu˜j (x))) 1) ein, wobei ζ : Ω → {−4, 0}, ζ(x) = −4 | 0 falls x ∈ ΩL | sonst mit ΩL ⊂ Ω (Gehirn-Ru¨ckenmarks-Flu¨ssigkeit). Der Trend der Funktion ξ ist fu¨r unterschiedliche γ ∈ R 0 in Abb. 1 (rechts) dargestellt.

Abschließend beleuchten wir die theoretischen Hintergru¨nde des vorgeschlagenen Deformationsmodells (Einer klaren Darstellung wegen verwerfen wir sowohl den oberen Index j als auch die Abha¨ngigkeit der Abbildung y von den Parametern ψ). Bestehende Ansa¨tze [ 4, 5 ] modellieren die Gewebedeformation basierend auf einem Kraftfeld b ∝ ∇u˜. Um einen Zusammenhang zu diesen herzustellen, bedienen wir uns der Variationsrechnung: Gegeben sei das Funktional P, Funktionen u˜, v˜ ∈ U , die Abbildung y ∈ Y sowie eine Testfunktion v ∈ Y. Gesucht ist das Gˆateaux-Differential δPu~,v~[y; v] = limε!0 1ε (Pu~,v~[y + εv] − Pu~,v~[y]) der Funktion Pu~,v~ an der Stelle y ∈ Y in Richtung der Testfunktion v ∈ V. Nutzen wir die Taylor-Entwicklung u˜(y(x) + εv) = u˜(y(x)) − ε⟨∇u˜(y(x)), v(x)⟩Rd + O(ϵ2) mit ε ∈ R 0 am Entwicklungspunkt y erhalten wir nach einigen elementaren Umformungen die Gˆateaux-Ableitung δPu~,v~[y; v] = ∫

Rd ⟨b(x, y(x)), v(x)⟩Rd dx (6) Hierbei bezeichnet b(x, y(x)) = (v˜(x) − w (u˜(x)) u˜ (y(x))) w(u˜(x))∇u˜ (y(x)) in Analogie zu variationellen Ansa¨tzen aus der Bildregistrierung [ 1 ] eine Kraft, womit die Verbindung zu bestehenden Ansa¨tzen [ 4 ] hergestellt ist. 2.2

Numerische Implementierung und Parameter

Fu¨r die Diskretisierung des Anfangsrandwertproblems (1) verwenden wir ein Euler-Cauchy-Verfahren. Die Zielfunktion (3) wird u¨ber ein Gradientenabstiegsverfahren (Schrittweite: (min., max.) = (0.1, 0.25 ∥hG∥2), max. 20 Iterationen) optimiert. Weiter gilt ∀x ∈ Rd : v˜j (x) = 0. Die experimentell bestimmten (fixen) Modellparameter sind κ = 1.0, γ = 20, α = 0.1, σ = 2.0, u˜S = 0.2, ρ = 0.012 d 1. Im Gegensatz zu [ 6 ] erfolgt die Berechnung der in der Nebenbedingung und der Optimierung auftretenden Ableitungen nicht u¨ber eine Approximation durch finite Differenzen sondern mittels analytischer Ausdru¨cke. 3

Ergebnisse

Ergebnisse fu¨r die Modellierung sind in Abb. 2 gezeigt (d = 2). Die berechneten Deformationsmuster yj sind diffeomorph. Die Parameter fu¨r die Wichtung des Potentials P wurden wie in Abb. 2 angegeben varriert (tj = 550 d). Neben diesen numerischen Experimenten stellen wir ein weiteres Simulationsergebnis einem Patientendatensatz gegenu¨ber (tj = 650 d). Die Abbildungen zeigen jeweils die berechnete Zelldichte u˜j und die berechneten Deformationsmuster yj als Gitter in u¨berlagerung. Die Farbskala fu¨r u˜j ist rechts angegeben. Die Pfeile markieren, ab welcher relativen Dichte der Tumor in T1- (orange) und T2-gewichteten (blau) Magnet-Resonanz-Aufnahmen als sichtbar angenommen werden kann.

Mang et al.

Abb. 2. Ergebnisse der Modellierung. Links: Wichtung p1 = 0:01, p2 = 0:1; Mitte: Wichtung p1 = 0:01, p2 = 0:5; Rechts: Vergleich zu realen Daten. 4

Diskussion

Die vorliegende Arbeit stellt Erweiterungen und eine theoretische Einordnung des in [ 6 ] beschriebenen Ansatzes vor. Es wurde experimentel gezeigt, dass die Einfu¨hrung nicht-linearer Wichtungsfunktionale (4), (5) eine sehr viel adaptivere Steuerung der modellierten Deformationsmuster erlaubt. Generell ist festzuhalten, dass das vorgestellte Modell nicht dafu¨r entwickelt wurde, in-vivo beobachtete Deformationsmuster exakt zu reproduzieren, sondern als bio-physikalischer Prior fu¨r eine modellbasierte Registrierung von Tumorbilddaten dient. Ein qualitativer Vergleich zu patientenindividuellen Bilddaten zeigt die Plausibilita¨t; eine quantitative Validierung des Modells ist nach dem derzeiteigen Stand der Forschung als große Herausforderung zu bewerten. Die Integration des vorgestellten Modells in einen Bildregistrieransatz ist aktueller Gegenstand unserer Forschung und bietet – durch den Vergleich zu existierenden Bildregistrieralgorithmen [ 1, 2, 3 ] – die Mo¨glichkeit einer ersten quantitativen Bewertung des vorgestellten Modells.

Literaturverzeichnis

References

1. Modersitzki J. Numerical Methods for Image Registration . Oxford University Press, New York; 2004 .

2. Fischer B , Modersitzki J . Ill-posed medicine: an introduction to image registration . Inverse Probl . 2008 ; 24 ( 034008 ): 1 - 19 .

3. Henn S , et al. Lesion preserving image registration with application to human brains . Pattern Recognit . 2004 ; 3175 : 496 - 503 .

4. Hogea C , et al. Brain-tumor interaction biophysical models for medical image registration . SIAM J Sci Comput . 2008 ; 30 ( 6 ): 3050 - 72 .

5. Hogea C , et al. A robust framework for soft tissue simulations with application to modeling brain tumor mass effect in 3D MR images . Phys Med Biol . 2007 ; 52 ( 23 ): 6893 - 908 .

6. Mang A , et al. Coupling tumor growth with brain deformation: a constrained parametric non-rigid registration problem . proc SPIE . 2010 ; 7623 : C1 - 12 .

7. Rohlfing T , et al. Volume-preserving nonrigid registration of MR breast images using free-form deformation with an incompressiblity constraint . IEEE Trans Med Imaging . 2003 ; 22 ( 6 ): 730 - 41 .