-
本文使用的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碳离子基础数据,只有α和β的深度分布不同,其余数据均相同。
数据序号 αion,βion约束 是否考虑细胞
存活率标准差考虑贴壁
效率的方式基础数据1 αion>0,βion≥0 是 仅考虑0 Gy 基础数据2 αion>0,βion≥0 否 仅考虑0 Gy 基础数据3 αion>0,βion≥0 是 考虑多剂量点 -
为了更方便地展示不同数据处理方式对生物剂量分布计算的影响,本研究的碳离子治疗计划使用一个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×106,射束点个数为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值对应的曲线拟合优度R2均大于0.97。LET = 100 keV/μm时,细胞存活曲线的拟合结果如图2所示。可以发现,是否考虑细胞存活率(SF和sf)的标准差拟合得到细胞存活曲线略有差异,但总体上差别不大。
-
在仅考虑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较大时,拟合时是否考虑细胞存活率标准差对α值影响较大,对β值影响较小。 -
将细胞存活率标准差纳入到细胞存活曲线拟合中,使用式(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%。将2.2的结果与2.3的结果相对比可以发现,针对本工作中的细胞克隆存活数据进行存活曲线拟合时是否考虑细胞存活率标准差会对α,β值产生较大的影响,而以两种不同方式考虑贴壁效率对α,β值的拟合结果影响较小。
-
为了评估对治疗计划的影响,以覆盖靶区95%体积的最小生物剂量D95,靶区最大生物剂量Dmax和靶区最小生物剂量Dmin作为治疗计划生物剂量分布的评估指标,如表2所列。为了评估在相同的物理计划参数下,不同的生物参数对计划生物剂量分布计算值的影响,采用基础数据1进行计划优化(Planset1),分别使用基础数据2(Planset2)和基础数据3(Planset3)重新计算生物剂量分布,以Planset1的D95,Dmax和Dmin为基准,分别计算Planset2和Planset3相对Planset1的偏差,如表3所列。可以发现,在处方剂量为2~6 Gy(RBE)时,使用基础数据2重新计算所得的结果与Planset1相比,总是高估了Dmax,Dmin和D95。使用基础数据3重新计算的结果于Planset1相比,并不总是低估Dmax,Dmin和D95,具体情况与处方剂量的大小有关。注意到Planset1和Planset2的Dmax相对偏差的绝对值总是明显大于Planset1和Planset3的相对偏差的绝对值,即是否考虑细胞存活率标准差对治疗计划中靶区最大生物剂量Dmax的影响比以不同方式考虑贴壁效率造成的影响更大,这与实验数据中得到的α,β数据的分析所得的结论相吻合。对于Dmin,两组相对偏差的绝对值较为接近,两者大小关系没有明显随处方剂量变化而变化的趋势。对于D95,随着处方剂量的上升,Planset1和Planset2之间的相对偏差绝对值逐渐下降,而Planset1和Planset3之间的相对偏差绝对值逐渐上升,两者趋势相反。与Planset1相比,在处方剂量为2~6 Gy时,Planset2总是高估D95;而Planset3在处方剂量为3~6 Gy时,总是低估D95。本研究中各个处方剂量下靶区各生物剂量参数Dmax,Dmin和D95中的相对偏差的绝对值最大为9.45%,与3.1和3.2中α,β值的相对偏差绝对值比,偏差范围减小。
单位:Gy(RBE) 处方剂量 Planset1 Planset2 Planset3 Dmax Dmin D95 Dmax Dmin D95 Dmax Dmin D95 2 2.13 1.69 1.89 2.33 1.77 1.97 2.19 1.76 1.90 3 3.20 2.48 2.85 3.38 2.55 2.93 3.26 2.57 2.82 4 4.28 3.27 3.80 4.37 3.37 3.87 4.31 3.37 3.70 5 5.32 4.16 4.74 5.39 4.27 4.79 5.33 4.28 4.54 6 6.39 5.01 5.69 6.41 5.13 5.71 6.35 5.12 5.37 处方剂量 100%×(Planset1-Planset2)/Planset1 100%×(Planset1-Planset3)/Planset1 Dmax Dmin D95 Dmax Dmin D95 2Gy(RBE) −9.45% −4.21% −4.41% −2.92% −4.15% −0.79% 3Gy(RBE) −5.66% −3.09% −2.76% −1.78% −3.92% 0.86% 4Gy(RBE) −2.04% −3.03% −1.94% −0.70% −3.05% 2.62% 5Gy(RBE) −1.28% −2.53% −1.04% −0.12% −2.78% 4.22% 6Gy(RBE) −0.35% −2.22% −0.30% 0.51% −2.14% 5.68% 图5展示的是射野中心轴线上,以处方剂量归一化的深度剂量分布。可以发现,随着处方剂量增加,在射野入口区Planset1-3的相对生物剂量都降低了。由于计划是基于基础数据1优化得到,所以Planset1在SOBP区域不同处方的相对剂量差异可以忽略。Planset1和Planset2在SOBP峰区的剂量均匀性与Planset1相比均出现了恶化。在SOBP末端,均出现局部的尖锐部分,且Planset3的峰区末端比Planset2的更尖锐。考虑到SOBP末端的LET较高(>200 keV/μm),其LET与α,β的关系是由实验数据的外推得到,因此外推拟合处理对α,β值的影响会体现在SOBP末端。值得指出的是,随着处方剂量的增加,Planset2和Planset3的靶区相对生物剂量均匀性逐渐变好,这可能是因为处方剂量增大,α,β值的差异作用减弱所致。
Influence of Cell Survival Experiment Data Processing on the Calculated Values of Biological Dose Distribution in Carbon Ion Radiotherapy
doi: 10.11804/NuclPhysRev.40.2022127
- Received Date: 2022-12-12
- Rev Recd Date: 2023-02-15
- Available Online: 2024-02-04
- Publish Date: 2023-12-20
-
Key words:
- carbon ion beam /
- RBE /
- LQ model /
- treatment planning /
- biological dose
Abstract: The influence of different data processing methods on the fitting values of linear square(LQ) model parameters and the calculated values of biological dose distribution used in treatment planning were investigated in the cell cloning survival experiment. Four sets of LET-α and LET-β data were obtained based on the LQ model fitted to the data from the clonogenic survival experiments of carbon ion irradiated cells from T98G cells, according to whether the standard deviation of the cell survival rate was taken into account, and the plating efficiency was taken into account in two different ways, respectively, and three of them were selected to build the basedata of three datasets of the matRad treatment planning system for the treatment planning study. When only the 0 Gy plating efficiency was considered, whether the fitting considered the standard deviation of survival rate had little influence on the α value when the LET was low, but had a greater influence on the β value; With the LET getting larger, it had a greater impact on the α value and a lesser impact on the β value. When considering the standard deviation of survival rate, the deviations of fitting results in 2 different ways to consider the plating efficiency were small. Except for 200 keV/μm, the relative deviation of α values between the two groups did not exceed ±10%. Planset2 always overestimated Dmax, Dmin and D95 compared to Planset1. The results illustrate that whether the standard deviation of survival rate were considered or not had a greater influence on the fitted α, β coefficients than considering the plating efficiency in 2 different ways. The relative deviation range of biological dose distribution’s calculation values among different datasets in the treatment planning system was smaller than the relative deviation range of α and β values of each group obtained by fitting the experimental data.
Citation: | Xingzhu DUAN, Yuanyuan MA, Hui ZHANG, Zhongying DAI, Qiang LI, Xinguo LIU. Influence of Cell Survival Experiment Data Processing on the Calculated Values of Biological Dose Distribution in Carbon Ion Radiotherapy[J]. Nuclear Physics Review, 2023, 40(4): 619-627. doi: 10.11804/NuclPhysRev.40.2022127 |