《1 前言》

1 前言

中子俘获治疗 ( neutron capture therapy , NCT ) 是在细胞层次上将具有高传能线密度 ( LET ) 的带电粒子有选择性地定位辐射到肿瘤区的一种现代治疗技术｡早在 1932 年 Chadwick 发现中子和 Goldhaber 于 1934 年发现了具有超常热中子吸收截面( 3 840 Barn ) 的天然同位素硼 ( 10B ) 之后不久, 人们便提出了 NCT 的概念｡ 10B 吸收热中子后迅速变为 11B , 紧接着分解为具有总动能为 2.33 MeV 的 α 粒子 ( 4He ) 和 7Li 粒子, 这两种粒子飞行方向相反, 在一条直线上｡粒子在肌体中的作用范围约为 12~ 13 μm ( 同细胞大小可比拟 )｡该反应构成了 BNCT 的物理基础[1] ｡ Gordon Locher 于 1936 年首先提出了 BNCT 的原理｡假定硼可有选择性地富集于肿瘤区, 由于硼具有超常的热中子吸收截面,硼俘获热中子反应主要发生在肿瘤区, 反应后的能量在约 10 μm 内沉积, 因而, 被热中子照射后, 肿瘤区吸收的剂量就远比正常组织高 , 就可以达到杀死癌细胞, 保护正常细胞的目的｡因此, 从理论上讲, 如果能找到理想的含硼药物，注入人体后使硼在肿瘤区的浓度比正常组织高 3 倍以上（目前的 BPA 已能满足此要求），治疗效果就会很好。

BNCT 用于临床治疗的历史已有 20 年，主要优点是：ａ. 当前的硼化合物是无毒的，且投放肿瘤区和正常组织的硼浓度富集度差可以做到最大化，使得中子与硼作用产生的高能量、短射程的 α 离子在杀死肿瘤细胞的同时，对正常组织的伤害达到最小；ｂ. 中子束技术日趋先进，可以来自加速器，也可来自微型反应堆；ｃ. 中子束强度可控，既可发出超热中子束，也可发出热中子束，根据肿瘤深度，确定中子束类型；ｄ. BNCT 适合治疗不宜手术的晚期癌症患者，且治疗的肿瘤类型由脑胶质瘤，向其他包括肝癌、肺癌、皮肤癌等发展；ｅ. BNCT 的最大好处还在于，它给患者带来的治疗副作用少，没有放疗、化疗那样的不适反应；ｆ. BNCT 治疗后的病人存活期通常要长于其他治疗方法，且成本不高。因此 BNCT 目前已成为肿瘤治疗的一个较好选择。

由于被治疗的病人总会带有一定的放射性，这会给医护工作者带来一定的心理负担和健康危害。另外，在拥有众多人口的医院内建医用反应堆治疗室，在许多国家获得批准的可能性较小，这也多少妨碍了 BNCT 的推广应用。我国 BNCT 起步较晚，但发展速度很快，在周永茂院士和王忠诚院士的推动下，这项工作已取得实质进展。经过不懈努力， 30 kW 的微型医用反应堆于 2010年在房山中国原子能科学研究院建成，目前正在进行临床前的动物实验。

世界上开展 BNCT 研究的国家有 20 多个，美、日、欧、阿根廷等国家和地区代表世界最高水平。然而，用于临床治疗的剂量计算软件目前只有６个，它们是 NCTPlan､MacNCTPlan[2] ､BNCT - rtpe[3] ､ SERA[3] ､JCDS[4]和笔者研究团队研制的 MCDB[5] ｡ MCDB 尚需临床应用的检验，目前只做了一些比对计算和效率测试。

《2 MCDB 核心算法介绍》

2 MCDB 核心算法介绍

MCDB 的前、后处理是完全自主开发的，中间剂量计算涉及中子－光子耦合输运问题，这部分是在引进蒙特卡罗程序 MCNP4C[6] 基础上开发研制的。由于 MCNP 程序是个通用程序，直接用于剂量计算，一来原程序没有剂量计算所需的响应截面库 ( Kerma 因子），二来计算效率较低，无法满足临床对计算精度和计算时间的要求，三来 BNCT 的网格数巨大，已超过 MCNP 的最大几何块 106 限制。临床计算要求每个体素（voxel）网格剂量的统计误差 5%，计算时间不超过 2 h，这对 MCNP 串行计算来说是做不到的。因此，需要对程序做多方改进。

当前 MCDB 相对 MCNP 的主要改进包括：a. 增加了前处理，即体素（voxel）模型的构造，把 CT 和 MRI 数据自动转化为 MCNP 的输入文件；b. 突破了 MCNP 程序对几何块、几何面及计数的限制，根据计算机内存，可以动态调整这些参数；c. 针对特殊体素网格模型，发展了一套快速粒子径迹算法；d. 计数矩阵和材料矩阵构造，克服了 MCNP 程序初始化时间随网格数增加呈指数增长的不足；e. MPI 并行计算。

下面给出剂量计算的步骤，总剂量 D 由硼､、中子、质子四部分剂量组合而成，可表示为：

式中：Wc 为硼剂量的复合生物效应值； 分别为光子、快中子、质子（氮俘获）剂量的相对生物效应值。

输运计算涉及中子－光子耦合输运，剂量通过通量和 Kerma 因子响应得到：

这里 Kerma ( j, E )为第 j 个体素网格、能量为 E 的通量剂量转换因子，数据来自 ICRU-632000[7] ； 为第 j 个体素网格、能量为 E 的通量；V ( j ) 为第 j 个体素网格的体积。

由于 BNCT 复杂的几何和对方向、能量的精度要求很高，因此，通常用 MC 方法进行模拟，MC 计算通量的公式如下：

式中：S0 为源强；N 为模拟的样本数（粒子数）；wm 为第 m 次碰撞的粒子权重（初值 w0 =1）；dm 为粒子在第 j 个网格走过的径迹长度；分别为粒子的宏观吸收和总截面[7] 。

对 dm 的计算，涉及粒子飞行线与被穿过几何块交点的计算，由于MC 输运计算时间的 70％花在处理碰撞和粒子径迹的计算上，缩短计算时间的关键就是要缩短计算径迹的时间。针对体素模型单一结构网格的特点，发展快速计算径迹 dm 的算法是完全可能的。

《2.1 体素模型的构造》

2.1 体素模型的构造

如图 1 所示，在 Goorley 的文章中，介绍了 Snyder 简化头模型[8]，该模型包括头皮、颅骨和大脑三部分。体素模型将基于头的大小，确定 x､y､z 方向的范围，形成一个能正好罩住头部的盒子，其定义域为  ×  ×  。在每个方向分别确定网格数目，假定 x､y､z 方向的网格数分别为 M､N､L，则总网格数为 M ×N ×L，每个网格大小相同，为一长方体。根据 CT 灰度值确定每个网格的材料成分。困难在两种物质交界面的处理，它含有两种以上的物质，混合后形成新的物质，解析模型只有 4 种基本物质 ( 见图 1 (a) )，转换为体素模型后，新增了 282 种物质，总物质数达到 286 种 ( 见图 1 (b)）。

这里重点介绍 MCDB 的两种算法。

《图1》

图1 Snyder基准模型

Fig.1 Snyder benchmark model

《2.2 中心点方法》

2.2 中心点方法

为了节约内存，降低模拟问题的复杂度，用每个网格中心点的材料作为整个网格的材料。这种处理只有在质量守恒前提下有效。经验证，当体素网格边长小于 5 mm 时，就是在最薄的皮肤层质量误差也小于 1 %，满足守恒要求。由于中心点方法不增加材料数目，因而可以节约大量内存，串行对计算时间影响不大[9]，但对并行计算来说，通讯量的减少意味着并行加速比的提高。

《2.3 快速粒子径迹算法》

2.3 快速粒子径迹算法

基于整数的运算速度最快，每个网格被映射到单位立方体网格，在单位立方体网格下，每个网格节点均可用整数 ( i , j , k )表示，网格定义域可表示为：，它的相邻 6个网格分别为，。建立一个与网格对应的三维材料矩阵和计数矩阵。对所有网格进行计数。由于采用半开半闭的区间表示，每个粒子的属性是唯一的，因此，不会出现 MCNP 那样的由于二异性导致的丢失粒子现象。由于粒子轨迹是在立方体网格中穿越的，粒子沿某个飞行方向穿过的几何块径迹不必按ＭCNP 的模式计算，利用等差关系就可快速求出粒子沿飞行线穿过的多个网格的径迹。同时，到碰撞点的距离也不必在每个网格抽取，已抽取的距离可以一直用下去，直到碰撞发生。其后再重新抽取到达碰撞点的距离。详细算法可参考文献[6,10]。

《3 有效性测试》

3 有效性测试

《3.1 剂量计算》

3.1 剂量计算

MCDB 用 BNCT 基准模型的测试已经做过[6,10] ,这里选择一个实例, 笔者等从天坛医院得到一个病人的 CT 图 ( 共 43 片, 见图 2 ), 首先用 MCDB 的前处理进行三维重构 ( 见图 3 (a) ), 然后生成体素模型 ( 见图 3 (b) ), 体素模型网格大小为两种: a. CT8模型,每个网格由 8×8 像素组成, 共分 64 ( x 方向 ) ×64 ( y 方向 ) × 43 ( z 方向 ) = 176 128 个网格, 每个网格大小为: 0. 370 3 mm × 0. 370 3 mm × 0.3 mm ; b. CT4模型, 每个网格由 4×4 像素组成, 共分 128 ( x 方向 )×128 ( y 方向 ) × 43 ( z 方向 ) = 704 512 个网格, 每个网格大小为: 0. 185 2 mm × 0.3 mm × 0.3 mm｡

《图2》

图2 CT图

Fig.2 CT pictures

《图3》

图3 MCDB 前处理示意图

Fig.3 Sketch of MCDB pre-processor

入射中子束使用宽谱超热中子束，其中 1 %为快中子 ( E 10 keV )，10 %为热中子( E 0.5 eV )，其余 89 %为超热中子 (0.5 eV E 10 keV )。在 3个能量区间中，能谱均为 1/E 分布，入射中子均匀分布在半径5 cm的圆面上，源强为 1010 cm-2 · s-1 ，模拟 1000 万粒子。为了确保热中子、快中子和剂量的同步收敛，采用源能量偏倚发射。图 4 给出 MCDB 和 MCNP5 热中子剂量、快中子剂量和剂量比较( MCNP5 的几何块､几何面仍受到 106 的限制, 但 mesh 计数没有限制 )｡图 5 给出两者的数值差距的比较, 由于数据量巨大, 只好通过二维平面来表示,  用 - 表示 MCNP5 剂量结果; -表示 MCDB 剂量结果 ( 剂量单位: Gy ), 表示网格号｡若 x =y, 则 ( x , y )正好在 45°线上｡最后表 1 给出两种体素模型下, MCNP5 和 MCDB 串行计算时间的比较｡图 6 分别给出中子､光子剂量面图｡

《图4》

图4 MCDB 和 MCNP5 热中子剂量、快中子剂量和剂量比较

Fig.4 The dose comparison of MCDB and MCNP5 about thermal neutron，fast neutron and photon

《图5》

图5 MCDB 和 MCNP5 剂量计算结果偏差比较

Fig.5 Difference comparison between MCDB and MCNP5 doses

《表1》

表1 MCDB与MCNP5计算时间比较

Table 1 Comparison of MCDB and MCNP5 in computational time

注: 计算机为 Pentium IV 3.0 GHz PC 机

《图6》

图6 照射部位的确定

Fig.6 Determination of irradiation location

《3.2 结果分析》

3.2 结果分析

从图 4 给出的结果，通过肉眼很难看出 MCDB 和 MCNP 两者的差距，但从图 5 可以看出两个程序的计算结果差异还是存在的，特别是 计数，由于其主要来自中子碰撞产生的次级光子，数目较少，因此统计误差较大，两个程序的计算结果都存在一定的不确定性。由于 MCNP 程序的权威性，它的计算结果通常被用来标定其他程序结果的正确性，因此可以初步肯定 MCDB程序的计算结果是正确的。 MCDB还发展了组合网格模型及其相关网格下的快速粒子径迹算法，这些算法的有效正确性也得到了验证[6] 。

《4 结语》

4 结语

通过一个实际模型全过程的模拟，验证了 MCDB治疗计划软件的正确有效性。下一步 MCDB 将用于临床实际，根据遇到的问题，再做进一步的改进。

致谢软件发展到今天，得到周永茂院士的精心指导，他提出了许多有益的建议。借此，向他表示深深地感谢！