Geant4不同物理模型對放療質子束模擬的影響

林 輝，謝 聰，張擁軍,2，熊楨宇，吳東升,*，曹瑞芬，FDS團隊

(1.合肥工業大學 電子科學與應用物理學院，安徽 合肥 230009；2.蚌埠醫學院 公共課程部，安徽 蚌埠 233000；3.中國科學院 核能安全技術研究所，安徽 合肥 230031)

?

Geant4不同物理模型對放療質子束模擬的影響

林 輝1，謝 聰1，張擁軍1,2，熊楨宇1，吳東升1,*，曹瑞芬3，FDS團隊3

(1.合肥工業大學 電子科學與應用物理學院，安徽 合肥 230009；2.蚌埠醫學院 公共課程部，安徽 蚌埠 233000；3.中國科學院 核能安全技術研究所，安徽 合肥 230031)

Monte Carlo方法常用于質子束劑量模擬以優化質子束臨床治療，合適的物理模型選擇關系到Monte Carlo模擬結果的可靠性及模擬效率。本工作選取9個可用于質子模擬的Geant4物理模型，計算比較了質子束在多種介質中的劑量及次級粒子產額。結果表明：標準電磁物理模型、低能電磁Penelope模型和Livermore模型雖在模擬質子整體劑量分布上可接受，但在微觀上缺失重粒子生成。參數化驅動模型LHEP的模擬時間最短,但與QGSP相同，均不能產生復雜重離子。QGSP_BIC_EMY模型較其他模型明顯低估Bragg峰-入口劑量比和絕對劑量。QGSP_BERT模型、QGSP_BIC模型和FTFP_BERT模型較適合放療質子束模擬。

物理模型；質子放療；微觀產額；劑量特征量

質子放療是一新興的腫瘤物理治療技術。質子劑量曲線具有尖銳Bragg峰，在穿越路徑時僅釋放出少許能量，而達到治療深度時釋放出大量能量[1]。放療質子的能量范圍通常為60～250 MeV，可治療從淺表層腫瘤(如眼部腫瘤)到深度腫瘤[2]。

蒙特卡羅(MC)方法通過建立求解問題的概率模型，對模型的多次抽樣給出計算問題的近似統計解[3]。MC程序可利用機器隨機數模擬粒子輸運發生的復雜物理過程，給出粒子出射的角分布特征[4]、非均勻介質中精確的能量分布[5]以及微觀分析關鍵物理作用的權重比例[6]等。但MC方法不同于實驗測量，需基于預先知道的理論模型對粒子輸運作用和方向等進行判斷。

目前模擬質子輸運的程序有MCNPX、Fluka、Geant4[7]等。其中開源程序平臺Geant4包含了描述粒子與物質相互作用的多種物理模型，其面向對象的模塊化編程設計增強了程序的可移植性。關鍵的問題是使用者需了解哪些物理模型適合所研究問題[8]。錯失粒子類型或物理作用均可造成結果偏差，但過于細致的物理模型又會造成模擬速度的降低。文獻[2,5,6,8]比較了不同物理模型對質子總沉積劑量的影響，但對不同物理模型對質子微觀產額的影響缺乏研究。

本工作通過開發Geant4用戶程序，基于Geant4類庫提供的9種物理模型模擬放療質子束在多種介質中的劑量及次級粒子微觀產額，擬找到平衡精度和速度適合模擬治療質子束能量范圍的物理模型。

1 材料和方法

1.1 程序發展

本用戶程序基于Geant4程序例extendedelectromagneticTestEm3進行了二次開發，借助其內外環套(沿y-z面)的層狀(沿x方向)幾何設計，模擬質子束照射模體(圖1)。Geant4程序sourcephysics_listslists目錄下提供了常用的模型物理文件，特別是一些字符型模型可進行強子輸運模擬。若用戶需要還可制作自己的模型物理文件，所基于的初級模塊可在sourcephysics_listsuilderssrc目錄下查找。本工作采用Geant4_9_5_p02版本。

圖1 照射模體Fig.1 Configuration of irradiation

1.2 Geant4物理模型

Geant4程序平臺提供了多種粒子和物理作用模型。Geant4提供了玻色子、輕子、介子、重子、短壽命粒子和離子6類粒子的定義。程序中需考察的初級粒子和次級粒子類型必須在PhysicsList中進行ConstructParticle()構建。Geant4包含的物理作用有電磁作用、強子相互作用、輸運、衰變、光學作用、光輕子-強子、特定體積參數化等。程序中需用到的作用類型必須在PhysicsList的ConstructProcess( )中構建。

Geant4提供的電磁作用模型有標準電磁模型(standard EM，約1 keV～1 TeV)、低能電磁模型(250 eV～1 PeV)和光量子模型。標準電磁模型提供了G4EmStandardPhysics、G4EmStandardPhysics_option1、G4EmStandard-Physics_option2和G4EmStandardPhysics_option3 4個物理清單。G4EmStandardPhysics_option3考慮了質子的各種電磁相互作用，適合醫學物理模擬。低能電磁模型LowE_Penelope和LowE_Livermore對于γ和e-/e+(電子/正電子)輸運在0～1 GeV之間使用自己的模型，但在1 GeV～10 TeV之間仍使用standard EM。雖然標準和低能電磁模型未直接模擬強子作用過程，但均包含了一些質子電離處理的替代模型，并加入次級粒子產生處理[8]。

Geant4提供以下強子相互作用模型。

1) 參數化驅動模型(PDM，parametrization driven model)。模型中僅強子-核子碰撞被詳細模擬，其他發生在核子內的相互作用是通過產生額外的強子，并當作次級粒子進行模擬。而額外的強子數目、類型和分布是通過擬合實驗數據獲取的函數決定的。因此，此模擬需依賴大量的參數獲取，以保證模型的合理性。目前有兩套PDM可用于模擬高能強子-核子相互作用：低能模型(1～25 GeV)和高能模型(25 GeV～10 TeV)。

2) 級聯模型(CM，cascade model)。模擬了基本粒子和次級粒子的級聯輸運過程，包括BIC(binary cascade model(約<10 GeV))和BERT(Bertini cascade model(約<10 GeV))。但級聯模型不適合低能，因為低能涉及蒸發和預平衡衰變。僅當粒子的動能高于閾值(75 MeV)時才能使用級聯模型。

3) 預復合核模型(PRECO，pre-compound model)。預復合核模型是強子動力學模型的擴展，它將強子動力學模型適用范圍擴展到核-核非彈性碰撞。考慮核反應的預復合核階段直到非平衡態。從激發核里釋放核碎片或光子的模擬采用平衡模型。PRECO適用于0～150 MeV范圍內的質子、中子模擬。

此外，適合高能強子模擬的模型還有QGS(quark gluon string model(約>20 GeV))和FTF(Fritiof model(約>10 GeV))，它們與低能模型PRECO結合成為QGSP和FTFP才適合治療質子束模擬。

本工作考察了以下物理模型對放療質子劑量和微觀產額的影響。

模型1，G4EmStandardPhysics_option3；模型2，LowE_Penelope；模型3，LowE_Livermore；模型4，LHEP(低能與高能PDM)+模型1；模型5，QGSP(QGS+PRECO)+模型1；模型6，QGSP_BERT(QGS+PRECO+BERT)+模型1；模型7，QGSP_BIC(QGS+PRECO+BIC)+模型1；模型8，FTFP_BERT (FTF+PRECO+BERT)+模型1；模型9，QGSP_BIC_EMY(QGSP_BIC模型+EMY處理)+模型1。

其中，模型1～3為電磁模型。為考察強子模型的影響，模型4～9在含有標準電磁模型選項3(即模型1)之外又增添了強子模型處理。模型9為Geant4advancedhadrontherapy例提供的物理模型，在使用QGSP_BIC模型之外添加了G4EmExtraPhysics( )、G4Hadron-ElasticPhysics( )、G4QStoppingPhysics( )、G4IonBinaryCascadePhysics( )、G4Neutron-TrackingCut()處理。除模型1外，其他模型均采用了放射活度衰變離子物理處理。

1.3 模擬幾何

由于質子側向散射較小，一般更關注其劑量曲線深度方向的分布。本模型采用一維長方體模型，模體大小為xcm×10 cm×10 cm，x為照射方向，可根據質子束能量調節(圖1)。計數探測器沿x方向(即照射中心軸)平行排列，在垂直于照射軸面上的大小為5 cm×5 cm，沿照射方向根據部位不同，分別取非關鍵深度1 mm以及入口和Bragg深度處1 mm、0.1 mm或0.01 mm。照射野為5 cm×5 cm，垂直入射。各粒子的截斷距離取1 mm。Geant4中截斷距離對應了產生δ電子和軔致輻射光子的射程閾值。Geant4中過小的射程截斷會顯著增加模擬時間和計算結果的漲落[9]。Geant4advancedhadrontherapy例中相對于0.2 mm的網格厚度也采用了1 mm的粒子射程截斷。為平衡精度和速度的矛盾，模擬前期進行了截斷值選擇的驗證工作，發現1 mm差異可接受。模體中分別充滿均勻介質水、脂肪、肺、肌肉和骨(密度分別為1、0.99、0.30、1.04、1.49 g/cm3，數據來自于advancedhuman_phantom例)。每種情況跟蹤2×105個粒子，數據不確定度<1%。整個計算是在2.67 GHz CPU/2G內存的計算機上進行。

1.4 模型校驗

2 結果與討論

2.1 PDD與次級粒子產額

為驗證模擬的可靠性，模擬的PDD與文獻[10]數據進行了對比(圖2)。從圖2可見，除物理模型9和文獻相差較大(約3.5%)外，其他模型幾乎完全一致(入口處(<5 cm)約1.0%，其他處約0.2%)。表1列出了各模型的次級粒子產額，其中模型1～3僅產生γ和e-，無重子產生，表明標準電磁模型和低能電磁模型中植入的替代強子輸運模型較簡單，未考慮重子生成。雖然三者的γ和e-產額有少許差異，但較相近。模型4和5可產生重子n、H+、2H、3H、α和輕子e+，且產額較接近，但不能產生一些重離子(如Li+、Be2+、B3+、C4+、N3-、O2-等及其同位素)。模型6～9可產生重離子，說明級聯模型BIC和BERT加入了重離子產生模擬。其中模型6的n、H+產額較模型7高20%～30%，但2H、3H、α產額較模型7降低了30%～50%。模型6和模型8的一致，這是因為在模擬的能量范圍內兩者均使用了BERT和PRECO模型。模型9的H+產額較模型7的增加約38%，但重離子產額增加約79%，反映EMY模型的影響。從模擬時間看，采用簡單參數化處理的LHEP模型的時間最短，QGSP_BIC_EMY時間最長(LHEP的170%)，其他模型時間相仿(約LHEP的120%)。

圖2 模擬與文獻的PDD對比Fig.2 Comparison of simulated and published PDDs

表1 高斯能譜次級粒子產額比較Table 1 Comparison of secondary yield for Gaussian spectrum

2.2 質子入射能量影響

表2列出了60 MeV和200 MeV的次級粒子產額，其總體特點與高斯能譜類似，但具體產額比有些差異，如n、H+產額方面，對于60 MeV，模型6較模型7高30%～50%，對于200 MeV僅高約20%；2H和3H產額方面，對于60 MeV，模型6較模型7降低了80%和5%，對于200 MeV降低了30%和27%，而α產額降低差不多(約45%),因此BERT模型(模型6)較BIC模型(模型7)對于低能傾向于產生n、H+，而對于高能更傾向于產生3H、α(但仍較BIC模型少)。對于60 MeV和200 MeV，模型9的H+產額較模型7增加了約49%和35%；而重離子產額增加了96%和69%，因此對比高斯能譜可見，EMY模型的影響隨入射質子能量的減小而增大，尤其是對于治療眼部腫瘤的60 MeV的影響較顯著。

表2 60 MeV和200 MeV次級粒子產額比較Table 2 Comparison of secondary yields for 60 MeV and 200 MeV

2.3 材料影響

表3列出了模型4、6、7 3個具有代表性模型在不同材料模體中的次級粒子產額。由于肺密度較小，射程較長(Bragg深度約60 cm)，為對比需要，表3中給出了肺模體厚30 cm和90 cm兩種情況，其他材料模體均厚30 cm。可見，除γ和e-外，其他與水的類似，特別是重子產額表現出參數化LHEP模型與所研究材料的依賴性不大。從模擬時間看，在模體厚度一定的情況下，材料的密度越大，模擬時間越短。部分原因是：密度越大，粒子射程越短(如肺-0.1 cm/30 cm(網格厚0.1 cm，總厚30 cm)粒子輸運達30 cm，但水、軟組織等僅約17 cm)，輸運時間也就降低了。此外，肺材料中易發生散射的e-產額特別高(是其他材料的116～260倍)。若增大肺模體的厚度和網格厚度(0.3 cm)，追蹤粒子模擬直至其“死亡”或離開模體，其模擬時間不及精細網格(厚度為0.1 cm) 30 cm厚的模體的長，這說明Geant4中粒子模擬在邊界跨越上花費的時間更長，較模體厚度對粒子輸運的影響大。QGSP_BIC模型和QGSP_BERT模型與材料的關聯性增加，其中骨中各重子產額相對較高。整體上看，在模體厚度一定的情況下，材料的密度越大，模擬時間越短。

2.4 網格厚度與Bragg峰-入口劑量比

第二，資金匱乏。我國的養老服務項目在目前主要依靠政府購買的方式來維持，由于政府購買服務的資金比重較為偏低，導致服務受眾群體范圍變窄，很多需要幫助的老人不能享受到相應的服務，從而使養老的福利性嚴重不足[9]。另外，由于智慧養老直接服務群體和間接服務群體的支付能力有限等因素的制約，再加上智能設備的購買和維護成本高等原因，限制了養老服務對象的積極參與；專業服務人才隊伍建設成本高；政府資金投入力度不夠，資金獲取渠道單一；系統編程制作與系統維護成本高，支撐系統軟件的硬件設施成本高等，都是當前智慧養老服務模式建設所面臨的挑戰。

為研究不同能量照射下，不同網格厚度對Bragg峰-入口劑量比的影響。計算了入口和Bragg峰前后1 cm處分別劃分為厚1 mm、0.1 mm或0.01 mm網格的深度劑量曲線和Bragg峰-入口劑量比。圖3示出了60 MeV和高斯能譜對模型4的比較(其他模型與之相似)。可見，對于高斯能譜，精細的網格(0.1 mm或0.01 mm)和1 mm厚網格對劑量曲線無影響。但對于60 MeV，0.1 mm與0.01 mm的較一致(差異<1%)，但1 mm的與它們的有偏移。>160 MeV的情況的與高斯能譜的相同，本文未給出。本工作也研究了其他能量1 mm和0.1 mm厚網格差異，發現直到160 MeV，兩者的差異才小于2%。因此計算Bragg峰-入口劑量比，對于低能質子(<160 MeV)，劑量模擬有必要劃分精細的網格(約0.1 mm)；而對于高能質子(≥160 MeV)，1 mm的網格可接受。由于一般蒙特卡羅程序算法處理邊界跨越均需耗費大量時間用以判斷，所以過于精細的網格劃分需考慮時間成本。

表3 高斯能譜模型4、6、7不同材料模體中的次級粒子產額Table 3 Secondary yields of model 4, 6 and 7 with Gaussian in different material phantoms

圖3 網格厚度對模型4 PDD的影響Fig.3 Influence of grid thickness on PDD of model 4

圖4為Bragg峰-入口劑量比隨入射質子能量和模體材料的變化。由于模型4～8均較相近，僅給出了模型5作代表。從圖可見，隨入射質子能量的增加，Bragg峰-入口劑量比單調降低，說明隨著質子射程的增加，更多的能量損失在傳播路徑上，使得Bragg峰處的治療優勢大為降低。其中，模體材料的變化對Bragg峰-入口劑量比的影響很小，遠小于模型9的EMY處理的影響(5%～13%)。圖4中也給出了Wroe等[6]模擬的結果，Wroe等使用了適用范圍0～150 MeV的PRECO模型。其中60 MeV時，Wroe等和模型5的結果較接近，但200 MeV時，Wroe等模擬的結果較本工作采用模型的高，這可能是Wroe等超限使用了PRECO模型。

圖4 入射質子能量和模體材料對Bragg峰-入口劑量比的影響Fig.4 Influence of H+ energy and phantom material on dose ratio of Bragg peak to entrance

圖5為Bragg峰深度隨入射質子能量和模體材料的變化。由于各模型的Bragg峰位置相同，圖中只代表性地給出了Bragg峰對應的深度。其中水、軟組織和肌肉的Bragg峰深度變化規律相近，骨的Bragg峰深度變化相對緩慢些。此Bragg峰深度可用1個二次多項式擬合(R2>0.99)，如對軟組織，Bragg峰深度=-17.475 88+0.536 97E+0.004 15E2；對骨，Bragg峰深度=-11.681 63+0.369 25E+0.002 98E2。其中，能量E的單位為MeV，深度的單位為mm。圖5中也給出了Wroe等[6]的Bragg峰深度，從圖可見，與本模擬的結果較接近。

2.5 絕對劑量

表4給出了Bragg峰處60 MeV和250 MeV單個質子沉積的絕對劑量。模型6和8在研究的能量范圍內均使用了PRECO+BERT模型，數據相同，在表中未進行區分。可見質子能量越小，Bragg峰處沉積的劑量越大。60 MeV的Bragg峰絕對劑量約為250 MeV的5倍，而60 MeV的Bragg峰-入口劑量比為250 MeV的2倍，因此60 MeV的入口劑量僅為250 MeV的2.5倍。靶區獲取的絕對劑量決定著腫瘤療效，在靶區劑量閾值一定的情況下，60 MeV質子束前端對正常組織的輻射損傷相對較好。

圖5 入射質子能量和模體材料對Bragg峰深度的影響Fig.5 Influence of H+ energy and phantom material on Bragg depth

表4 Bragg峰處單個入射質子沉積劑量Table 4 Deposit dose for single H+ at Bragg peak

但從射程角度看，60 MeV的Bragg峰深度在水及等效介質中僅3 cm，而250 MeV的高達37.5 cm，因此在治療深度腫瘤方面，高能量的質子束具有優勢。此外，通過旋轉分割照射，也可減少高能質子對相同部位正常組織的輻射損傷。

模型4～8的結果較接近，模型9的絕對劑量明顯偏低，且此差異隨質子能量的增加而增大(如60 MeV降低4%，250 MeV降低16%)。模型9使用QGSP_BIC模型外，添加了G4EmExtraPhysics( )、G4HadronElasticPhysics( )、G4QStoppingPhysics( )、G4IonBinaryCascade-Physics( )、G4NeutronTrackingCut( )處理，致使其模擬結果偏離其他物理模型及文獻發表數據。盡管Geant4hadrontherapy例物理清單中指出此為簡潔版物理模型，但INFN使用的復雜物理模型不易獲取。此例題與Geant4新版本反復發布的簡潔物理模型，易使模擬者錯誤采納，應引起模擬者注意。

3 結論

本工作選取9個Geant4常用于治療質子束能量范圍(60～250 MeV)的物理模型進行比較研究。結果發現：若僅使用電磁物理模型，雖整體劑量分布上可接受，但微觀重子產額有缺失；參數化驅動模型LHEP的模擬時間最短，但不能產生復雜重子；級聯模型BERT和BIC的生物損傷輻射防護的預警性更強，尤以BIC模型更甚，因此QGSP_BERT模型、QGSP_BIC模型和FTFP_BERT模型較適合放療質子束模擬。從平衡網格劃分與模擬效率的角度看，<160 MeV的質子需0.1 mm厚網格，≥160 MeV的質子1 mm厚網格即可。低能量(如60 MeV)質子束的Bragg峰-入口劑量較高，前端對正常組織的輻射損傷相對較好。然而，由于能量與射程的二次多項式增長關系，在治療深度腫瘤方面，高能量的質子束具有優勢。通過旋轉照射也可減少高能量的質子束對正常組織的輻射損傷。

Geant4中hadrontherapy例提供的物理模型QGSP_BIC_EMY，與Geant4新版本多次發布，易被錯誤采納。本工作驗證了hadrontherapy例物理模型與其他物理模型的差別，以提醒使用者。

[1] 唐勁天. 腫瘤質子放射治療學[M]. 北京：中國醫藥科技出版社，2004.

[2] 劉云鵬，湯曉斌，謝芹，等. 質子放療過程中射束與人眼體作用的蒙特卡羅模擬計算[J]. 原子能科學技術，2010，44(增刊)：602-607.

LIU Yunpeng, TANG Xiaobin, XIE Qin, et al. Monte-Carlo simulation of proton radiotherapy for human eye[J]. Atomic Energy Science and Technology, 2010, 44(Suppl.): 602-607(in Chinese).

[3] 許淑艷. 蒙特卡羅方法在實驗核物理中的應用[M]. 2 版. 北京：原子能出版社，2006.

[4] 鞠志萍，曹午飛，劉小偉. 質子散射角分布的蒙特卡羅模擬[J]. 物理學報，2009,58(1):174-177.

JU Zhiping, CAO Wufei, LIU Xiaowei. Study of scattering angular distribution of proton using Monte-Carlo method[J]. Acta Physics Sinica, 2009, 58(1): 174-177(in Chinese).

[5] 謝朝，鄒煉，侯氫，等. 質子束治療中非均勻組織的等效水厚度修正研究[J]. 物理學報，2013,62(6):068701-5.

XIE Chao, ZOU Lian, HOU Qing, et al. Study of inhomogeneous tissue equivalent water thickness correction method in proton therapy[J]. Acta Physics Sinica, 2013, 62(6): 068701-5(in Chinese).

[6] WROE A J, CORNELIUS I M, ROSENFELD A B, et al. The role of nonelastic reactions in absorbed dose distributions from therapeutic proton beams in different medium[J]. Med Phys, 2005, 32(1): 37-41.

[7] Geant4 is a toolkit for the simulation of the passage of particles through matter[R/OL]. http:∥geant4.cern.ch/.

[8] CIRRONE G A P, CUTTONE G, ROSA F D, et al. Validation of Geant4 physics models for the simulation of the proton Bragg peak[C]∥Conf Rec 2006 IEEE Nucl Sci Symp. [S. l.]: [s. n.], 2006: 788-792.

[9] TOSHITO T, BAGULYA A L, IVANCHENKO A V, et al. Validation of new Geant4 electromagnetic physics models for ion therapy applications[J]. Progress in Nuclear Science and Technology, 2011, 2: 918-922.

[10]PAGANETTI H. Nuclear interactions in proton therapy: Dose and relative biological effect distributions originating from primary and secondary particles[J]. Phys Med Biol, 2002, 47: 747-764.

Effect of Different Geant4 Physical Models on Simulation of Radiotherapy Proton Beam

LIN Hui1, XIE Cong1, ZHANG Yong-jun1,2, XIONG Zhen-yu1,WU Dong-sheng1,*, CAO Rui-fen3, FDS Team3

(1.SchoolofElectronicScience&ApplicationPhysics,HefeiUniversityofTechnology,Hefei230009,China; 2.DepartmentofBasicCourses,BengbuMedicalCollege,Bengbu233000,China;3.InstituteofNuclearEnergySafetyTechnology,ChineseAcademyofSciences,Hefei230031,China)

The Monte Carlo method is usually used to simulate the feature of the external radiotherapy of proton beam to optimize the clinical therapy scenario. The selection of a suitable physical model is very critical to ensure the correction and the efficiency of the simulation. This work studied nine Geant4 physical models, which are often used in the simulation of radiotherapy of proton beam with energy of 60-250 MeV. The dosage features and the microcosmic secondary yields in different materials were compared. The result shows that, although the standard EM, the low energy Penelope model and the low energy Livermore model can output the right dose distribution, they can’t simulate the baryon yields. The parametrization driven model LHEP and the Quark gluon string-pre-compound model both can’t simulate the complex ion transportation. Thus the above five models are not enough for the research calculation in the radiation protection and the radiation damage. The QGSP_BIC_EMY model in Geant4 hadron therapy example is very bad for underestimating the dose ratio of Bragg peak to entrance and absolute dose. The QGSP_BERT model, QGSP_BIC model and FTFP_BERT model are suitable for the radiotherapy proton beam with energy of 60-250 MeV.

physical model; proton radiotherapy; microscopic yield; dosage characteristic

2014-03-10;

2014-06-25

中國科學院戰略性先導科技專項資助項目(XDA03040000)；合肥工業大學科學發展研究基金資助項目(2013HGXJ0193，J2014HGXJ0093)；合肥工業大學自主創新項目資助(2012HGZY0007)

林 輝(1973—)，女，安徽阜陽人，副教授，博士，核技術醫學應用專業

*通信作者：吳東升，E-mail: wudongsheng@hfut.edu.cn

TL84

A

1000-6931(2015)07-1290-08

10.7538/yzk.2015.49.07.1290