Optimal control model for radiation therapy inverse planning applying the Boltzmann transport ...

20 Pages·2008·0.24 MB·English

View metadata, citation and similar papers at core.ac.uk brought to you by CORE
provided by Elsevier - Publisher Connector
Available online at www.sciencedirect.com
LinearAlgebraanditsApplications428(2008)1230–1249
www.elsevier.com/locate/laa
Optimal control model for radiation therapy inverse
planning applying the Boltzmann transport equation
J. Tervo a,∗, M. Vauhkonen b, E. Boman b
a DepartmentofMathematicsandStatistics,UniversityofKuopio,P.O.Box1627,FIN-70211,Kuopio,Finland
b DepartmentofPhysics,UniversityofKuopio,Kuopio,Finland
Received31July2006;accepted4March2007
Availableonline24March2007
SubmittedbyY.Censor
Abstract
We consider an inverse problem related to external radiation therapy treatment planning. The dose
calculation (the forward problem) is based on the Boltzmann Transport Equation (BTE) which models
exactlythetransportofchargedparticlesintissue.Theinverseplanning(theinverseproblem)isformulated
asanoptimalboundarycontrolproblem.Theoptimalcontrolvariableistheincoming(external)fluxand
theoutputisthedosedistributioninpatientdomain.Bothphysicalandbiologicalcostfunctionsaredefined
buthereweconcentrateonthephysicallybasedoptimization.Thediscretizationisdonebyfiniteelement
approximationsforwhichtheBTEisexpressedinitsvariationalform.Intheoptimizationprocessofthe
discretemodelweapplysocalledparametrization(whichmayessentiallydiminishthedecisionvariablesin
optimization)ofthematrixequations.Parametrizationisbasedoncertainlinearalgebraicdecompositions
ofmatrices.Onesimulationispresentedwhichshowthefunctionalityofthedevelopedmethods.
©2007ElsevierInc.Allrightsreserved.
Keywords: Boltzmanntransportequation;Radiotherapy;Optimalcontrol;Parametrization;Globaloptimization
1. Introduction
In the field of radiation therapy much research has been done to obtain new innovative
techniques for treating cancer patients with radiation. New equipments and new computational
approacheshavebeendevelopedbuttherestillisalackofautomatictreatmentplanningsystem.
Treatmentplanningisaninverseprobleminnature.Insolvingtheinverseproblemonealways
∗ Correspondingauthor.
E-mailaddress:Jouko.Tervo@uku.fi(J.Tervo).
0024-3795/$-seefrontmatter(2007ElsevierInc.Allrightsreserved.
doi:10.1016/j.laa.2007.03.003
J.Tervoetal. / LinearAlgebraanditsApplications428(2008)1230–1249 1231
needsadosecalculationmodel,thatis,aphysical“doselaw”whichdeterminesthedoseinthe
patientbasedontheknowledgeoftheradiationfluxincidentonthesurfaceofthepatientdomain.
Inthemathematicalterminology,thismeansthesolutionoftheforwardproblem.
For solving the forward problem one has developed and applied various dose calculation
models. Widely used models are variants of so-called pencil beam models which require the
applicationofadosedepositionkernel.Themostaccurate,physicallyrigorousandrobustmodels
arebasedontheBoltzmanntransportequation(BTE)[1,6,8,9,12,15].Transportequationbased
dosecalculationmodelsarevalidininhomogeneousmaterialandtheyalsotakerigorouslyinto
accountthescatteringeffects.Mainlyduetocomputationalproblems,BTEbasedmodelsarenot
yetextensivelyusedinsolvingtheforwardproblem.
Inoncologicalsocietythesolutionoftheinverseproblem,i.e.,theconceptof“thebestplan”
seems not to be absolutely unique. In some cases it may be difficult to measure the optimality
criteria.However,theprinciplesoftheoptimalplancanbefoundedeitheronthephysicalcriteria
oronthebiologicalcriteria.Nowadaysinpracticethecriteriaaremainlybasedonphysicalones.
Inthecaseofphysicalcriteriaonetriestocreateaprescribeddosedistributioninthetarget
(tumor) such that the dose in the healthy tissue is under certain prescribed value. In addition,
one takes special care of critical organs by demanding that the dose in vulnerable regions is
under certain prescribed level. Some additional constraints (e.g. dose volume constraints) may
be included. Using physical criteria the objective functions are relatively easy to construct and
a lot of various alternatives can be found in literature (e.g. [23,24] and the references therein).
Anotherpossibilityintheoptimizationapplyingphysicalcriteriaistoseekonlyfeasible(nearly
optimal)solutions[5,16].Recentlywehavedevelopedoptimizationalgorithmswhichusedirectly
themultileafcollimator(MLC)systemparametersinoptimization[16–19].Thecorresponding
deliveryquality(objective)functionsareusuallyhighlynonlinearandmultiextremal.Themodel
constraintsaremainlylinearinequalitiesandsorelativelysimple.
Intheoptimizationwhichisbasedonbiologicalcriteriaoneutilizesstatisticaldatacollected
fromclinicaltrials.Inthecaseofbiologicalcriteria,twodifferentprobabilities,namelythetumor
controlprobability(TCP)andthenormaltissuecomplicationprobability(NTCP)areemployedin
thetreatmentplanning[3,4,24].TheaimistoprovideashighTCPaspossibleandatthesametime
maintaintheNTCPatanacceptablelevelincertainorgansofthepatient.Thebiologicalobjective
functions are more difficult to construct and as a rule one can say that the biological objective
functionsaremathematicallymorecomplextohandlethanthephysicalobjectivefunctions.
Recentlywehavedevelopedanoptimalcontrolapproachforsolvingtheoptimizationproblem
[20].InthisapproachtheBTEhasbeenutilizedandtheoptimizationproblemissolvedbythe
boundarycontroloptimizationtheory.Inthispaperweintroducearadiationtreatmentplanning
approachthatisbasedontheBoltzmanntransportequationandoptimalcontrol.Anovellinear
algebraicparametrizationapproachisutilizedforthediscretizedmodelinordertoreducethenum-
berofoptimizedparametersintheoptimizationprocedure.Solvingeffectivelylargedimensional
linearsystemsandthecalculationofthesingularvaluedecompositionsforthemaredifficulties
inourapproach.Thecontrolparametersinthisstudyaretheintensitiesoftheappliedradiation
fields, though it would be fairly straightforward to formulate the problem also for the case of
MLC.Inthatcase,however,theobjectivefunctionsbecomehighlynonlinearandthenumberof
decisionparametersincreases.Numericalresultsfortwo-dimensionalpatientbodyareshownin
thepaper.Innumericalconsiderationswehaveconcentratedontheuseofphysicaloptimization
criteriabuttheoptimizationschemeapplyingthebiologicaloptimizationcriteriaisalsodefined
inthepaper.Theoptimizationproblemsforrealistic(3-dimensionalpatientbody)casesarevery
largedimensionalandtheyarestillcomputationallytooexpensive.
1232 J.Tervoetal. / LinearAlgebraanditsApplications428(2008)1230–1249
2. Calculationofdose
2.1. Boltzmanntransportequationmodel
PhysicallytheBoltzmanntransportequationisbasedontheparticleequilibriumininfinitesimal
small voxels of tissue. The modeling requires the knowledge of the differential and total cross
sectionswhichdependonthenatureofinteractionsofparticles.Thetotaldoseisobtainedasa
superpositionofdifferentfieldsaswillbedescribedindetailbelow(Section2.3).
We consider only the stationary transport of particles. Assume that ψ =ψ (x,E,(cid:2)),j =
j j
1,2,3 are the phase space densities for photons, electrons and positrons, respectively. x =
(x ,x ,x )isthepointinthepatientdomainV ⊂R3.(cid:2)isapointontheunitsphereS ofR3.
1 2 3
Inthecasewhereelasticcollision,inelasticcollisionandbremsstrahlungaretakenintoaccount
themodelconsistsofthefollowingcoupledsystemofpartialdifferential–integralequations
(cid:2)·∇ψ +K (ψ ,ψ ,ψ )=Q (x,E,(cid:2)),
1 1 1 2 3 1
(cid:2)·∇ψ +K (ψ ,ψ ,ψ )=Q (x,E,(cid:2)), (2.1)
2 2 1 2 3 2
(cid:2)·∇ψ +K (ψ ,ψ ,ψ )=Q (x,E,(cid:2)).
3 3 1 2 3 3
Themodel(2.1)islinear.Itneglectssomenonlinearinteractionswhicharenotessentialindose
calculation.ThefunctionsK (ψ ,ψ ,ψ ),j =1,2,3arecollisiontermsresultingfromdifferent
j 1 2 3
kindofinteractions(mentionedabove).K ,K ,K arelinearfunctionsofψ:=(ψ ,ψ ,ψ ).The
1 2 3 1 2 3
interactionscanbedescribedbytheintegrals(j =1,2,3)
(cid:3) (cid:3)
(cid:2)3 (cid:2)3
K ψ = σ (x,E(cid:5),E,(cid:2)(cid:5),(cid:2))ψ (x,E(cid:5),(cid:2)(cid:5))d(cid:2)(cid:5)dE(cid:5)− (cid:3) (x,E)ψ .
j k,j k k,j k
k=1 I S k=1
Aboveσ (x,E(cid:5),E,(cid:2)(cid:5),(cid:2))aredifferentialcrosssectionsand(cid:3) (x,E)aretotalcrosssections.
k,j k,j
Iftheinteractionfromparticlej toparticlekisnotpossible,σ =0.Fortotalcrosssections
k,j
(cid:3) =0, j =/ k. (2.2)
k,j
Finally,Qj(x,E(cid:4),(cid:2))arethesourceterms.Thesemaydescribesourcesintissue(interiortherapy).
The integral f((cid:2))d(cid:2) denotes the standard surface integral. We denote by I the energy
S
interval,say[E ,E ].Thesolutionψisdefinedinthe6-dimensionalstatespaceG:=V ×I ×S.
1 2
Inthefollowingwedenotethesystem(2.1)moresimplyby
((cid:2)·∇+K)ψ =Q, (2.3)
where
(cid:2)·∇ψ =((cid:2)·∇ψ ,(cid:2)·∇ψ ,(cid:2)·∇ψ )
1 2 3
and Kψ =(K ψ,K ψ,K ψ), Q=(Q ,Q ,Q ). We consider exterior radiation therapy.
1 2 3 1 2 3
Hence we have Q=0 and the boundary condition is used to model the extrageneous particle
fluxes.
2.2. Photoninflow.Theboundarycondition
We consider the photon inflow in the stationary case. This corresponds the exterior photon
radiation. Other modalities are similarly considered. We assume that the boundary (cid:2)V is a
Lipschitz-boundary. Then the outward normal ν(x) exists and is it continuous on (cid:2)V in spite
ofpossiblyasetwithsurfacemeasurezero.
J.Tervoetal. / LinearAlgebraanditsApplications428(2008)1230–1249 1233
Totakeintoaccounttheincomingextrageneousfluxwemustputsome(boundary)conditions
forthesolution.Intheexteriortherapythetypicalconditionforthesolutionψ isoftheform
ψ (x,E,(cid:2))=ψ (x,E,(cid:2))=0 for(x,E,(cid:2))∈(cid:2)V ×I ×S
2 3
suchthat(cid:2)·ν(x)<0,
ψ (x,E,(cid:2))=u(x,E,(cid:2)) for(x,E,(cid:2))∈(cid:2)V ×I ×S
1
suchthat(cid:2)·ν(x)<0. (2.4)
uisthephotonfluxdensityincidenton(cid:2)V.Weassumethatu∈L ((cid:2)V ×I ×S).Thecondition
2
ψ =ufor(cid:2)·ν(x)<0,x ∈(cid:2)V meansthatthebeam(thefluxu)isinwardlyonthepatch(cid:2)V
1
andtheconditionψ =0,j =2,3for(cid:2)·ν(x)<0meansthatnootherparticlesgenerateinward
j
fluxes.
Remark1. Fortheexteriorelectronradiationtherolesofψ andψ arechangedintheboundary
1 2
condition.Otherwisetheboundaryconditionissameasin(2.4).
2.3. Thedefinitionofdose
In practical situations the total dose distribution is computed as follows. Let the incoming
(initial)fluxdensityofthelthfieldS beu.Weassumethatu ∈L ((cid:4) ×I ×S),where(cid:4) is
l l l 2 l l
apatchof(cid:2)V throughwhichtheradiationisenteringintodomainV.Letψl =(ψl,ψl,ψl)be
1 2 3
thefluxdensitycorrespondingtothefieldS,thatisψl isthesolutionoftheequation
l
((cid:2)·∇+K)ψl =0 (2.5)
withtheboundarycondition
ψl(x,E,(cid:2))=ψl(x,E,(cid:2))=0 for(x,E,(cid:2))∈(cid:2)V ×I ×S
2 3
suchthat(cid:2)·ν(x)<0,
ψl(x,E,(cid:2))=0, for(x,E,(cid:2))∈((cid:2)V \(cid:4))×I ×S
1 l
suchthat(cid:2)·ν(x)<0, (2.6)
ψl(x,E,(cid:2))=u(x,E,(cid:2)) for(x,E,(cid:2))∈(cid:4) ×I ×S
1 l l
suchthat(cid:2)·ν(x)<0.
Above
(cid:2)·∇ψl:=((cid:2)·∇ψl,(cid:2)·∇ψl,(cid:2)·∇ψl), Kψl:=(K ψl,K ψl,K ψl).
1 2 3 1 2 3
ThedosecontributionD(x)fromthefieldS atapointxofthepatientdomainV isobtained
l l
fromthe(measurement)integral
(cid:3) (cid:3)
(cid:2)3
D(x)= κ (x,E)ψl(x,E,(cid:2))dEd(cid:2), (2.7)
l j j
j=1 S I
where κ (x,E) are known (stopping power) factors (κ (x,E)=0). The total dose is obtained
j 1
from
1234 J.Tervoetal. / LinearAlgebraanditsApplications428(2008)1230–1249
(cid:3) (cid:3)
(cid:2)L (cid:2)L (cid:2)3
D(x)= D(x)= κ (x,E)ψl(x,E,(cid:2))dEd(cid:2). (2.8)
l j j
l=1 l=1j=1 S I
Thecomputationoftotaldosecanalsobeformulatedasfollows.Defineu∈L ((cid:2)V ×I ×S)
2
suchthat
(cid:2)L
u= uχ, (2.9)
l l
l=1
whereχ :(cid:2)V ×I ×S →Rarethecharacteristicfunctionsof(cid:4) ×I ×S.Letψ =(ψ ,ψ ,ψ )
l l 1 2 3
bethesolutionoftheproblem
((cid:2)·∇+K)ψ =0 (2.10)
withtheboundarycondition
ψ (x,E,(cid:2))=ψ (x,E,(cid:2))=0 for(x,E,(cid:2))∈(cid:2)V ×I ×S
2 3
suchthat(cid:2)·ν(x)<0,
ψ (x,E,(cid:2))=u, for(x,E,(cid:2))∈(cid:2)V ×I ×S
1
suchthat(cid:2)·ν(x)<0, (2.11)
whereuisdefinedby(2.9).
Thesolutionoftheproblem(2.10)–(2.11)is
(cid:5) (cid:6)
(cid:2)L (cid:2)L (cid:2)L
ψ = ψl, ψl, ψl , (2.12)
1 2 3
l=1 l=1 l=1
whereψl(l =1,...,L)arethesolutionsof(2.5)–(2.6).Theprooffollowsfromtheuniqueness
ofsolutions(2.3)–(2.4)whichwewillformulateinSection2.4.Nowthetotaldoseis
(cid:3) (cid:3)
(cid:2)3
D(x)= κ (x,E)ψ (x,E,(cid:2))dEd(cid:2), (2.13)
j j
j=1 S I
whereψ isthesolutionof(2.10)–(2.11).Wefindthat(2.10)–(2.11)isexactlytheproblem(2.3)–
(2.4)withu∈L ((cid:2)V ×I ×S)givenby(2.9).
2
Remark2. Theenergydepositioninapatientisduetothechargedparticles.Henceinpractice
wecanassumethatκ (x,E)=0.Moregenerallythefactorsκ maydependalsoon(cid:2).
1 j
2.4. Variationalformofequations
Let G be as above and let L (G) be the Lebesgue space of (real valued) square integrable
2
functionsonGwiththeusualinnerproduct.Furthermore,let
(cid:7)
L ((cid:2)V ×I ×S,|(cid:2)·ν|dσdEd(cid:2)) := g :(cid:2)V ×I ×S →R|gismeasurableand
2
(cid:8)
(cid:3) (cid:3) (cid:3)
|(cid:2)·ν|g2dσdEd(cid:2)<∞ ,
S I (cid:2)V
J.Tervoetal. / LinearAlgebraanditsApplications428(2008)1230–1249 1235
whereσ isthesurfacemeasureon(cid:2)V.ThenL ((cid:2)V ×I ×S,|(cid:2)·ν|dσdEd(cid:2))isaHilbertspace
2
equippedwiththeinnerproduct
(cid:3) (cid:3) (cid:3)
(cid:9)g1,g2(cid:10)L2((cid:2)V×I×S,|(cid:2)·ν|dσdEd(cid:2)):= |(cid:2)·ν|g1g2dσdEd(cid:2).
S I (cid:2)V
Denote
D(V ×I ×S):={f| |f ∈C∞(R3×R×S)}.
V×I×S
LetH bethecompletionofD(V ×I ×S)withrespecttotheinnerproduct
1
(cid:9)f1,f2(cid:10)H1:=(cid:9)f1,f2(cid:10)L2(G)+(cid:9)f1,f2(cid:10)L2((cid:2)V×I×S,|(cid:2)·ν|dσdEd(cid:2)).
Furthermore,letH bethecompletionofD(V ×I ×S)withrespecttotheinnerproduct
2
(cid:9)f ,f (cid:10) =(cid:9)f ,f (cid:10) +(cid:9)(cid:2)·∇f ,(cid:2)·∇f (cid:10) .
1 2 H2 1 2 L2(G) 1 2 L2(G)
Finally,letH:=H ∩H isequippedwiththestandardHilbertspaceinnerproduct
1 2
(cid:9)f ,f (cid:10) =(cid:9)f ,f (cid:10) +(cid:9)f ,f (cid:10) .
1 2 H 1 2 H1 1 2 H2
Foranyf ∈H1 thetracef|(cid:2)V×I×S) iswelldefinedinthesensethatthereexistsasequence
{f }⊂D(V ×I ×S)suchthat
n
fn →f|(cid:2)V×I×S inL2((cid:2)V ×I ×S,|(cid:2)·ν|dσdEd(cid:2)).
One knows that for each f ∈H2 the restriction f|(cid:2)V×I×S ∈L2(K) where K is a compact
subsetof
{(x,E,(cid:2))∈(cid:2)V ×I ×S||(cid:2)·ν(x)|>0}
but f|(cid:2)V×I×S is not necessarily in L2((cid:2)V ×I ×S,|(cid:2)·ν|dσdEd(cid:2)) [8, pp. 220–221]. For
ψ,v ∈H theGreen’sformula
(cid:3) (cid:3) (cid:3)
(cid:9)ψ,(cid:2)·∇v(cid:10) +(cid:9)(cid:2)·∇ψ,v(cid:10) = ((cid:2)·ν)ψvdσdEd(cid:2) (2.14)
L2(G) L2(G)
S I (cid:2)V
isvalid[8,p.225].
IntheproductspacesH3,i =1,2weusetheusualinnerproducts
i
(cid:2)3
(cid:9)f,h(cid:10) = (cid:9)f ,h (cid:10)
H3 j j Hi
i
j=1
for f =(f ,f ,f ),h=(h ,h ,h )∈H3. In the similar way we define the inner product in
1 2 3 1 2 3 i
H3.
Inthefollowingthesubscript“−”referstothenegativepartofafunctionandthesubscript“+”
referstothepositivepartofafunction.Thevariationalformulationoftheproblem(2.3)–(2.4)is
givenby[21].
Theorem1. Assumethat
1. (cid:3)k,j ∈L∞(V ×I), (2.15)
2. σ ∈C(V ×I2×S2), (2.16)
k,j
3. Q ∈L (G),u∈L ((cid:2)V ×I ×S). (2.17)
j 2 2
1236 J.Tervoetal. / LinearAlgebraanditsApplications428(2008)1230–1249
Then the variational form of the equation (2.10) with the stated boundary condition (2.11) is
givenby
B(ψ,v)=F(v), v ∈H3, (2.18)
whereB(·,·):H3×H3 (cid:12)→Risthebilinearform
(cid:3) (cid:3) (cid:3)
(cid:2)3
B(ψ,v)=−(cid:9)ψ,(cid:2)·∇v(cid:10)L2(G)3 +j=1 S I (cid:2)V((cid:2)·ν)+ψjvjdσdEd(cid:2)+(cid:9)Kψ,v(cid:10)L2(G)3
(2.19)
and
(cid:3) (cid:3) (cid:3)
F(v)=(cid:9)Q,v(cid:10)L2(G)3 + S I (cid:2)V((cid:2)·ν)−uv1dσdEd(cid:2). (2.20)
Weformulatethefollowinglemma[21]:
Lemma1. Assumethat
1. (cid:3)k,j ∈L∞(V ×I), (2.21)
2. σ ∈C(V ×I2×S2). (2.22)
k,j
3.Thereexitsκ >0suchthatforψ ∈H3
(cid:9)Kψ,ψ(cid:10) (cid:2)κ(cid:13)ψ(cid:13)2 . (2.23)
L2(G)3 L2(G)3
ThenthebilinearformB(ψ,v)satisfies
B(ψ,v)(cid:3)C(cid:13)ψ(cid:13) (cid:13)v(cid:13) (boundedness) (2.24)
H3 H3
1
forψ,v ∈H3and
B(ψ,ψ)(cid:2)c(cid:13)ψ(cid:13)2 (H3−coercitivity) (2.25)
H3 1
1
forψ ∈H3.
Inaddition,F ∈(H3)∗(herethesuperscript“∗”referstotheadjointspace)andthereexists
1
C >0suchthat
(cid:13)ψ(cid:13) (cid:3)C(cid:13)F(cid:13), (2.26)
H3
1
where
(cid:9)
(cid:3) (cid:3) (cid:3)
(cid:13)F(cid:13):=(cid:13)Q(cid:13)L2(G)3 + S I (cid:2)V |((cid:2)·ν)−|u2dσdEd(cid:2).
In [2,21] we have proved the following sufficient algebraic criterion for the coercitivity
assumption(2.23)
J.Tervoetal. / LinearAlgebraanditsApplications428(2008)1230–1249 1237
Lemma2. Supposethat(cid:3)j,j ∈L∞(V ×I)andσk,j ∈C(V ×I2×S2).Furthermore,suppose
that(2.2)isvalidandthatthereexistsα >0suchthatalmosteverywhere(x,E,(cid:2))∈G
(cid:3) (cid:3)
(cid:2)3
(cid:3) (x,E,(cid:2))− σ (x,E(cid:5),E,(cid:2)(cid:5),(cid:2))dE(cid:5) d(cid:2)(cid:5) (cid:2)α (2.27)
j,j k,j
S I k=1
and
(cid:3) (cid:3)
(cid:2)3
(cid:3) (x,E,(cid:2))− σ (x,E,E(cid:5),(cid:2),(cid:2)(cid:5))dE(cid:5)d(cid:2)(cid:5) (cid:2)α (2.28)
j,j j,k
S I k=1
forj =1,2,3.Thentheassumption
(cid:9)Kψ,ψ(cid:10) (cid:2)α(cid:13)ψ(cid:13)2 (2.29)
L2(G)3 L2(G)3
isvalid.
Aphysicalbackgroundcanbeeasilyfoundforthecondition(2.27),wheretheintegrationsover
energyandangulardomainsbasicallydescribethetotalscatteringcrosssections.Thecondition
(2.27)statesthatthesumofthesescatteringcrosssectionscannotbelargerthanthetotalcross
section itself. This is always true, because the total cross section is the sum of scattering and
absorptioncrosssections.Thephysicalmeaningofthesecondcondition(2.28)isnotsoclear.It
statesthatthoseintegrated‘inverse’scatteringcrosssections,whichchangeparticlesfromj tok
cannotbelargerthanthetotalcrosssectionoftheparticlej.However,itseemsthatthissecond
condition(2.28)isalsophysicallyrelevant.
By the boundedness (2.24) the bilinear form B(·,·) can be extended on H3×H3 (here we
1
denotetheextensionagainbyB(·,·)).Theextensionsatisfiestheestimates
|B(ψ,v)|(cid:3)C(cid:13)ψ(cid:13) (cid:13)v(cid:13) , ψ ∈H3, v ∈H3 (2.30)
H3 H3 1
1
and
B(ψ,ψ)(cid:2)c(cid:13)ψ(cid:13)2 , ψ ∈H3. (2.31)
H3
1
Weformulatethefollowingexistenceresultofsolutions:
Theorem2. Suppose that the assumptions of Theorem 1 are valid and that (2.2) and (2.27)–
(2.28)hold.Thenthevariationalequation
B(ψ,v)=F(v), v ∈H3 (2.32)
hasoneandonlyonesolutionψ ∈H3.
Inaddition,F ∈(H3)∗and
1
1
(cid:13)ψ(cid:13) (cid:3) (cid:13)F(cid:13), (2.33)
H13 c
where
(cid:9)
(cid:3) (cid:3) (cid:3)
(cid:13)F(cid:13)(cid:3)(cid:13)Q(cid:13)L2(G)3 + S I (cid:2)V |((cid:2)·ν)−|u2dσdEd(cid:2). (2.34)
Proof. The proof follows from estimates (2.30), (2.31) and the generalized Lax–Milgram
Theorem[22,p.403].Weomitthedetailshere. (cid:4)
1238 J.Tervoetal. / LinearAlgebraanditsApplications428(2008)1230–1249
Theorem2givestheexistenceofweaksolutionsfortheproblem(2.10)–(2.11).
3. Optimalcontrolprobleminradiationtherapy
3.1. Inverseproblembasedonphysicalobjectives
ThepatientdomainV ⊂R3containstumor’s(target’s)regionT,criticalorgans’regionCand
normaltissue’sregionNandsowehavetheunionV =T∪C∪N.LetusassumethatwehaveL
fieldsS,l =1,...,L.Thismeansthatgantry,couchandcollimatoranglesaredeterminedand
l
thewholetreatmentcontainsLdifferentanglesettings.
Abovewehavefoundthatthedosecanbeobtainedfromthefunctional
(cid:3) (cid:3)
(cid:2)3
D(x)= κ (x,E)ψ (x,E,(cid:2))dEd(cid:2), (3.1)
j j
j=1 S I
whereψ isthesolutionof(2.10)–(2.11).In(3.1)κ (x,E)areknownfactors(so-calledstopping
j
powers).Weassumethatκj ∈L∞(V ×I).Consideringtheexternaltherapywehavenointernal
sources.SoQ =0foreachfieldS.Applyingtheaboveconcepts,theinverseradiationtreatment
l l
planningproblemstates:
SupposethatD istheprescribed(uniform)doseintumorTandthatD andD aretheupper
0 C N
bounds of dose in the critical organ C and in the normal tissue N, respectively. Furthermore,
supposethatthenumberLandthegantry,couchandcollimatoranglesα,β andθ offieldsS
l l l l
aregiven.
Determinetheincomingfluxu∈L ((cid:2)V ×I ×S)suchthat
2
D(x)=D , x ∈T,
0
D(x)(cid:3)D , x ∈C,
C
D(x)(cid:3)D , x ∈N. (3.2)
N
andthat
u(cid:2)0. (3.3)
Tumor, critical organ and normal tissue may be divided into many separate parts and dose
limitsmaybedifferentintheseparts.
Besidestherequirements(3.2)oneoftendemandsthatsocalleddosevolumeconstraints[3,23]
are fulfilled. Dose volume constraints may be necessary for certain structures, for example for
criticalorgans.WedescribethisrequirementshortlyforcriticalorganC.Letv(D)bethevolume
fractionofCthatreceivesadosegreaterthanD.Wedemandthat
v(D)(cid:3)v whenD (cid:2)d , (3.4)
0 C
where v is a given volume fraction and d is a given dose. Since the function v =v(D) is a
0 C
decreasingfunctionofDthecondition(3.4)isequivalentto
v(d )(cid:3)v . (3.5)
C 0
Thecondition(3.5)meansthat
μ({x ∈C|D(x)(cid:2)d })
C (cid:3)v , (3.6)
0
μ(C)
J.Tervoetal. / LinearAlgebraanditsApplications428(2008)1230–1249 1239
(cid:10)
whereμistheLebesguemeasure.LetH :R→RbetheHeavisidefunctionH(x)= 1, x(cid:2)0.
0, x<0
Thenthedosevolumeconstraint(3.6)canbeexpressedas
(cid:3)
1
H(D(x)−d )dx (cid:3)v . (3.7)
C 0
μ(C)
C
In computer programs the Heaviside function H can be accurately replaced by the error func-
(cid:4)
tion erf (x)= √1 x e−s2/(cid:10)2ds. Furthermore, the requirement (3.6) can be approximated
(cid:10) π(cid:10) −∞
accuratelybytherequirement
(cid:2)
1
erf (D(x )−d )(cid:3)v , (3.8)
|I | (cid:10) p C 0
C p∈IC
whereI isaselectedfinitesetofC(forexample,thesetofnodesinfiniteelementcomputations).
C
Whenthedosevolumeconstraintsaretakenintoaccountrequirementslike(3.7)(or3.8))areadded
totherequirements(3.2).
3.2. Optimalcontrolproblem
Toclarifytheu-dependenceofvariableswedenoteψ =ψ(u),ifneeded.LetF :L ((cid:2)V ×
2
I ×S)(cid:12)→(H3)∗betheoperatordefinedby
(cid:3) (cid:3) (cid:3)
(Fu)(v)= ((cid:2)·ν)−uv1dσdEd(cid:2).
S I (cid:2)V
NotethatF(v)=(Fu)(v)(orF =Fu).Usingthesenotationsψ =ψ(u)satisfiesthevariational
equation
B(ψ(u),v)=(Fu)(v), v ∈H3. (3.9)
LetL:L (G)3 →L (V)bethefunctional
2 2
(cid:3) (cid:3)
(cid:2)3
(Lψ)(x)= κ (x,E)ψ (x,E,(cid:2))dEd(cid:2). (3.10)
j j
j=1 S I
Thenwehaveforthedose
D =Lψ =Lψ(u). (3.11)
ForlaterneedswenoticethattheadjointL∗ :L (V)→L (G)3ofLisgivenby
2 2
L∗v =(κ v,κ v,κ v), v ∈L (V). (3.12)
1 2 3 2
The dose D(x) must be as near as possible to the described dose D in the tumor and the
0
upperboundsofdoseincriticalorgansandnormaltissuemaynotbeviolated.Hencewetryto
optimizetheincomingfluxusothatthisholds.Theconcreteimplementationofthisleadstothe
followingkindofoptimizationproblems,forexample.
Defineacostfunctionalby
J(u)=c (cid:13)D −Lψ(u)(cid:13)2
1 0 L2(T)
+c2(cid:13)(DC−Lψ(u))−(cid:13)2L2(C)+c3(cid:13)(DN−Lψ(u))−(cid:13)2L2(N)
+c4(cid:13)(u)−(cid:13)2L2((cid:2)V×I×S)+c5(cid:13)u(cid:13)2L2((cid:2)V×I×S), (3.13)