FLUKA综述：粒子治疗多功能技术仿真实验室
首先，蒙特卡罗程序的广泛应用得益于它可以描述射线与物质相互作用的诸多细节，比如对加速器、束流输运系统以及扫描系统的研究；但应用到粒子治疗领域时，不仅物理模型很重要，生物剂量的准确计算更为重要，尤其是比质子重的离子，比如4He、12C离子治疗，因此，就要求FLUKA计算剂量时能够把相对生物有效性(RBE)和生物学模型考虑进去。其次，准确预测次级射线对于射程程验证的研究有极大的帮助，粒子治疗中产生的正电子和瞬发伽马都是射程在线验证的优秀选手。对于用于粒子治疗的能量段，FLUKA改进的核反应模型能够更精细地描述治疗时的混合射线场，对于低能的部分也做了精细的处理。另外，为了适应医院的应用情景，最近更新的Flair版本提升了DICOM处理的功能，能够更好的读取和显示DICOM CT、RTDOSE，轻松转化数据格式，也可以将大量RTPlan的束流参数一键转换成FLUKA可以读取的输入文件[1][2]。对于连续多重物理过程的硼中子俘获治疗(BNCT)，FLUKA也可以对其进行建模和分析。
本文将从以下几个方面进行讲述：
核反应模型以及生物剂量；
射程验证研究 ；
DICOM数据处理；
加速器和束流线；
硼中子俘获治疗；
其他领域的应用研究；
讨论总结。
质子与物质的相互作用类型总结如下表[3]：
可以看出，粒子治疗发生的反应主要是与核外电子和原子核的库伦散射。跟核外电子的库伦散射是能量损失和剂量沉积的主要形式，决定了初始带电粒子在患者体内的射程；与核子的库伦散射则主要影响束流的横向运动，决定了束斑的散射情况，会增加侧向伴影；非弹性核反比如(p, n)(p, γ)(p, x)应则产生丰富的二次射线，包括中子、正电子、伽马射线等，这些次级射线为间接的在线射程验证研究带来了多种可能性，是目前质子重离子研究的一个热点领域。
FLUKA中，精细描述重带电粒子的方程是：
其中Tmax是传递给电子的最大能量：
不同于通常的近似计算，FLUKA将计算上述公式的所有的项。其中C项称壳修正项，粒子能量低的时候尤为重要，为使结果准确，FLUKA采用了ICRU49号报告中的值[4]。上述公式中的所有项FLUKA都做了详细考究，采用了权威报告中的数据和结果，并及时更新至最新结果。这保证了计算精度，并与大量实验结果吻合的很好，下图是FLUKA计算的深度剂量曲线与实验的验证结果。实线和虚线是计算曲线，点是实验值。
质子和重离子在反应模型上还是有很大区别的，主要体现在重离子能够产生更多碎片，使得在布拉格峰跌落后沿呈现尾巴剂量，并且由于碎片有横向速度，使得横向剂量也需要多重考虑。
Glauber-Gribov级联和高能碰撞；
(广义)核内级联；
预平衡发射；
蒸发反应/削裂反应/裂变反应和最终的退激发。
对于能量为亚GeV/n附近的离子，FLUKA处理这类离子的模型是RQMD-2.4，RQMD是相对论量子分子动力学模型(relativistic quantum molecular dynamics model)，它可以嵌入到核内级联中运行。但是RQMD只能处理快核反应过程，其激发出来的碎片能量较低，需要回到PEANUT模型中处理。下图是FLUKA+RQMD的计算不同能量下135 MeV/n的碳离子打靶反应截面和实验结果的对比，曲线是FLUKA计算值，圆点是实验结果，可见符合的相当好。
FLUKA中用玻耳兹曼主方程(BME)模型来处理低能部分，当离子能量低于150 MeV/n时用BME来处理，高于150 MeV/n时用RQMD来处理。通过BME方程数值积分抽样，FLUKA也用它来处理两核子的热核聚变反应。最新版本的FLUKA也把BME的事件产生器与PEANUT的预平衡发射模块接入，用于处理所有核子的第一能级的退激发。这种最新处理对于核反应的退激发过程描述极为重要，如下图所示100 MeV/的He核轰击石墨靶反应不同能量中子的产额计算结果，与实验结果吻合的很好。
这些复杂的核反应模型与商业TPS中的近似公式处理相比将更为准确，这也就是为什么蒙特卡罗被认为是剂量计算的金标准。
目前版本的FLUKA计算生物剂量时采用的生物学模型是LEM模型[6]，程序计算每一个Voxel体素的α和β参数，计算考虑了不同射线i对参数的总和贡献：
这样FLUKA就可以计算生物剂量，下图是C离子的生物剂量计算结果。
将来FLUKA也可以计算LET剂量，将更加有利于临床的研究，让我们拭目以待吧。
虽然每个医院的放疗都有严格的质量保证QA流程以达到精确治疗的目的。但实际上人们确实无法知道每一束质子、每一个射野是否都准确按照计划输送到肿瘤区域，因为现有技术还无法让人们“看到”质子轰击癌细胞的过程，虽然我们有很多先进的影像系统，如CT、PET、MRI等都无法做到。
所以目前射程验证的研究是粒子放疗研究中的一个热点，很多科研小组在这方面做了大量探索工作[7]。其中正电子发射成像和瞬发伽马成像是比较有前景的两个方向。蒙特卡罗也是这方面研究的有力工具，因为正电子和瞬发伽马都是核反应的产物，这就要求蒙特卡罗程序能够非常准确的计算次级射线。
正电子和瞬发伽马射线都是核反应的最后阶段发射的次级射线，因此它们对核反应的细节过程非常敏感。放射性核素的产额以及光子的发射产额不仅与核子能级有关，还跟自旋和宇称有关。尤其在布拉格峰区，入射粒子的能量已经低于核反应阈能，因此，出于对强子放疗的考虑，FLUKA低能核反应模型都做了特别的处理。
下图是质子的12C(p, x)11C和16O(p, x)15O反应截面计算结果与数据库EXFOR的实验值对比。生成同位素核11C 和15O 都是正电子发射体；x是反应的发射粒子1H、2H、3H和He。可见计算值与实验值符合的很好。
FLUKA&Flair工具包中提供了各种PET矩阵，用户可以根据自身需要对矩阵参数进行修改。
下图是用FLUKA计算的不同同位素正电子发射的分布。可见11C和15O是主要的贡献核素。该图的情况是用质子照射20 s后0~2 min的平均分布。
FLUKA也可以抽取一个截面，查看不同部位正电子发射情况。
用FLUKA研究瞬发伽马一般先建立与实际情况对应的几何模型，比如常用的Slit相机。
然后分析探测器的伽马能谱，如下图所示：
最后分析深度方向的伽马产额，将随着深度剂量曲线出现相应的分布变化，尤其是布拉格峰下降沿，伽马光子的产额迅速跌落，这也就是瞬发伽马用来监测质子重离子的射程的物理基础。
近些年，射程验证这个领域的研究涌现出很多新思路，比如2020年Changran报道了用核物理中的多普勒漂移的方法来计算射程。一些核物理的基础方法正向核技术临床应用方向发展，例如多普勒漂移法、飞行时间方法等。
简化和增强功能的DICOM接口；
新的DICOM编辑器可以执行简单的更改，例如匿名化DICOM文件；
体素生成可以覆盖执行布尔操作(如在几何图形中)的ROI上的材质；
增强的RTViewer与RTPLAN进行交叉检查计算；
RTPLAN处理为离子和PHOTON计划创建必要的光束；
将USRBIN转换为RTDOSE DICOM格式；
自动导出到RTDOSE和USRBIN的DVH；
下图显示了CT导入FLUKA的体素模型，同时也显示了RTDOSE。这一节将简单介绍它的处理过程。
Flair的DICOM模块主要基于Pydicom编码，pydicom是一个读取DICOM文件的开源程序语言，也是python的一个子模块。另外也调用了python的数值模块NumPy来处理DICOM文件。下图是Pydicom读取DICOM文件的一个简单范例。当然Flair的图形界面已经做好，用户一般情况下无需深入的代码内部。但是，由于是非商业软件，有些数据格式会出现报错的情况，这就需要读者进入到程序本身来查看原因，这需要有一定的编程经验，至少能读代码。
CT值也就是通常所说的Hounsfield单位，定义如下：
把CT图像转换成体素矩阵Voxel；
每个Voxel对应一个HU值；
根据Schneider的参数方法，将HU分成24个材料区间；
给每一个HU范围定义材料属性，包括组织成分和密度。
治疗计划系统生成的RTPLAN包含了大量的束流相关信息：束流类型(光子、质子、重离子)、射野方向、扫描点分布、束流强度(MU)、能量等等。手工生成FLUKA需要的source非常困难，新版Flair通过一键提取RTPLAN的重要参数生成BEAMSPOT文件，大大简化了束流模型编辑的工作量。但是FLUKA只认粒子数，不识别MU，所以每个MU对应多少粒子数需要用户根据治疗系统本身的情况而定，然后转换成粒子数；当然有的TPS直接给出的是粒子数，就不用转换。这一点需要特别注意。
下图左上是计划系统给出的剂量分布，右上图是FLUKA根据导出的数据验算的结果。右下图DVH。所以FLUKA&Flair可以用于病人治疗计划的质量保证QA。详情请见质子中国往期报道《FLUKA高级应用——病人QA中的MC剂量验证计算》。目前患者QA的程序通常是通过水箱实验完成，这种方式需要占用束流时间和人力资源，因此，独立的剂量验证计算是病人QA的高效替代方案。
当研究新的粒子治疗系统时，硬件的开发极为昂贵，蒙卡软件则可以用来初步研究各种参数匹配的结果，以初步判定使用哪种方案，这将极大地为企业和研究所节约时间和经费。2018年Cuccagna等人就提出了一个非常亮眼的质子直线加速器方案，目标是设计一种可以替代现有光子治疗的单室机型，可以把质子加速到232 MeV[8]。
该系统分上下两部分，下部为固定段，可以将质子加速到70 MeV，上部旋转段为高能段，可以提高到232 MeV。在他们的设计任务中采用了三种工具：RT Track、MADX-PTC、FLUKA&Flair。前两个用于加速管和输运线的设计，FLUKA和FLAIR则用于剂量配送系统(即扫描治疗头)的设计研究。
将设计好的扫描头结构，材料等参数输入到FLUKA中，比较重要的SAD，以及扫描场大小等。
然后从MADX的输出结果中获取束流的相空间参数。相空间参数就作为FLUKA的束流的源source。两个扫描磁铁还需要添加合适的磁场。关于扫描头方面的信息，请见质子中国往期报道《FLUKA高级应用：扫描治疗头的蒙特卡罗模型研究》。
FLUKA可以输出深度剂量曲线IDD，ISO处的束斑大小等临床需要的结果，这样就可以进一步判断系统设计是否合理，是否满足临床需求。
使用蒙特卡罗的好处是，研究人员和工程师可以方便的调节设计参数，使最终结果趋近甚至达到设计目标值，而不需要硬件的迭代。FLUKA也可以替代MADX完成前两步的工作，但目前来说还没有MADX用的方便。主要还是建模和计算量方面还有待改进。
最后，整个粒子治疗系统加速器物理和主要部件设计完成后，FLUKA还可以用于粒子治疗系统的辐射防护设计，尤其是辐射剂量的计算，将比经验公式更是更加可靠，结果更加直观。详情请见质子中国往期报道《FLUKA实操——质子治疗瞬发辐射剂量计算》。
对于有磁场的情况，FLUKA也可以通过添加特定的卡片和用户程序激活磁场。详情请见质子中国往期报道《FLUKA高级应用——质子在磁场中的运动》。
FLUKA也可用于硼中子俘获治疗(BNCT)的研究。Pazirandeh等人用FLUKA创建了基于一个电子加速器的中子源模型，并用于BNCT的研究。
为了得到最高的中子流强，他们对各种把材料及几何参数进行了优化，同时周围的屏蔽体材料及厚度也进行了优化计算。模型如下图，绿色的phantom是为得到不同处的剂量设置的。得到中子通量最高的是钨靶，可以达到1011(n/cm2s)。
在这项研究中，采用的是瓦里安2300 C/D电子加速器，电子能量可达20 MeV，23 mA。高能电子轰击到钨靶产生韧致辐射，韧致辐射的高能光子在通过轰击铅和铍靶，发生(γ, n)反应产生高能中子，热中子的能量最高可到5 MeV，计算的中子通量分布如下，对于不同的治疗，可选择不同的剂量率对应的位置。
经过调整之后，在Phantom的部位中子的能谱如下图。在大约0.05 eV处峰值为1.02×108±5%(n/cm2s)，且高能的部分比较干净。这种能谱的中子就特别适合黑色素瘤的硼中子俘获治疗。
以上综述了FLUKA在粒子放疗方面的应用研究，由于其丰富和精确的粒子物理及核反应模型，FLUKA也广泛用于宇宙射线、中微子、核探测器、辐照损伤等方面的研究。下图是2015年Kim等人用FLUKA计算GEM探测器的中子探测效率：
2020年Mazziotta等人用FLUKA研究了宇宙射线与太阳大气的相互作用，产生可以到达地球的次级粒子。在这项工作中，他们使用FLUKA综合计算大量次级粒子(如伽马射线、电子、正电子、中子和中微子)的产额。还通过将它们与撞击在太阳表面的宇宙射线的强度折合，来估算这些次级粒子在太阳下的强度以及地球上的通量。
本文简要综述了FLUKA在粒子治疗各个子领域的使用情况。Flair丰富的工具包更加扩展了其应用空间，从生物剂量计算，到射程验证研究，到DICOM文件处理，再到BNCT研究，FLUKA都可以完成；宏观上加速器的设计、辐射防护设计，FLUKA也是有力的工具。FLUKA&Flair就像提供了一个虚拟的实验室，里面有丰富的工具，能够模拟核物理范围下混合物理场的几乎所有环节。
但它也不是万能的，FLUKA的模型大都只能是静态模型，对于动态演化方面的研究几乎很难做到。例如Voxel就无法模拟呼吸运动，无法创建一个准临床环境，对于复杂的核反应过程例如快中子过程的模拟，FLUKA目前也是无能为力的。(质子中国 编辑报道)