Optimalizace léčebného plánu protonové terapie

Optimalizace l´eˇcebn´eho pl´anu protonov´e terapie Soutˇeˇz o cenu akademika Baˇzanta Katedra mechaniky ˇ Stavebn´ı fakulta CVUT

Marek Tyburec [email protected] Vedouc´ı projektu: Doc. Ing. Matˇej Lepˇs Ph.D. duben 2014

Obsah ´ 1 Uvod

2

2 Popis Braggovy kˇ rivky v 1D

3

3 Vliv nehomogenit na Braggovu kˇ rivku

6

4 Popis Braggovy kˇ rivky v prostoru

8

5 Optimalizace 5.1 Popis a ˇreˇsen´ı u ´loh line´ arn´ıho programov´an´ı simplexovou metodou . . . . . . . . . . . 5.2 Pˇr´ıklady optimalizace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Sloˇzitˇejˇs´ı pˇr´ıklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11 11 11 14

6 Z´ avˇ er

16

7 Pˇ r´ılohy 7.1 Numerick´e ˇreˇsen´ı urˇcit´ ych integr´ al˚ u. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17 17

1

1

´ Uvod

ˇ Tato pr´ace vznikla v r´ amci studentsk´e vˇedeck´e ˇcinnosti na katedˇre mechaniky Stavebn´ı fakulty CVUT. Jej´ım c´ılem je postupnˇe vytvoˇrit zjednoduˇsen´ y model protonov´eho paprsku, tzv. beamletu, a pot´e sestavit optimalizovan´ y l´eˇcebn´ y pl´ an. Rakovina je nejˇcastˇejˇs´ı pˇr´ıˇcinou smrti nejen na u ´zem´ı Evropy. K jej´ı l´eˇcbˇe lze pˇristupovat r˚ uzn´ ymi zp˚ usoby, mezi kter´e se ˇrad´ı tak´e radioterapie. Radioterapie je zaloˇzena na oz´aˇren´ı zhoubn´eho n´adoru pomoc´ı ionizuj´ıc´ıho z´ aˇren´ı, pˇriˇcemˇz se minimalizuj´ı n´asledky pro okoln´ı zdravou tk´an ˇ. Nejˇcastˇeji se pouˇz´ıvaj´ı svazky elektron˚ u a foton˚ u. V posledn´ı dobˇe se ale zaˇc´ınaj´ı pouˇz´ıvat tak´e svazky hadron˚ utedy proton˚ u a lehk´ ych iont˚ u. Ionty jsou vyzaˇrov´ any z urychlovaˇce ˇc´ astic smˇerem po ose z, hloubce. D´avka z´aˇren´ı, kterou ionty odevzd´avaj´ı okol´ı, z´ avis´ı i na pouˇzit´ ych ˇc´astic´ıch. Z nich nejl´epe vych´azej´ı protony, kter´e pˇred´ avaj´ı vˇetˇsinu d´avky z´ aˇren´ı v kr´ atk´em intervalu, tzv. Braggovˇe peaku.

Obr´azek 1: Porovn´ an´ı odevzdan´e d´ avky protonov´eho a fotonov´eho paprsku v z´avislosti na hloubce. Protonov´ y paprskem m´ a v´ yrazn´ y extr´em - Bragg˚ uv peak. [15] L´ekaˇri pracuj´ıc´ı v protonov´ ych centrech maj´ı spolu se spolupracuj´ıc´ımi fyziky za u ´kol sestavit takov´ y l´eˇcebn´ y pl´ an, aby doˇslo k destrukci n´adoru pacienta a z´aroveˇ n co nejm´enˇe poˇskodili okoln´ı tk´an ˇ. Vych´az´ı pˇritom z prostorov´eho modelu sestaven´eho pomoc´ı CT (Computed Tomography), kam umist’uj´ı jednotliv´e Braggovy kˇrivky. Pot´e se optimalizuje intenzita jednotliv´ ych kˇrivek (resp. mnoˇzstv´ı proton˚ u vyj´adˇren´ ych fluenc´ı energie). Jelikoˇz se tedy nez´avisle na sobˇe ovlivˇ nuje fluence jednotliv´ ych kˇrivek, naz´ yv´a se tato metoda pencil beam, metoda tuˇzkov´ ych svazk˚ u. Tato pr´ace se postupnˇe zab´ yv´ a nejprve analytick´ ym popisem Braggovy kˇrivky v 1D, ve smˇeru ˇs´ıˇren´ı po ose z, pot´e zav´ ad´ı vliv nehomogenit prostˇred´ı, kter´ ymi se paprsek ˇs´ıˇr´ı. D´ale rozˇsiˇruje popsan´ y model do 2D, tedy uvaˇzuje prostorov´e ˇs´ıˇren´ı protonov´eho paprsku do okol´ı. Posledn´ı kapitola ukazuje jednoduch´ y zp˚ usob optimalizace 1D i 2D l´eˇcebn´eho pl´anu pˇri vyuˇzit´ı line´arn´ıho programov´an´ı.

2

2

Popis Braggovy kˇ rivky v 1D

V roce 1997 aproximoval Thomas Bortfeld tvar Braggovy kˇrivky pro protonov´e paprsky [2]. Vyˇsel pˇritom z pˇribliˇzn´eho vztahu, tzv. Bragg-Kleemanova pravidla aproximujic´ı z´avislost dosahu protonov´eho paprsku R0 na poˇc´ ateˇcn´ı energii E0 . R0 = αE0p .

(1)

Dosah R0 je definov´ an jako vzd´ alenost, kde pˇresnˇe polovina proton˚ u ztrat´ı svoji energii. Koeficienty α a p jsou materi´ alovˇe z´ avisl´e konstanty vych´azej´ıc´ı z aproximace namˇeˇren´ ych hodnot. Pro pˇresnˇejˇs´ı v´ ypoˇcet je moˇzn´e interpolovat tabulkov´e hodnoty. Vlivem ˇcasov´e i prostorov´e z´ avislosti energie ˇc´astic nemus´ı m´ıt ani protony se stejnou poˇc´ateˇcn´ı energi´ı E0 ve stejn´e hloubce z stejnou d´ avku z´aˇren´ı. Proto zav´ad´ı Bortfeld smˇerodatnou odchylku σ, kter´a ud´av´a hloubkovou nejistotu ztr´ aty energie a z´aroveˇ n zahrnuje i nejistotu dosaˇzen´ı pˇresn´e hodnoty E0 pˇri v´ ystupu ˇc´ astic z urychlovaˇce. s  2 p3 α2/p 3−2p σ= α0 R0 + (0.01E0 )2 α2 p2 E02p−2 , (2) 3p − 2 kde α0 je materi´ alov´ a konstanta z´ avisej´ıc´ı na elektronov´e hustotˇe [2]. V´ ysledn´a kˇrivka d´ avky D je sumou vˇsech Gaussov´ ych rozdˇelen´ıch d´avek z´aˇren´ı po ose z. Po u ´pravˇe integr´alu dostal Bortfeld tvar: 2

e−ζ(z,E0 ) /4 σ(E0 )1/p Γ(1/p) × D(z, E0 ) = Φ0 √ 2πρα1/p (1 + βR0 (E0 ))     1 β  P−1/p−1 (−ζ(z, E0 )) . × P (−ζ(z, E0 )) + + γβ + σ(E0 ) −1/p p R0 (E0 )

(3)

V t´eto rovnici se vyskytuje hned nˇekolik druh˚ u veliˇcin. Prvn´ı jsou konstanty z´avisl´e na materi´ alu, kam kromˇe jiˇz zm´ınˇen´ ych veliˇcin α a p patˇr´ı tak´e ρ, hustota materi´alu, , β a γ. D´ale se zde vyskytuj´ı veliˇciny z´avisl´e na hodnotˇe E0 , konkr´etnˇe zm´ınˇen´a smˇerodatn´a odchylka σ, ale tak´e ζ, kter´ a je definov´ana dle vztahu R0 (E0 ) − z . (4) σ Veliˇcina ζ je pouˇzita tak´e ve funkci parabolick´eho v´alce P [20]. Pro samotnou optimalizaci je ovˇsem nejd˚ uleˇzitˇejˇs´ı veliˇcinou Φ0 , prim´arn´ı fluence energie, popisuj´ıc´ı energii ˇc´astic, kter´e proˇsly dan´ ym m´ıstem prostoru, v tomto pˇr´ıpadˇe poˇc´atkem (z = 0). Zde je d˚ uleˇzit´e zm´ınit, ˇze d´avka D(z) je na fluenci energie line´arnˇe z´avisl´a. Vzhledem k tomu, ˇze pˇri sestavov´ an´ı l´eˇcebn´eho pl´anu je vhodn´e Braggovu kˇrivku uchopit”za jej´ı ” extr´em, tzv. peak, byla sestavena z´ avislost poˇc´ateˇcn´ı energie E0 na poloze vrcholu zmax . Nejprve byla aproximov´ana z´ avislost R0 − z na E0 , jak je vidˇet na Obr´azku 3. Do aproximace byla zahrnuta oblast poˇc´ateˇcn´ı energie E0 v rozmez´ı 1MeV aˇz 300MeV v intervalu 1MeV, celkem 300 hodnot zmax . Poloha zmax byla stanovena s pˇresnost´ı 0.1mm. Vyjdeme-li z rovnice (1) a zjiˇstˇen´e analytick´e aproximace ve tvaru ζ(z, E0 ) =

R0 − z = a1 E04 + a2 E03 + a3 E02 + a4 E0 + a5 ,

(5)

kde a1 aˇz a5 jsou ˇc´ıseln´e konstanty, lze v´ yslednˇe rovnice upravit na tvar zmax = αE0p − a1 E04 − a2 E03 − a3 E02 − a4 E0 − a5 .

(6)

Porovnaj´ı-li se takto spoˇcten´e hodnoty zmax s hodnotami v´ yˇse zjiˇstˇen´ ymi s pˇresnost´ı 0.1mm, je zde maxim´aln´ı hodnota chyby 0.0646mm a jej´ı smˇerodatn´a odchylka 0.0288mm. V´ ysledky tedy odpov´ıdaj´ı zvolen´e pˇresnosti.

3

Takto je jiˇz moˇzn´e rozm´ıstit v 1D prostoru Braggovy kˇrivky tak, aby mˇely v urˇcit´ ych bodech extr´em. Konkr´etn´ı hodnoty konstantn´ıch veliˇcin platn´ ych pro vodu jsou uvedeny opˇet ve zm´ınˇen´e pr´aci T. Bortfelda [2]. 100

80

60

40

20

0

0

2

4

6

8

10

12

14

16

18

20

Obr´azek 2: Zn´azornˇen´ı z´ akladn´ıch parametr˚ u Braggovy kˇrivky s E0 = 158.6MeV - dosahu R0 a polohy extr´emu zmax .

1 10−6 (0.000027E04 − 0.020207E03 + 14.126506E02 + 170.637718E0 + 1725.877297) Fitované hodnoty 0.8

R0 − zmax [cm]

0.6

0.4

0.2

0

−0.2

0

50

100

150 E0 [MeV]

200

250

300

Obr´azek 3: Aproximovan´ a z´ avislost R0 − zmax na E0 se zobrazen´ ymi fitovan´ ymi hodnotami R0 − zmax . Z´aˇren´ı ve vodˇe.

4

−9

6

x 10

100 140 180 200

5

MeV MeV MeV MeV

D´ avka D [Gy]

4

3

2

1

0

0

5

10

15 z [cm]

20

25

30

Obr´azek 4: Vliv rozd´ıln´e poˇc´ ateˇcn´ı energie E0 na v´ yslednou d´avku1 . Z´aˇren´ı ve vodˇe.

100 100 140 180 200

90

Pomˇern´ a d´ avka D/Dmax [%]

80

MeV MeV MeV MeV

70 60 50 40 30 20 10 0

0

5

10

15 z [cm]

20

25

30

Obr´azek 5: Vliv rozd´ıln´e poˇc´ ateˇcn´ı energie E0 na pomˇernou d´avku. Dle obr´azku je zˇrejm´e, ˇze pro co nejmenˇs´ı zasaˇzen´ı okoln´ı tk´ anˇe je nejv´ yhodnˇejˇs´ı ozaˇrovat co nejkratˇs´ı cestou.

1

Jeden Gray (Gy) odpov´ıd´ a absorbci energie 1J v 1kg l´ atky.

5

3

Vliv nehomogenit na Braggovu kˇ rivku

Prostˇred´ı n´adoru v lidsk´em organismu nen´ı ovˇsem obklopeno vodou. Tak´e nen´ı homogenn´ı. Z toho vypl´ yv´a d˚ uleˇzitost sestaven´ı v´ ypoˇctov´eho modelu tak´e pro heterogenn´ı materi´aly a smˇesi. Fyzici zde vyuˇzili jedn´e vlastnosti - podobnosti ˇs´ıˇren´ı protonov´eho z´aˇren´ı ve vodˇe a v tk´ani. D´ıky tomu je moˇzn´e zav´est veliˇcinu W ET , tzv. vodn´ı ekvivalent, ud´avaj´ıc´ı takovou tlouˇst’ku vodn´ı vrstvy, kter´a by mnoˇzstv´ım ztr´ aty energie nahradila tlouˇst’ku dan´e nehomogenity. Jak uˇz bylo zm´ınˇeno dˇr´ıve, ztr´ata energie nen´ı konstantn´ı. Proto i hodnota vodn´ıho ekvivalentu z´avis´ı nejen na materi´alu, ale i na energii v m´ıstˇe poˇc´ atku nehomogenity a na jej´ı tlouˇst’ce.

D [Gy]

3 2 1 0

0

5

10

15

20

25

30

35

20

25

30

35

20

25

30

35

z [cm]

D [Gy]

3 2 1 0

0

5

10

15 z [cm]

D [Gy]

3 2 1 0

0

5

10

15 z [cm]

Obr´azek 6: Vliv nehomogenity na Braggovu kˇrivku. Prvn´ı kˇrivka ukazuje Braggovu kˇrivku s vrcholem ve vzd´alenosti 25cm ve vodˇe, bez vlivu nehomogenit. Druh´a kˇrivka zn´azorˇ nuje opˇet vrchol ve vzd´alenosti 25cm, ale mezi 5 a 15cm je um´ıstˇena nehomogenita, kost. Tˇret´ı kˇrivka m´a stejnou poˇc´ateˇcn´ı energii E0 jako kˇrivka druh´ a, ale nen´ı zde vliv nehomogenity. Na prav´e i lev´e stranˇe od ˇsed´ ych oblast´ı jsou druh´a a tˇret´ı kˇrivka stejn´e. Vyjde-li se opˇet z Bragg-Kleemanova pravidla (1), lze z nˇej vyj´adˇrit pˇribliˇzn´ y vztah pro funkci energie za pˇredpokladu, ˇze hodnota dosahu R0 se po z mˇen´ı line´arnˇe. Potom:  R0 − z 1/p . (7) R0 − z = αE(z) ≡ E(z) = α Tento vztah je znaˇcnˇe zjednoduˇsen´ y. Zanedb´av´a zejm´ena hloubkovou nejistou danou dle rovnice (2). Nen´ı tedy platn´ y v oblasti vrcholu. Pokud je jiˇz zn´ am´ a energie v dan´e pozici, pomoc´ı speci´aln´ı teorie relativity Alberta Einsteina lze z´ıskat hodnotu rychlosti protonu [6]: v v u 1 β= =u (8) 2 , t1 −  c 1 + mEp c2 

p

6

kde mp je hmotnost protonu, E je kinetick´a energie a c rychlost svˇetla. Potom Lorentz˚ uv faktor γ je: s γ=

r

1 1−

v2 c2

=

1 . 1 − β2

(9)

Pro vypoˇcten´ı vodn´ıho ekvivalentu byla pouˇzita Bethe-Blochova rovnice [21] a upravena na tvar:    m Zef f 2me c2 γ(E)2 β(E)2 2 , W ET (E) ≈ tw = tm ρ ln − β(E) (10) Aef f Ief f w kde tw oznaˇcuje ekvivalentn´ı tlouˇst’ku vody, tm re´alnou tlouˇst’ku nehomogenity, ρ hustotu, Aef f efektivn´ı nukleonov´e ˇc´ıslo, Zef f efektivn´ı atomov´e (protonov´e) ˇc´ıslo, me hmotnost elektronu a Ief f efektivn´ı excitaˇcn´ı energii. Pokud se nehomogenita skl´ ad´ a z materi´alu sloˇzen´eho pouze z jednoho chemick´eho prvku, potom se m´ısto efektivn´ıch hodnot dosazuj´ı konkr´etn´ı hodnoty dotyˇcn´eho prvku. V pˇr´ıpadˇe smˇesi se pro zjednoduˇsen´ı v´ ypoˇct˚ u zav´ ad´ı efektivn´ı atom”, tj. atom, kter´ y m´a stejnou elektronovou hustotu. Potom ” se efektivn´ı nukleonov´e ˇc´ıslo spoˇc´ıt´ a jako [12]: Pn Ni A2i , (11) Aef f = Pi=1 n i=1 Ni Ai kde n je poˇcet prvk˚ u, ze kter´ ych se smˇes skl´ad´a a N je poˇcet atom˚ u dan´eho prvku. Efektivn´ı protonov´e ˇc´ıslo se spoˇc´ıt´ a [12]: Zef f = Aef f

n X Ni Zi . Ni Ai

(12)

i=1

Hodnota excitaˇcn´ı energie je pomˇernˇe nejist´a, jelikoˇz ji nelze pˇresnˇe stanovit. V ˇcl´anku Zhanga a Newhausera [21] je uveden vztah: I = kZ,

(13)

  14.5, pro Z ≤ 8 k = 13, pro 8 < Z ≤ 13 .   11, pro Z > 13

(14)

Pro smˇesi plat´ı vztah [3]: ln Ief f

Pn i=1 Ni Zi ln Ii = P . n i=1 Ni Zi

(15)

Vzhledem k line´ arn´ı z´ avislosti W ET na tlouˇst’ce materi´alu tm (10) a pˇredpokladu, ˇze se po ose z mˇen´ı energie, je zˇrejm´e, ˇze Bethe-Blochova rovnice plat´ı pouze pro infinitezim´aln´ı tlouˇst’ky. V mnou proveden´e implementaci bylo vˇetˇs´ı pˇresnosti dosaˇzeno rozdˇelen´ım intervalu dle zadan´e pˇresnosti. Materi´al Al

E [MeV] 200 200 200 100 100

tm [mm] 19.73 14.90 4.83 14.90 4.83

W ETexp [mm] 42.30 32.20 10.04 31.50 10.30

W ETBK [mm] 41.7 31.5 10.2 31.1 10.1

Chybaexp−BK [%] -1.42 -2.17 +1.59 -1.27 -1.94

W ET [mm] 41.861 31.614 10.249 31.357 10.166

Chybaexp−W ET [%] -1.038 -1.820 +2.082 -0.454 -1.301

Tabulka 1: Porovn´ an´ı pˇresnosti v´ ypoˇctu W ET . W ETexp ud´av´a experiment´alnˇe zmˇeˇrenou hodnotu vodn´ıho ekvivalentu pro zadan´e charakteristiky, W ETBK je hodnota vypoˇcten´a podle [22], hodnota W ET je vypoˇctena podle implementace popsan´e v t´eto pr´aci. Chybaexp−BK ud´av´a rozd´ıl mezi W ETexp a W ETBK , Chybaexp−W ET ud´ av´ a rozd´ıl mezi W ETexp a W ET . 7

4

Popis Braggovy kˇ rivky v prostoru

V letech 1947 a 1948 napsal Moli´ere z´ asadn´ı teorii, kter´a popisuje ˇs´ıˇren´ı ˇc´astic prostorem. Jednalo se nejprve o teorii jednoduch´eho [16] a pozdˇeji mnohon´asobn´eho rozptylu [17]. Tato teorie ovˇsem vyˇsla pouze nˇemecky. V angliˇctinˇe se probl´emem zab´ yval H.A. Bethe [1], kter´ y v roce 1952 vydal ˇcl´ anek Molier’s Theory of Multiple Scattering, kde ovˇsem zanedbal nˇekter´a zobecnˇen´ı platn´a pro smˇesi [18]. Moli´erova teorie je obecnˇe br´ ana jako nejpˇresnˇejˇs´ı, ovˇsem tak´e nejsloˇzitˇejˇs´ı. Nastala proto potˇreba teorii zjednoduˇsit. Tohoto kroku se ujal v roce 1975 Virgil L. Highland ve sv´e pr´aci Some practical Remarks on Multiple Scattering [10], kde zavedl jednoduchou Gaussovskou aproximaci, ˇc´ımˇz v´ ypoˇcet znaˇcnˇe zjednoduˇsil a zachoval dostateˇcnou pˇresnost. Nicm´enˇe tato aproximace platila pouze pro mal´e tlouˇst’ky, a tak byla v roce 1993 zobecnˇena B. Gottschalkem [8]. V´ ysledn´ y tvar aproximace je tento:   "Z t  2 0 #1/2 1 1 t dt θ0 = 14.1z 1 + log10 × , 9 LR pv LR 0

(16)

kde θ0 je charakteristick´ y u ´hel mnohon´asobn´eho rozptylu. Veliˇcina z v tomto pˇr´ıpadˇe oznaˇcuje n´aboj ˇc´astice, v pˇr´ıpadˇe proton˚ u je to 1eV, t znaˇc´ı tlouˇst’ku materi´alu a LR radiaˇcn´ı d´elku. p je hybnost ˇc´astice a v jej´ı rychlost. Rychlost se spoˇc´ıt´a z kinetick´e energie podle rovnice (8). Hybnost ˇc´astice se podle speci´ aln´ı teorie relativity spoˇc´ıt´a jako [6]: p = γmp v,

(17)

kde γ je Lorentz˚ uv faktor spoˇcten´ y dle rovnice (9) a mp je hmotnost protonu. Posledn´ı nezn´ amou veliˇcinou je radiaˇcn´ı d´elka LR . Tu lze podle [9] vyj´adˇrit takto: LR =

716.4 287 . Z(Z + 1) ln √ Z

(18)

Pro v´ ypoˇcet radiaˇcn´ı d´elky pro deskovou skladbu se pouˇz´ıv´a vztah: t0 ρ0 t1 ρ1 t2 ρ2 = + , LR LR,1 LR,2

(19)

kde ρ je hustota materi´ alu. Z´ aroveˇ n plat´ı vztah: t0 ρ0 = t1 ρ1 + t2 ρ2 .

(20)

Pro smˇesi se ud´ av´ a vztah: n

A0 N0 X Ai Ni = , LR LR,i

(21)

i=1

kde n je poˇcet prvk˚ u, ze kter´ ych se smˇes skl´ad´a, N je mol´arn´ı mnoˇzstv´ı a A je nukleonov´e ˇc´ıslo. Pot´e je jiˇz moˇzn´e pˇristoupit k v´ ypoˇctu integr´alu v rovnici (16), k jehoˇz v´ ypoˇctu byla pouˇzita Newtonova metoda 3/8 (viz Pˇr´ıloha). V literatuˇre [8] se uv´ad´ı pouˇzit´ı Simpsonova pravidla a rozdˇelen´ı tlouˇst’ky t pravidelnˇe podle dosahu R0 . Tvar Gaussovsk´eho rozdˇelen´ı je d´ an vztahem: h

1 f (θ) = 2πθ02

−1 θ 2 2 θ0

i

,

(22)

kde θ je u ´hel mnohon´ asobn´eho rozptylu dle Obr´azku 7. Pˇredpokl´adejme, ˇze v prostoru se protonov´ y paprsek ˇs´ıˇr´ı tedy kromˇe smˇeru z tak´e do smˇer˚ u kolm´ ych, konkr´etnˇe osy y a x. V pˇr´ıpadˇe 2D modelu budeme uvaˇzovat jen smˇery z a y, jak je patrno na Obr´azku 7. Z nˇej rovnˇeˇz vypl´ yv´ a geometrick´a z´avislost veliˇcin: tan θ = 8

y . z

(23)

y θ

z

Obr´ azek 7: Z´avislost u ´hlu θ, hloubky z a y. Takto lze potom jednoduˇse pro kaˇzd´e m´ısto v prostoru dopoˇc´ıtat d´avku z´aˇren´ı: "

D(y, z) = D(z) ×

1 2πθ02

−1 2



y arctan z θ0

2 #

.

(24)

−4

y [cm]

−2

0

2

4

0

5

10

15

20

25

z [cm]

Obr´ azek 8: Grafick´e zn´ azornˇen´ı vlivu mnohon´asobn´eho rozptylu. V Tabulce 2 je uvedeno porovn´ an´ı v´ ypoˇctu hodnot podle Moli´era, Highlandovy aproximace podle [8] a hodnot vypoˇcten´ ych dle zde popsan´eho postupu. Lze si povˇsimnout, ˇze hodnoty θ0 a θH se po vˇetˇsinˇe intervalu v´ yraznˇe neodliˇsuj´ı, probl´em nast´av´a v pˇr´ıpadˇe, kdy se tlouˇst’ka pˇribliˇzuje dosahu R0 , jelikoˇz Bragg-Kleemanova rovnice (1) neuvaˇzuje hloubkovou nejistotu danou smˇerodatnou odchylkou (2), tud´ıˇz vˇsechny protony ztr´ ac´ı energii v m´ıstˇe dosahu R0 . Zlepˇsen´ı by mohlo pˇrin´est pouˇzit´ı interpolace z tabulkov´ ych hodnot.

9

D [Gy]

4 2 0 0

2

4

6

8

10

12

14

−2 16

18

0 20

z [cm]

2 22

y [cm]

Obr´ azek 9: Vliv mnohon´ asobn´eho rozptylu na intenzitu d´avky z´aˇren´ı. Materi´ al Al

Tlouˇst’ka [g/cm2 ] 0.2160 0.8170 2.1729 3.3500 7.0960 11.9570 13.5690 17.7230 21.2450 21.9150 22.1100 22.3300

θM [mrad] 3.701 8.051 13.880 16.920 28.357 42.065 42.422 61.129 91.129 92.504 98.021 98.256

θH [mrad] 3.534 7.670 13.104 16.258 26.931 39.986 40.534 58.230 87.103 88.657 93.645 94.390

ChybaM −H [%] -4.512 -4.732 -5.591 -3.913 -5.029 -4.942 -4.451 -4.742 -4.418 -4.159 -4.464 -3.935

θ0 [mrad] 3.500 7.428 13.024 16.823 26.976 39.607 44.231 59.329 87.279 104.860 117.185 chyba

ChybaM −0 [%] -5.431 -7.738 -6.167 -0.573 -4.870 -5.843 +4.271 -3.092 -4.225 +16.600 +19.551 chyba

Tabulka 2: Hodnoty mnohon´ asobn´eho rozptylu. θM je nejpˇresnˇeji stanoven´a hodnota dle Moli´erovy teorie mnohon´asobn´eho rozptylu. Hodnota u ´hlu θH je uv´adˇena v [8] jako v´ ysledek Highlandovy apro´ ximace. Uhel θ0 je poˇc´ıtan´ y dle popsan´e implementace. ChybaM −H ud´av´a procentu´aln´ı rozd´ıl hodnot θM a θH , ChybaM −0 ud´ av´ a procentu´ aln´ı rozd´ıl hodnot θM a θ0 .

10

5

Optimalizace

C´ılem optimalizace je minimalizace d´ avky z´aˇren´ı D tak, aby doch´azelo k co nejmenˇs´ımu poˇskozen´ı okoln´ı tk´anˇe a z´ aroveˇ n byl zniˇcen n´ ador. Optimalizovanou veliˇcinou je fluence energie Φ0 popisuj´ıc´ı jednotliv´e Braggovy kˇrivky, na kter´ ych je d´avka D line´arnˇe z´avisl´a. Je tedy moˇzn´e pouˇz´ıt line´ arn´ı programov´an´ı.

5.1

Popis a ˇ reˇ sen´ı u ´ loh line´ arn´ıho programov´ an´ı simplexovou metodou

´ Ulohou line´arn´ıho programov´ an´ı se naz´ yv´ au ´loha, kde je d´ana line´arn´ı n-dimenzion´aln´ı u ´ˇcelov´a funkce f (x) = c1 x1 + c2 x2 + ... + cn xn ,

(25)

´ celov´a funkce se bud’ minimalizuje nebo maximalizuje. kde c1 aˇz cn se naz´ yvaj´ı cenov´e koeficienty. Uˇ D´ale se zav´ad´ı line´ arn´ı omezuj´ıc´ı podm´ınky pro jednotliv´e promˇenn´e. Obecnˇe maj´ı tvar:

a1,1 x1 + a1,2 x2 + ... + a1,n xn ≤ b1 a2,1 x1 + a2,2 x2 + ... + a2,n xn ≤ b2 .. .

(26)

am,1 x1 + am,2 x2 + ... + am,n xn ≤ bm . ˇ sen´ı u Reˇ ´loh line´ arn´ıho programov´ an´ı m˚ uˇze prob´ıhat napˇr. pomoc´ı simplexov´e metody, metody vnitˇrn´ıch bod˚ u a elipsoidov´e metody. Zde se zab´ yv´am pouze simplexovou metodou. Simplexov´a metoda stav´ı na pˇredpokladu, ˇze m omezuj´ıc´ıch podm´ınek oˇreˇze”n-dimension´aln´ı pro” stor na simplex o m hraniˇcn´ıch bodech. Uvnitˇr a na hran´ach simplexu se nach´az´ı vˇsechna pˇr´ıpustn´ a ˇreˇsen´ı. Optimum pak leˇz´ı v nˇekter´em z hraniˇcn´ıch bod˚ u, pˇr´ıpadnˇe na cel´e hranˇe simplexu. C´ılem simplexov´e metody je postupovat od jednoho bodu pˇr´ıpustn´eho ˇreˇsen´ı tak, aby doch´azelo k minimalizaci c´ılov´e funkce (pokud hled´ ame minimum) po hranˇe k dalˇs´ımu. V pˇr´ıpadˇe, ˇze ˇz´adn´ y sousedn´ı hraniˇcn´ı bod nenab´ız´ı zlepˇsen´ı hodnot u ´ˇcelov´e funkce, nach´az´ıme se v glob´aln´ım optimu. Pro v´ ypoˇcet u ´loh line´ arn´ıho programov´an´ı je pouˇzit software Matlab, prvn´ı pˇr´ıklad je pro n´azornost ˇreˇsen graficky, ostatn´ı pomoc´ı funkce linprog.

5.2

Pˇ r´ıklady optimalizace

Aby bylo moˇzn´e uv´est jednoduch´ y pˇr´ıklad optimalizace, je nutn´e nejprve zav´est pˇredpoklady, se kter´ ymi se konkr´etnˇe do tˇechto u ´loh vstupuje. Tk´an ˇ je rozdˇelena v podstatˇe do tˇrech podmnoˇzin. Prvn´ı tvoˇr´ı T ARGET , coˇz je m´ısto n´ adoru. Zde je l´ekaˇri pˇredeps´ana nˇejak´a minim´aln´ı d´avka oz´ aˇren´ı T ARGETmin , kter´ a mus´ı b´ yt splnˇena, aby doˇslo k destrukci n´adoru. Z´aroveˇ n zde b´ yv´a ud´av´an jeˇstˇe horn´ı limit T ARGETmax , kter´ y by nemˇel b´ yt pˇrekroˇcen. Druhou podmnoˇzinou je OAR (organ at risk), tedy takov´a prostorov´a ˇc´ast tk´anˇe, kde je um´ıstˇen d˚ uleˇzit´ y org´an, kter´ y chceme zas´ ahnout co nejm´enˇe. Opˇet je zde pˇredepsan´ y urˇcit´ y limit oz´ aˇren´ı OARmax , kter´ y se nem´ a pˇrekroˇcit. Posledn´ı oblast´ı je potom zbytek tk´anˇe, kter´a by mˇela b´ yt zasaˇzen´ a samozˇrejmˇe co nejm´enˇe, ale nen´ı zde pˇredeps´an ˇz´adn´ y limit. Vzhledem k vyuˇzit´ı line´ arn´ ıho programov´ a n´ ı byla zvolena jako u ´ˇcelov´a funkce minimalizace mnoˇzP stv´ı vyz´aˇren´ ych proton˚ u min Φ0 . Z v´ yˇse zm´ınˇen´eho lze jiˇz obecnˇe formulovat tvar u ´lohy liner´arn´ıho programov´an´ı:

11

min

n X

Φ0,i

i=1

s.t.:

n X i=1 n X

Di (Φ0,i , zOAR,i ) ≤ DOAR,max Di (Φ0,i , zT ARGET,i ) ≥ DT ARGET,min

i=1 n X

Di (Φ0,i , zT ARGET,i ) i=1 ∀ni=1 Φ0,i ≥ 0,

≤ DT ARGET,max (27)

kde D• je d´ avka v m´ıstˇe • a n znaˇc´ı poˇcet Braggov´ ych kˇrivek. Jako jednoduch´ y pˇr´ıklad optimalizace poslouˇz´ı u ´loha o dvou promˇenn´ ych Φ1 a Φ2 , viz Obr´azek 11. Umist’uj´ı se pouze dvˇe Braggovy kˇrivky s vrcholy ve vzd´alenostech 24 a 28cm. Ve vzd´alenosti 10-15cm je zad´an OAR a hodnota OARmax = 5Gy. T ARGET se nach´az´ı mezi 20 a 28cm. T ARGETmax = 10Gy a T ARGETmin = 5Gy. Ovˇeˇrovat splnˇen´ı podm´ınek budeme po 1cm. Nejprve se z´ıskaj´ı funkˇcn´ı hodnoty (d´ avky) jednotliv´ ych Braggov´ ych kˇrivek v bodech, kde se ovˇeˇruj´ı podm´ınky: Kˇrivka 1 2

OAR(15) 1.089 0.978

T(20) 1.330 1.066

T(21) 1.447 1.098

T(22) 1.623 1.139

T(23) 2.073 1.195

T(24) 2.940 1.271

T(25) 0.798 1.383

T(26) 0.007 1.568

T(27) 0.000 1.978

T(28) 0.000 2.657

Tabulka 3: Vypoˇcten´e d´ avky Braggov´ ych kˇrivek v bodech, kde se ovˇeˇruj´ı podm´ınky. Pot´e se sestav´ı jednotliv´e omezuj´ıc´ı podm´ınky:

OAR: 1) 1.089Φ1 + 0.278Φ2 ≤ 5 TARGET: 2) 1.330Φ1 + 1.066Φ2 ≥ 5

11) 1.330Φ1 + 1.066Φ2 ≤ 10

3) 1.447Φ1 + 1.098Φ2 ≥ 5

12) 1.447Φ1 + 1.098Φ2 ≤ 10

4) 1.623Φ1 + 1.139Φ2 ≥ 5

13) 1.623Φ1 + 1.139Φ2 ≤ 10

5) 2.073Φ1 + 1.195Φ2 ≥ 5

14) 2.073Φ1 + 1.195Φ2 ≤ 10

6) 2.940Φ1 + 1.271Φ2 ≥ 5

15) 2.940Φ1 + 1.271Φ2 ≤ 10

7) 0.798Φ1 + 1.383Φ2 ≥ 5

16) 0.798Φ1 + 1.383Φ2 ≤ 10

8) 0.007Φ1 + 1.568Φ2 ≥ 5

17) 0.007Φ1 + 1.568Φ2 ≤ 10

9) 0.000Φ1 + 1.978Φ2 ≥ 5

18) 0.000Φ1 + 1.978Φ2 ≤ 10

10) 0.000Φ1 + 2.657Φ2 ≥ 5

19) 0.000Φ1 + 2.657Φ2 ≤ 10

Nez´apornost fluenc´ı: 1.000Φ1 + 0.000Φ2 ≥ 0 0.000Φ1 + 1.000Φ2 ≥ 0.

(28)

´ celov´a funkce je: Uˇ f = min (Φ1 + Φ2 ).

(29)

Takto jednoduch´ au ´loha by ˇsla ˇreˇsit pouze graficky. Z´akladn´ı mnoˇzinu tvoˇr´ı rovina x ∈ R a y ∈ R. Tato mnoˇzina je pot´e oˇrez´ ana plnˇen´ım omezuj´ıc´ıch podm´ınek. Jak je vidˇet z obr´ azku, v tomto pˇr´ıpadˇe je simplex ohraniˇcen ˇctyˇrmi aktivn´ımi omezuj´ıc´ımi podm´ınkami (1, 10, 11 a 17), jejichˇz pr˚ useˇc´ıky vytv´aˇr´ı ˇctyˇri hraniˇcn´ı body:

12

9 8

2

7

4

7 8

6

5 9

5

Φ2

3

6

1

4

10 11

3

17

12

18

2

16

19

13

14

1

15

0 −1 0.5

1

1.5

Φ1

2

Obr´azek 10: Zn´ azornˇen´ı omezuj´ıc´ıch podm´ınek. Znaˇcen´ı jednotliv´ ych omezuj´ıc´ıch podm´ınek na tomto ˇ obr´azku odpov´ıd´ a jejich znaˇcen´ı v rovnici (28). Sedou plnou ˇsrafou je vyznaˇcena oblast, kter´a tvoˇr´ı ˇreˇsen´ı dan´e soustavy nerovnic. Oznaˇcen´ı bodu A1 A2 A3 A4

Φ1 0.744 1.210 1.212 1.735

Φ2 3.764 3.183 3.764 3.181

Funkˇcn´ı hodnota 4.508 4.393 4.976 4.916

Tabulka 4: Krajn´ı body simplexu a jejich funkˇcn´ı hodnota u ´ˇcelov´e funkce. Dle hodnot u ´ˇcelov´e funkce ve vypoˇcten´ ych bodech je zˇrejm´e, ˇze se glob´aln´ı optimum nach´ az´ı v bodˇe A2 s hodnotou u ´ˇcelov´e funkce 4.393. Znamen´a to tedy, ˇze glob´alnˇe-optim´aln´ım ˇreˇsen´ım, kter´e splˇ nuje vˇsechny omezuj´ıc´ı podm´ınky, je Φ1 -n´asobek Braggovy kˇrivky s vrcholem ve 24cm a Φ2 -n´asobek Braggovy kˇrivky s vrcholem v 28cm. 10 OAR Target Suma d´ a vek D´ a vky

D [Gy]

8 6 4 2 0

0

5

10

15

20

25

30

35

z [cm]

Obr´azek 11: Jednoduch´ a optimalizaˇcn´ı u ´loha - zn´azornˇen´ı zad´an´ı um´ıstˇen´ı vrchol˚ u a v´ ysledn´eho optimalizovan´eho ˇreˇsen´ı. Bˇeˇznˇe se ovˇsem u ´lohy line´ arn´ıho programov´an´ı neˇreˇs´ı graficky, jelikoˇz maj´ı v´ıce neˇz 2 nebo i 3 dimenze, tud´ıˇz je nelze takto snadno zn´ azornit. 13

100

Histogram TARGET

Objem [%]

80 60 40 20 0

0

1

2

3

4

5

6

7

8

9

D [Gy]

Obr´azek 12: Jednoduch´ a optimalizaˇcn´ı u ´loha - histogram. Zobrazen´ı, v jak´em procentu´aln´ım mnoˇzstv´ı T ARGET u bylo dosaˇzeno alespoˇ n urˇcit´e d´avky.

5.3

Sloˇ zitˇ ejˇ s´ı pˇ r´ıklady

Druh´ ym pˇr´ıkladem, kter´ y zde uv´ ad´ım, je u ´loha s vlivem nehomogenit a ozaˇrov´an´ım z obou stran. Jeden zdroj z´aˇren´ı je um´ıstˇen v z = 0cm, druh´ y zdroj v z = 50cm. Ve vzd´alenosti 10-15cm je uvaˇzov´an OAR s maxim´ aln´ı d´ avkou OARmax = 1Gy. T ARGET ve vzd´alenosti 20-30cm je omezen T ARGETmax = 10Gy a T ARGETmin = 5Gy. Mezi 10-15cm a 35-40cm se nach´az´ı kost.

D [Gy]

10 OAR TARGET Suma d´ a vek D´ a vky 5

0

0

5

10

15

20

25 z [cm]

30

35

40

45

50

Objem [%]

100

50 Histogram TARGET 0

0

1

2

3 D [Gy]

4

5

6

Obr´ azek 13: Sloˇzitˇejˇs´ı u ´loha s vlivem nehomogenit a ozaˇrov´an´ım z obou stran. Posledn´ım uveden´ ym pˇr´ıkladem je 2D u ´loha s vlivem mnohon´asobn´eho rozptylu a ozaˇrov´an´ım z obou stran. Rozmˇery jsou zˇrejm´e z Obr´ azku 15 - p˚ udorysu. OARmax = 1Gy, T ARGETmin = 5Gy, T ARGETmax = 10Gy

14

D [Gy]

10

5

0 50 45 40 35 30 25 20 15 z [cm]

10

10 5

5 0

0

y [cm]

Obr´azek 14: 2D optimalizaˇcn´ı u ´loha s ozaˇrov´an´ım z obou stran, bez vlivu nehomogenit.

12

y [cm]

10 8

OAR TARGET Obrys poˇc´ıtan´e oblasti

6 4 2 0 50

45

40

35

30

25 z [cm]

20

15

10

5

Obr´azek 15: P˚ udorys 2D optimalizaˇcn´ı u ´lohy s vyznaˇcen´ım T ARGET u a OAR.

15

0

6

Z´ avˇ er

C´ılem pr´ace bylo vytvoˇrit zjednoduˇsen´ y model ˇs´ıˇren´ı protonov´eho z´aˇren´ı prostˇred´ım. Implementov´ ano bylo prozat´ım ˇs´ıˇren´ı proton˚ u ve 2D prostoru za pˇr´ıtomnosti nehomogenit a ozaˇrov´an´ı z obou stran. Zjednoduˇsen´ y model splˇ nuje moˇzn´e pˇredem stanoven´e odchylky ve vˇetˇsinˇe pˇr´ıpadech. V pˇr´ıpadˇe v´ ypoˇctu vodn´ıho ekvivalentu doˇslo i k m´ırn´emu zlepˇsen´ı pˇresnosti oproti ˇcl´anku, ze kter´eho bylo ˇcerp´ano. Dalˇs´ı zlepˇsen´ı by jistˇe pˇrineslo pouˇzit´ı interpolace tabulkov´ ych hodnot Bragg-Kleemanova pravidla. N´asleduj´ıc´ım krokem bude vytvoˇren´ı plnˇe prostorov´eho modelu, zaveden´ı r˚ uzn´ ych smˇer˚ u z´ aˇren´ı au ´pravy v optimalizaci tak, aby v´ıce odpov´ıdala pouˇz´ıvan´ ym metod´am. Konkr´etnˇe u ´prava u ´ˇcelov´e funkce a typu nezn´ am´ ych veliˇcin. V´ ysledkem line´ arn´ıho programov´ an´ı je hodnota glob´aln´ıho optima, coˇz se v pˇr´ıpadˇe protonov´e l´eˇcby velice hod´ı, jelikoˇz nelze nal´ezt lepˇs´ı pl´an. Nev´ yhodou ovˇsem je, ˇze m˚ uˇze nastat situace, kdy v z´avislosti na omezuj´ıc´ıch podm´ınk´ ach ˇreˇsen´ı v˚ ubec neexistuje. Potom line´arn´ı programov´an´ı selˇze. ˇ sen´ım je zaveden´ı nezn´ Reˇ am´ ych v podobˇe odchylkov´ ych promˇenn´ ych δ [13]. Potom m´a u ´loha jiˇz vˇzdy ˇreˇsen´ı. Aby bylo moˇzn´e ovlivnit d˚ uleˇzitost jednotliv´ ych mez´ı, sestav´ı se v´ıcekriteri´aln´ı model line´arn´ıho programov´an´ı s u ´ˇcelov´ ymi funkcemi dan´ ymi tˇemito promˇenn´ ymi. Pot´e se sestav´ı Pareto-mnoˇzina, podle n´ıˇz se vˇzdy vybere ten nejlepˇs´ı pl´ an.

16

7 7.1

Pˇ r´ılohy Numerick´ eˇ reˇ sen´ı urˇ cit´ ych integr´ al˚ u

V pˇr´ıpadech, kdy se ned´ a ˇreˇsen´ı urˇcit´eho integr´alu vyj´adˇrit symbolick´ ym z´apisem, je moˇzn´e pˇr´ısluˇsnou hodnotu spoˇc´ıtat numericky. Jedin´e, co je tˇreba zn´at, jsou funkˇcn´ı hodnoty poˇzadovan´e funkce v urˇcen´ ych bodech. Potom, dle poˇzadovan´e pˇresnosti, se tyto body proloˇz´ı polynomi´aln´ı funkc´ı, jej´ıˇz integr´al se spoˇc´ıt´ a snadno. Obecnˇe se, dle poˇzadovan´e pˇresnosti, rozdˇel´ı interval na nˇekolik ˇc´ast´ı. V kaˇzd´e ˇc´asti intervalu se spoˇc´ıt´a n funkˇcn´ıch hodnot, kde n znaˇc´ı stupeˇ n polynomu, kter´ ym se funkce aproximuje. Bˇeˇznˇe se pouˇz´ıvaj´ı polynomy do 4. stupnˇe, v´ ypoˇcet s n ≥ 8 je uˇz ale problematick´ y, jelikoˇz koeficienty c nab´ yvaj´ı uˇz i z´aporn´ ych hodnot. [14] Z

b

f (x)dx ≈ a

n 0 1 2 3 4

b−a [c1 f (x0 ) + ... + cn f (xn )]. n

N´ azev pravidla Obd´eln´ıkov´e Lichobˇeˇzn´ıkov´e Simpsonovo Newtonovo 3/8 Milneovo

Koeficienty c 1 1 2

14 45

17

1 2

1 4 1 3 3 3 3 9 9 3 8 8 8 8 64 24 64 45 45 45

14 45

(30)

Reference [1] H. Bethe, Moliere’s theory of multiple scattering, Physical Review, 89 (1953), p. 1256. [2] T. Bortfeld, An analytical approximation of the bragg curve for therapeutic proton beams, Medical physics, 24 (1997), pp. 2024–2033. [3] J. Coderre, Principles of radiation interactions. http://ocw.mit.edu/courses/ nuclear-engineering/22-55j-principles-of-radiation-interactions-fall-2004/ lecture-notes/energy_depos_hcp.pdf. [Online; pˇr´ıstup 08-04-2014]. [4] P. de Vera, I. Abril, and R. Garcia-Molina, Water equivalent properties of materials commonly used in proton dosimetry, Applied Radiation and Isotopes, 83 (2014), pp. 122–127. [5] J. Demel, Operaˇcn´ı v´yzkum. http://kix.fsv.cvut.cz/~demel/ped/ov/ov.pdf. [Online; pˇr´ıstup 08-04-2014]. [6] N. Evans, Phys3016: lecture 28th february 2008. http://www.southampton.ac.uk/~evans/PHYS3017/Rel.pdf. [Online; pˇr´ıstup 08-04-2014]. [7] B. Gottschalk, On the scattering power of radiotherapy protons, Medical physics, 37 (2009), pp. 352–367. [8] B. Gottschalk, A. Koehler, R. Schneider, J. Sisterson, and M. Wagner, Multiple coulomb scattering of 160 mev protons, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 74 (1993), pp. 467–490. [9] M. Gupta, Calculation of radiation length in materials, tech. report, 2010. [10] V. L. Highland, Some practical remarks on multiple scattering, Nuclear Instruments and Methods, 129 (1975), pp. 497–499. [11] L. Hong, M. Goitein, M. Bucciolini, R. Comiskey, B. Gottschalk, S. Rosenthal, C. Serago, and M. Urie, A pencil beam algorithm for proton dose calculations, Physics in medicine and biology, 41 (1996), p. 1305. [12] E. Hussein, Handbook on Radiation Probing, Gauging, Imaging and Analysis: Volume II Applications and Design, Basics and techniques, Springer, 2003. ´ , V´ıcekriteri´ [13] J. Jablonsky aln´ı a c´ılov´e programov´ an´ı. http://nb.vse.cz/~JABLON/doc/vkr.pdf. [Online; pˇr´ıstup 08-04-2014]. [14] R. Kress, Numerical Analysis, Graduate Texts in Mathematics, Springer New York, 1998. [15] S. Lang and O. Riesterer, Modern techniques in radiation oncology. http://www.sps.ch/artikel/progresses/modern_techniques_in_radiation_oncology_36. [Online; pˇr´ıstup 08-04-2014]. `re, Theorie der streuung schneller geladener teilchen i. einzelstreuung am [16] v. G. Molie abgeschirmten coulomb-feld, Zeitschrift Naturforschung Teil A, 2 (1947), p. 133. [17]

, Theorie der streuung schneller geladener teilchen ii. mehrfach-und vielfachstreuung, Zeitschrift Naturforschung Teil A, 3 (1948), p. 78.

[18] H. Paganetti, Proton Therapy Physics, Series in Medical Physics and Biomedical Engineering, CRC Press/Taylor & Francis, 2012. [19] D. Pflugfelder, Risk-adapted optimization in intensity modulated proton therapy (impt), (2008).

18

[20] E. W. Weisstein, Parabolic cylinder function. http://mathworld.wolfram.com/ParabolicCylinderFunction.html. [Online; pˇr´ıstup 08-04-2014]. [21] R. Zhang and W. D. Newhauser, Calculation of water equivalent thickness of materials of arbitrary density, elemental composition and thickness in proton beam irradiation, Physics in medicine and biology, 54 (2009), p. 1383. [22] R. Zhang, P. J. Taddei, M. M. Fitzek, and W. D. Newhauser, Water equivalent thickness values of materials used in beams of protons, helium, carbon and iron ions, Physics in medicine and biology, 55 (2010), p. 2481.

19