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