与光子放疗相比，碳离子束用于肿瘤放射治疗具有物理学和放射生物学两方面优势：物理学上，它具有Bragg峰和较小的横向散射等特性，可以获得高的靶区剂量适形度；碳离子束Bragg峰区具有高的传能线密度(LET)，因而生物学上，它具有较高的相对生物效应(RBE)和较低的氧增比(OER)。碳离子放射治疗(CIRT)先后在日本和德国的几个治疗中心得到成功应用，并积累了丰富的科研成果和临床经验^[1–6]，随后中国也开启了碳离子放疗的设施建设、研究和临床工作^[7–9]。

截至到2019年底，全球超过22万患者接受了粒子治疗，其中重离子治疗的患者达2.8万人，占12.7%^[10]。据统计，中国大陆2020年共有147.53万人接受了放射治疗，而同年仅有5072名患者接受了碳离子放射治疗，可见与接受光子放疗的患者数目相比，接受碳离子放射治疗的患者数目极其稀少。光子放射治疗已具有百余年的发展历史，积累了大量临床经验，并掌握了针对不同肿瘤的光子处方剂量和正常组织器官的剂量限值。相比光子放射治疗，碳离子放射治疗的相关经验还很缺乏，这使得在应用碳离子束进行肿瘤治疗时缺乏可靠的肿瘤处方剂量和正常组织器官的剂量限值可供参考^[11]。为解决这个问题，相对生物学效应(RBE)的概念被引入到离子放射治疗中。RBE被定义为产生相同生物学效应所需参考射线剂量和高LET射线剂量的比值(参考射线一般为220 kV的X射线，ICRU 40, 1986)^[12]。通常，RBE值的大小受许多因素影响，例如细胞辐射敏感性、LET值、离子种类、剂量和生物学终点等。由于依赖关系复杂，因此对大量的数据进行汇总研究得出其变化规律是必要的^[13]。在碳离子治疗计划系统中，对RBE的计算是必须的，有的治疗计划系统直接使用细胞存活实验数据拟合线性平方(LQ)模型中的α和β参数计算RBE值^[14]，有的使用RBE理论计算模型^[15-16]。需要指出的是，在计划系统中使用RBE理论模型依然需要计算各剂量网格处对应的α，β参数，从而得到各剂量网格处的RBE值^[17]。此外，调试RBE计算模型通常也需要与实验获得的α，β参数进行比较^[18-19]，因此获得准确的α，β参数非常重要。

2013年Friedrich等^{[13, 20]}对大量文献中报道的重离子细胞存活实验据进行了汇总，建立了PIDE数据库，其中涵盖了不同实验室来源的实验数据及不同处理方式得到的α和β数据，从中可以发现对于相同的细胞，不同文献中的α和β值差异很大。此外，Friedrich等^{[13, 20]}对文献作者提供的细胞克隆存活原始实验数据进行了再次处理，其得到的α和β值与原文献公布的值也存在差异，如Furusawa等^[21]提供的人类唾液腺(HSG)细胞数据，图1所示。日本NIRS研究所使用了该HSG实验数据进行了脊形过滤器的设计^[22-23]，治疗计划设计^[14]和RBE微剂量动力学模型(MKM)参数的校准^[15]。从图1中可发现，Kanai等^[22-23]将实验获得的α和β数据进行了某种函数拟合和外推，获得了LET = 1~1 000 keV/μm范围内的α和β曲线。同时图1也说明，对于同一套细胞存活原始实验数据，不同作者给出的α和β值也不同，这可能是因为对数据处理方法或所用的数据处理软件造成的^{[13, 24]}。对细胞存活数据LQ模型拟合处理常见的方式有：细胞存活率标准差是否被考虑以及对细胞贴壁效率的不同考虑方式等。此外，将实验获得的α和β值用于治疗计划系统时，为了保证治疗计划生物剂量分布计算结果的均匀性符合临床要求，故需要采用平滑拟合等方法对实验获得的α和β值进行二次处理，正如Kanai处理Furusawa等^[21]的HGS细胞实验数据的方法^[22-23]。这些数据处理过程中的差异无疑会最终导致对治疗计划生物剂量分布估计的差异。为了评估这些差异，本文利用脑胶质瘤T98G细胞克隆存活实验原始数据研究不同的数据处理方式对LQ模型α，β系数及碳离子治疗计划的影响。

图  1  相同细胞克隆存活实验原始数据不同数据处理方法下的LET-α(a)和LET-β(b)关系(在线彩图)

本文使用的T98G细胞克隆存活实验原始数据，分别来自于利用甘肃武威重离子医院国产重离子加速器提供的260 MeV/u碳离子均匀扫描束进行的辐照实验和利用甘肃省人民医院瓦里安加速器提供的6 MV X射线进行的辐照实验。原始数据包含常氧状态下T98G细胞在各培养皿的存活细胞克隆数、各培养皿细胞种植数、照射剂量和LET值。在碳离子辐照细胞存活实验中，使用水等效系数为1.16的聚甲基丙烯酸甲酯(PMMA)对260 MeV/u的碳离子均匀扫描束进行降能，通过改变PMMA的厚度，从而获得LET约为15, 20, 25, 30, 35, 40, 50, 65, 80, 100, 125, 158, 200 keV/μm的碳离子束，并使用以上各LET的碳离子束对T98G细胞进行照射，对各LET的碳离子束分别照射了5个剂量点以获得细胞存活曲线。在X射线实验中，照射了6个剂量点。为了得到每个剂量点细胞存活率的标准偏差，实验设计了3个平行样。

通常细胞存活实验的数据处理方法如下，首先，通过0 Gy剂量点的实验数据得到细胞贴壁效率(PE)，即，$ \mathrm{P}\mathrm{E} = {\overline{n}}_{0}/{N}_{0p} $，其中，$ {\overline{n}}_{0} $，$ {N}_{0p} $分别为0 Gy剂量点平行样的平均细胞存活克隆数和细胞种植数。然后，通过各剂量点平行样中的平均克隆数、种植数和PE值计算得各剂量点的细胞存活率$ {\mathrm{S}\mathrm{F}}_{i} $，即$ {\mathrm{S}\mathrm{F}}_{i} = {\overline{n}}_{i}/\left({N}_{ip}\boldsymbol\cdot \mathrm{P}\mathrm{E}\right) $，其中$ {\mathrm{S}\mathrm{F}}_{i} $为归一化细胞存活率，$ {\overline{n}}_{i} $，$ {N}_{ip} $分别为i Gy剂量点平行样的平均细胞存活克隆数和细胞种植数。最后将各剂量点的细胞存活率和相应的剂量代入式(1)中进行拟合，即可得到LQ模型α和β参数：

由于细胞在辐射后可能会导致贴壁效率的改变，因而仅使用0 Gy剂量点的数据计算PE值不能反映其他剂量点的细胞贴壁情况。为了总体考虑所有剂量点的细胞贴壁效率，需将贴壁效率代入到拟合公式中，即：

其中：$ \mathrm{s}\mathrm{f} $为未归一化细胞存活率，即$ {\mathrm{s}\mathrm{f}}_{i} = {\overline{n}}_{i}/{N}_{ip} $。曲线拟合后可以获得细胞贴壁效率$ p $与LQ模型α和β参数。

本研究使用Python语言scipy工具包中optimzie.curve_fit函数来拟合细胞存活曲线，其中，sigma是optimzie.curve_fit函数的一个可选输入参数。当sigma输入为“none”时，表明拟合时不考虑细胞存活率的标准差，从而得到两组不考虑细胞存活率标准差的存活曲线和对应的α，β系数值；当sigma输入为SF或sf的标准差时，表明拟合时会考虑细胞存活率的标准差(SF为经过0 Gy归一化后的细胞存活率，sf为未经过归一化的细胞存活率)，从而得到两组考虑细胞存活率标准差的存活曲线和对应的α，β系数值。在此基础上，最终获得4组α和β值随LET变化的关系曲线。

为了评估不同处理方式获得的α，β系数对治疗计划的影响，本研究选用matRad计划系统进行碳离子生物剂量计算。matRad是一个基于Matlab语言的开源非商业治疗计划系统，主要用于治疗计划相关算法的开发与测试，支持光子、质子和碳离子束等多种射线的调强计划设计^[17]。在matRad中，碳离子的基础数据包含每个能量的深度剂量分布，深度LET分布及相应的α和β深度分布。每个能量的α和β深度分布由局部效应模型(LEM)计算获得。在碳离子剂量优化计算前，matRad除了计算每个笔形束的剂量影响矩阵，同时还计算每个笔形束的α影响矩阵和β影响矩阵。在剂量优化计算时，根据混合束模型计算每一个剂量网格中的剂量平均α和剂量平均β，从而得到对应的生物剂量。这一过程与日本NIRS的治疗计划^[14]及武威重离子医院使用的治疗计划系统类似，因此，可以将实验获得α、β与LET的关系引入到matRad治疗计划系统中。

由于实验数据中LET最大值为200 keV/μm，而matRad碳离子基础数据中LET范围大于200 keV/μm，因此需要通过数据外推的方式得到超出200 keV/μm的α和β值。此外，如果依据实验获得的α和β值随LET变化关系直接进行数据插值，往往会导致治疗计划系统的基础数据中出现局部跳跃性的LET-α和LET-β的分布，从而导致治疗计划系统计算出的生物剂量分布不满足临床要求。因此，还需要对实验获得的α，β值进行曲线拟合，从而得到光滑连续的LET-α和LET-β曲线，以保证治疗计划系统计算出的生物剂量分布符合临床要求。这一数据处理过程与Kanai处理Furusawa等^[21]的实验数据类似。首先，对实验中获得的4组LET-α和LET-β数据中的3组，利用Python脚本进行分段拟合，得到LET-α和LET-β值的关系曲线，范围在0~1 000 keV/μm。然后，根据matRad中碳离子深度与LET的分布关系通过线性插值得到3组matRad碳离子基础数据(基础数据1~3)，如表1所列。这3组matRad碳离子基础数据，只有α和β的深度分布不同，其余数据均相同。

为了更方便地展示不同数据处理方式对生物剂量分布计算的影响，本研究的碳离子治疗计划使用一个30 cm × 30 cm × 30 cm均质立方体水模体，其中只有一个直径为5 cm的水均质球作为靶体，其球心位于15, 10, 15 cm处，是该模体的唯一的感兴趣区域(VOI)。

治疗计划系统matRad的基础数据包括121个碳离子束对应水模体中的深度剂量分布曲线、横向束斑尺寸深度变化曲线以及LET随深度变化的曲线，能量层间隔设置为2 mm，剂量网格大小为3 mm × 3 mm × 3 mm。计划采用机架角和治疗床角为0°的单照射野，采用不同能量的碳离子笔形束进行照射，其能量范围为184.03~236.22 MeV/u。计划的体素点总数为8×10⁶，射束点个数为2 461。使用基础数据1作为matRad基础数据，设定靶区的目标函数为matRad中的square deviation形式，按照1 Gy(RBE)为步长，依次设定靶区分次处方剂量为2~6 Gy(RBE)，得到一组治疗计划，记为Planset1。将基础数据2和3分别设定为matRad基础数据，使用Planset1的点扫描权重重新计算生物剂量分布，分别记为Planset2和Planset3，并进行对比。考虑到三组计划使用的是同一组扫描点权重，三组基础数据具有相同的物理剂量分布和LET分布，因此，三组计划的生物剂量分布差异及差异的大小取决于α和β的分布差异及处方剂量的变化。

根据1.1节中的式(1)和式(2)，分别对各LET的T98G细胞克隆存活实验的数据进行了细胞存活曲线拟合，得到了不同数据处理方式对应的LQ模型的α，β参数与LET的关系。各数据处理方式下，所有LET值对应的曲线拟合优度R²均大于0.97。LET = 100 keV/μm时，细胞存活曲线的拟合结果如图2所示。可以发现，是否考虑细胞存活率(SF和sf)的标准差拟合得到细胞存活曲线略有差异，但总体上差别不大。

图  2  基于LQ模型的不同数据处理方式拟合T98G细胞克隆存活实验数据的细胞存活曲线结果(在线彩图)

在仅考虑0 Gy贴壁效率的条件下，处理T98G细胞存活实验的原始数据时，是否将细胞存活率标准差纳入到LQ模型存活曲线拟合，可得到两组α，β值、对应的拟合标准差以及两组α，β值的相对误差，如图3所示。两种处理方法得到的α值都在100 keV/μm附近出现峰，β值随LET增加而变大。在LET<50 keV/μm时，是否考虑细胞存活率标准差拟合细胞存活曲线时得到的两组α，β值之间差异较小；而在LET>50 keV/μm时，差距明显。定义两组α，β值的相对偏差$ {\Delta }_{\alpha }/\alpha $和$ {\Delta }_{\beta }/\beta $分别为${\Delta }_{\alpha }/\alpha = \left({\alpha }_{\mathrm{w},\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\mathrm{s}} - {\alpha }_{\mathrm{w}\mathrm{o},\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}s}\right) / {\alpha }_{\mathrm{w},\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{s}}$和${\Delta }_{\beta }/\beta = \left({\beta }_{\mathrm{w},\mathrm{ }\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\mathrm{s}} - {\beta }_{\mathrm{w}\mathrm{o},\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\mathrm{s}}\right) / {\beta }_{\mathrm{w},\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{s}}$。在图3(c)当LET<100 keV/μm时，两组α值相对偏差在±15%范围内；当LET值增大，两组α值相对偏差剧烈变化。除200 keV/μm外，两组α值相对偏差的绝对平均值小于10%。由于β值相对较小，特别是拟合时采用了对β系数的非零约束，导致有些数据点的β值为0，且β值随LET增大呈现增大趋势。因此，如图3(f)所示，在低LET区β值的相对偏差非常大，而在高LET区β值的相对偏差反而减小。综上，对T98G细胞，在只考虑0 Gy贴壁效率的前提下，LET较低时，拟合时是否考虑细胞存活率标准差对α值影响较小，对β值影响较大；而LET较大时，拟合时是否考虑细胞存活率标准差对α值影响较大，对β值影响较小。

图  3  使用0 Gy贴壁效率是否将细胞存活率标准差纳入曲线拟合的结果(在线彩图)

将细胞存活率标准差纳入到细胞存活曲线拟合中，使用式(1)和式(2)进行曲线拟合，分别得到仅考虑0 Gy贴壁效率和总体考虑贴壁效率这两种情况下的α，β值、对应的拟合标准差以及两组α，β值的相对误差，如图4所示。两种贴壁效率考虑方式所得的LET-α和LET-β曲线基本重合。将两种拟合方式所得α，β值的相对偏差定义为${\Delta }_{\alpha }/\alpha = \left({\alpha }_{eq1}-{\alpha }_{eq2}\right)/{\alpha }_{eq1}$和${\Delta }_{\beta }/\beta = \left({\beta }_{eq1}- {\beta }_{eq2}\right)/ {\beta }_{eq1}$。除200 keV/μm外，两组α值的相对偏差不超过±10%，相对偏差的绝对平均值小于4%。相比于$ {\Delta }_{\alpha }/\alpha $，两组β值的相对偏差$ {\Delta }_{\beta }/\beta $的范围则大得多。由于曲线拟合时对β值进行了非负约束，因此有部分β值为0，导致β值相差很小但相对偏差无穷大。随着LET增大，β值逐渐增大，两组β值的相对偏差绝对值降低到约4%。

图  4  以两种不同方式考虑贴壁效率的拟合结果(在线彩图)

将2.2的结果与2.3的结果相对比可以发现，针对本工作中的细胞克隆存活数据进行存活曲线拟合时是否考虑细胞存活率标准差会对α，β值产生较大的影响，而以两种不同方式考虑贴壁效率对α，β值的拟合结果影响较小。

为了评估对治疗计划的影响，以覆盖靶区95%体积的最小生物剂量D₉₅，靶区最大生物剂量D_max和靶区最小生物剂量D_min作为治疗计划生物剂量分布的评估指标，如表2所列。为了评估在相同的物理计划参数下，不同的生物参数对计划生物剂量分布计算值的影响，采用基础数据1进行计划优化(Planset1)，分别使用基础数据2(Planset2)和基础数据3(Planset3)重新计算生物剂量分布，以Planset1的D₉₅，D_max和D_min为基准，分别计算Planset2和Planset3相对Planset1的偏差，如表3所列。可以发现，在处方剂量为2~6 Gy(RBE)时，使用基础数据2重新计算所得的结果与Planset1相比，总是高估了D_max，D_min和D₉₅。使用基础数据3重新计算的结果于Planset1相比，并不总是低估D_max，D_min和D₉₅，具体情况与处方剂量的大小有关。注意到Planset1和Planset2的D_max相对偏差的绝对值总是明显大于Planset1和Planset3的相对偏差的绝对值，即是否考虑细胞存活率标准差对治疗计划中靶区最大生物剂量D_max的影响比以不同方式考虑贴壁效率造成的影响更大，这与实验数据中得到的α，β数据的分析所得的结论相吻合。对于D_min，两组相对偏差的绝对值较为接近，两者大小关系没有明显随处方剂量变化而变化的趋势。对于D₉₅，随着处方剂量的上升，Planset1和Planset2之间的相对偏差绝对值逐渐下降，而Planset1和Planset3之间的相对偏差绝对值逐渐上升，两者趋势相反。与Planset1相比，在处方剂量为2~6 Gy时，Planset2总是高估D₉₅；而Planset3在处方剂量为3~6 Gy时，总是低估D₉₅。本研究中各个处方剂量下靶区各生物剂量参数D_max，D_min和D₉₅中的相对偏差的绝对值最大为9.45%，与3.1和3.2中α，β值的相对偏差绝对值比，偏差范围减小。

图5展示的是射野中心轴线上，以处方剂量归一化的深度剂量分布。可以发现，随着处方剂量增加，在射野入口区Planset1-3的相对生物剂量都降低了。由于计划是基于基础数据1优化得到，所以Planset1在SOBP区域不同处方的相对剂量差异可以忽略。Planset1和Planset2在SOBP峰区的剂量均匀性与Planset1相比均出现了恶化。在SOBP末端，均出现局部的尖锐部分，且Planset3的峰区末端比Planset2的更尖锐。考虑到SOBP末端的LET较高(>200 keV/μm)，其LET与α，β的关系是由实验数据的外推得到，因此外推拟合处理对α，β值的影响会体现在SOBP末端。值得指出的是，随着处方剂量的增加，Planset2和Planset3的靶区相对生物剂量均匀性逐渐变好，这可能是因为处方剂量增大，α，β值的差异作用减弱所致。

图  5  深度-相对生物剂量分布($ {D}_{\mathrm{p}} $：处方剂量)(在线彩图)

在碳离子放射治疗中，决定RBE值的α，β参数的不确定性以及其带来的影响一直是被关注的焦点之一。本工作以T98G细胞克隆存活实验的原始数据为对象，探究了在数据处理过程中是否考虑细胞存活率标准差拟合和两种不同方式考虑贴壁效率拟合对LQ模型的α，β参数的影响，并进一步探究利用这些参数进行治疗计划设计导致的对生物剂量分布计算值的影响。对细胞存活实验获得的数据使用不同的数据处理方式会影响α，β值，进行存活曲线拟合时，考虑细胞存活率标准差以及考虑各剂量点的贴壁效率可以获得更准确的结果。虽然由实验获得的α，β值存在较大偏差，但是将其导入到治疗计划系统，最终获得治疗计划生物剂量分布计算值之间的差异将减小，且与处方剂量相关。此结果一方面是因为RBE计算公式的误差传递作用。另一方面可能是在计划系统中不同的能量相互混合使得α，β值正向和负向偏差相互抵消。此外，α和β值的偏差在低LET区域较小，在高LET区域较大，在碳离子治疗计划中高LET剂量被低LET剂量混合使得整个计划的剂量平均LET值并不是很高，因而最终导致生物剂量偏差减小。针对这一现象，我们对图1(a)中的LET-α曲线的高LET部分重新拟合，进行计划生物剂量分布计算值的差异比较，得到了类似结果。

在细胞存活曲线数据拟合中，两种不同方式考虑细胞贴壁效率似乎对本工作中的实验数据影响较小，这可能与细胞存活实验的操作流程有关。本工作使用的T98G细胞辐照实验数据，其操作流程是细胞在未经过辐照时就将培养瓶中的细胞转移到各剂量点的平行样的培养皿中，细胞在被辐照之前就完成了细胞贴壁。因此，理论上各剂量点的细胞贴壁效率应该一致，所以两种不同方式考虑贴壁效率对结果影响较小。然而，有些细胞存活实验的操作步骤是，先对细胞进行照射，然后再将辐照后的细胞转移到平行培养皿中，例如Kase等^[25]的实验。这种情况下使用未辐照细胞的贴壁效率去归一化不同照射剂量的细胞存活率，会引入误差。这一结果也提示，在使用其他实验室公布的α，β值时需要格外注意其数据处理方式和实验流程。

在碳离子放射治疗中，即使使用相同的细胞实验数据，仅是拟合细胞存活曲线时所采用的数据处理方式不同也会导致计划生物剂量出现明显的差异，这表明实际的碳离子放射治疗中生物剂量分布存在非常大的不确定性。由于碳离子放射治疗，不同的治疗中心使用不同的RBE计算模型，因此人们希望不同RBE模型之间剂量学的对比和处方转换实现临床经验的交换^[26-28]。由于RBE模型本身的不确定性会使得不同RBE模型之间的剂量转换工作变得困难，而且其意义受到影响。针对不同的RBE模型，如何制定更加安全可靠的计划将是碳离子放疗计划的方向之一。此外，由于RBE存在不确定性，在质子放射治疗领域已经提出了使用剂量平均LET取代治疗计划中生物剂量优化的方案^[29-30]。在碳离子治疗计划中，不依赖于RBE的计划方案还有待研究^[31]。

最后需要指出的是，由于本研究仅使用一种细胞数据实验数据说明数据处理对治疗计划生物剂量分布计算值及α，β系数的影响，取得的定量结果不能推广到其他细胞种类或模型数据。此外，研究只使用了简单的模体进行治疗计划的比较，在复杂的病例中各项指标如何变化还有待进一步研究。

在碳离子放疗中，由细胞克隆存活实验数据拟合得到的LQ模型α，β参数会因不同的实验数据处理方法而产生差异，进而对治疗计划系统(TPS)的生物剂量分布计算结果造成影响。因此，在使用其他实验室公布的α，β值进行研究时，需要注意其细胞存活曲线拟合方法和数据处理方式。