目前,基于逆向优化方法的调强放疗(intensity-modulated radiation therapy,IMRT)已成为放射治疗的主要技术手段[1]。医学物理师通过调整多个参数(如射束数量、射束角度、靶区和危及器官的剂量学指标和权重等)进行优化,在给予肿瘤致死剂量的同时尽可能减少周围正常组织的剂量[2]。但是,在临床实践中逆向优化也存在着一些挑战,比如物理师需要在优化过程中不断地进行试错操作,以找到具有患者特异性的优化目标参数及相应权重值;治疗计划质量也会因物理师的个人经验和治疗计划的复杂程度出现较大的差异。为了解决上述问题,一些研究者提出了基于方法论的自动计划方法:基于经验的计划设计(knowledge-based planning,KBP)、基于协议的自动迭代优化(protocol-based automatic iterative optimization,PB-AIO)、多标准优化(multi-criteria optimization,MCO)等[3]。
本研究提出一种基于优化参数树搜索算法(Optimization Parameter Tree Search Algorithm,OPTSA)的自动计划方法,该方法模拟人工计划的设计流程,通过算法迭代调整优化目标参数来逐步改善IMRT计划质量,直到满足终止条件。优化参数的调整过程则通过Eclipse TPS 15.6(美国Varian Medical Systems, Palo Alto)自带的脚本应用程序接口(Eclipse Scripting Application Programming Interface,ESAPI)来实现。本研究选取了直肠癌病例对上述方法进行评估,并比较了临床上的人工计划和基于OPTSA的自动计划在剂量学上的差异。
资料与方法1.病例资料:选择北京大学肿瘤医院放疗科2015年7月至2019年4月期间已经完成放疗的直肠癌术前患者20例。年龄40~70岁,中位年龄为61岁。所有患者均无放疗禁忌证。对患者的治疗部位进行CT扫描,层厚设置为5 mm。
2.勾画信息:在直肠癌IMRT计划中,既往已勾画的感兴趣区(region of interest,ROI)主要包括如下两类:靶区(GTV为肿瘤靶区,CTV为临床靶区,PGTV和PTV是经GTV和CTV外扩得到的计划肿瘤靶区和计划靶区)和OAR[股骨头、膀胱和辅助结构[BF](Avoidance)],其中,Avoidance定义为PTV前侧和后侧凹型区域(其原因是考虑小肠在腹部及盆腔中位置的不确定性,定位CT中的位置与治疗时可能存在较大差异,准确勾画的小肠难以准确评估其剂量),而且因为计划系统中提供的工具——正常组织目标(normal tissue objective,NTO)是各向同性的,采用Avoidance这种半自动生成的结构来控制小肠等正常组织的剂量。本次试验的所有病例均采用相同的临床处方剂量标准,PGTV 50.6 Gy/25次,PTV 41.8 Gy/25次,要求靶区95%以上的体积接受的剂量应达到处方剂量。
3. IMRT自动计划创建:本研究以资深物理师经过试错优化后得到的人工计划为参考,评价基于OPTSA创建的IMRT自动计划的优劣。自动计划的整体设计过程是在Eclipse TPS 15.6上利用其自带的ESAPI脚本接口,借助基于C#语言的Visual Studio 2019的脚本编译软件完成。通过“VMS.TPS”命令空间中的“Application”接口打开指定计划号的患者并创建新的IMRT计划。所有计划采用统一的参数设置要求,包括:射野能量为10 MV,射野采用等角7野均匀分布,角度分别为0°、51°、102°、153°、204°、255°、306°,计划等中心点选取在PTV体积中心,准直器角度与靶区外轮廓相适应,即铅门Y与PTV的长轴保持平行,铅门开口由系统自动设定。
4.优化参数树搜索算法:计划优化作为计划设计中最为关键的一环,该方法通过OPTSA模拟人工计划在优化中的试错过程,以寻找最理想的优化目标参数集。目前OPTSA方法是基于贪心算法的概念来探索参数空间,其基本工作流程如图 1所示。OPTSA的探索空间为一个层次树结构,其中每一个树节点都对应一个特定的优化目标参数集合。
(1) 初始优化目标参数:在OPTSA自动计划优化的初始阶段,将一组预定义的优化目标参数加载到光子优化器中,执行第1次优化,如表 1所示。初始优化参数参考临床计划经验并选择相对宽松的OAR体积剂量约束条件以扩展待搜索的整体参数空间。
(2) 综合OAR得分:对于未扩展的树节点,需要一个标准来指导参数空间的探索。根据人工优化的经验,自动计划中靶区的目标通常会在一个单独的优化循环中率先被满足。为了控制树节点的探索顺序,直观地定义了一个综合OAR得分,即Evasum,具体定义如下:
$ Ev{a_{{\rm{sum}}}} = \sum\limits_{i = 0}^n {P_i^k} /\sum\limits_{i = 0}^n {P_i^0} $ | (1) |
式中,n为OAR的数量;Pik和Pi0分别表示在第k次迭代和第1次迭代时,第i个OAR的平均剂量及2%体积和50%体积接受的绝对剂量的加权和,具体描述如下:
$ P_i^k = \alpha D_{{\rm{mean }}(i)}^k + \beta \left( {D_{2\% (i)}^k + D_{50\% (i)}^k} \right) $ | (2) |
$ P_i^0 = \alpha D_{{\rm{mean }}(i)}^0 + \beta \left( {D_{2\% (i)}^0 + D_{50\% (i)}^0} \right) $ | (3) |
式中,Dmean(i)k和Dmean(i)0分别为第k次迭代和第1次迭代时第i个OAR的平均剂量;D2%(i)k、D50%(i)k和D2%(i)0、D50%(i)0分别为第k次迭代和第1次迭代时第i个OAR的2%和50%体积所接受的最小剂量。算法中根据器官的属性可将OAR分类为串行器官和并行器官,其中利用α和β对两类器官的平均剂量和剂量体积参考点剂量权重进行调节[4]。对于串行器官,α和β设为1和2;对于并行器官,α和β设为2和1。本研究中涉及病例的危及器官均考虑为并行,故α和β按并行器官设置。综合OAR得分作为不同剂量分布之间的排序标准,其中得分越低的节点被认为是更有前景的分支。
(3) 参数空间探索:图 1中所展示的迭代过程是由以下3个主要的操作组成:①选择。是指选择要探索的下一个节点,即从优先队列中选取Evasum值最低的节点。在初始化阶段,根节点的初始优化目标参数集会被放入优先队列并进行第一次迭代。②扩展。描述了基于父节点生成多个子节点的过程。在目前的实现中,相互连接的两个树节点中仅有一个优化目标参数被改变。改变优化目标参数后的子节点数量与直肠癌计划中勾画的OAR数量相等。目前,在直肠癌病例中共选择了3种OAR(股骨头、膀胱和Avoidance)。在每次迭代中,从3个OAR中选择一个,降低此方案下相应OAR平均剂量限制的10%,生成新的优化目标参数集。③仿真。是根据上一节点的优化结果继续对新生成的优化目标参数集进行常规优化。优化结束后,得到OAR的DVH信息并计算当前子节点的综合OAR得分(Evasum值)。根据每个节点对应的Evasum值对子节点进行由大到小排序并放入优先队列。由于树节点数目随树深度的变化呈指数级扩展,因此,OPTSA进行了一个修剪过程来降低问题的复杂性,不符合临床标准的节点被排除在探索过程之外,修剪标准定义如下:①靶区按照临床处方剂量来限制,包括:PTV,D95%≥41.8 Gy;PGTV,D95%≥50.6 Gy;CTV,D100% > 41.8 Gy;GTV,D100% > 50.6 Gy。任一靶区的剂量数值不满足上述标准,则执行丢弃子节点操作。②OAR剂量则期望尽可能低。因此,如果子节点的Evasum值低于父节点,则将其添加到优先队列中,否则选择丢弃子节点。
(4) 迭代和终止条件:重复选择、扩展、仿真3个步骤进行迭代优化,直到满足以下两个条件之一时停止迭代:优先队列为空,无可继续探索的节点;达到最大迭代次数后,跳出循环。为了确定最大迭代次数,选取20例计划进行60次迭代收敛测试。Evasum值在迭代循环中被记录,以量化每次迭代的改进程度。在测试案例中,进行60次迭代循环的自动计划平均运行时间为38 min,这也是本单位高年资物理师对直肠癌案例常规所需的计划时间。跳出迭代后,得到OPTSA探索到的最理想计划。
利用解析各向异性算法(analytical anisotropic algorithm,AAA)计算出自动计划的剂量分布和剂量体积直方图(dose volume histogram,DVH)[5]。按照临床靶区处方剂量标准判断自动计划是否符合临床要求。如果符合临床剂量,则获得一个合格的IMRT自动计划,否则,OPTSA会自动执行二次优化循环,直到满足处方剂量标准:对不满足处方剂量标准的靶区,提高其优化目标参数权重为原先的105%,而满足剂量标准的ROI优化目标参数权重保持不变。重复进行二次优化和剂量判断,直到所有靶区均满足临床处方剂量的要求。
5.临床评估:利用ESAPI接口获得自动计划和人工计划中相关的剂量学指标[6],并对两者的计划质量进行对比评估。具体的剂量学指标包括:最大剂量Dmax;平均剂量Dmean;靶区接受大于107%处方剂量的相对体积V107%;靶区95%体积所接受的绝对剂量D95%;OAR体积的2%所接受的绝对剂量D2%;OAR体积的50%所接受的绝对剂量D50%;OAR接受大于20 Gy剂量的相对体积V20;OAR接受大于30 Gy剂量的相对体积,V30;靶区适形性指数(conformity index,CI)[7]。CI=(TVRI/TV)×(TVRI/VRI),其中,TVRI是处方剂量覆盖的靶区体积,TV是靶区体积,VRI是处方剂量覆盖的体积;靶区均匀性指数(homogeneity index,HI)[8]。HI=|D2-D98-Dp。其中,D2和D98分别为靶区体积2%和98%所接受的绝对剂量;Dp为靶区处方剂量。
6.统计学处理:采用SPSS 22.0软件,计量资料符合正态分布,用x±s表示。对自动计划和人工计划的所有剂量学指标采用配对t检验进行统计分析。P < 0.05为差异有统计学意义。
结果1.最大迭代次数的确定:迭代次数与Evasum值和运行时间的关系如图 2所示。记录所有病例从第1次到第60次迭代过程中每一次迭代的Evasum值和运行时间。每5次迭代的Evasum值的中位数分别为0.819、0.749、0.713、0.695、0.692、0.686、0.684、0.667、0.654、0.654、0.652、0.653。每10次迭代的运行时间的中位数分别为5.9、12.2、18.0、22.9、28.8、37.9 min。第45次迭代后,Evasum值趋于收敛。根据平均收敛速度,最终允许的最大迭代次数设置为50次。在本研究中,自动计划平均的运行时间为(28.1±3.6)min,人工计划为(36.7±4.6)min。
2.剂量覆盖和靶区适形性:自动计划和人工计划关于靶区覆盖情况的剂量学指标定量比较如表 2所示。其中,自动计划的PTV和PGTV平均D95%分别为41.9和50.8 Gy,而人工计划的PTV和PGTV平均D95%分别为42.7和51.1 Gy。所有的自动计划和人工计划都满足临床处方剂量的要求。对于PTV的平均V95%指标,在自动计划和人工计划中有相近的统计结果,分别为(99.0±0.26)%和(99.9±0.09)%。与人工计划中PTV和PGTV的平均Dmax(55.1和54.9 Gy)相比,在自动计划中PTV和PGTV的平均Dmax下降到53.9和53.7 Gy。相比于人工计划,自动计划中PTV的平均CI值和HI值升高(t=14.11、10.73,P<0.05)。PGTV的平均V107%自动计划与人工计划比较,差异有统计学意义(t=-4.36,P<0.05)。对于PGTV的CI和HI指标差异无统计学意义(P > 0.05)。在机器总跳数方面,自动计划的平均MU值比人工计划降低20.2%。
3.危及器官的保护:自动计划和人工计划有关OAR保护的剂量学指标如表 3所示。与人工计划相比,除Avoidance的平均Dmean指标外,自动计划的其他剂量学指标均有不同幅度的下降。其中,股骨头的平均V20(t=-5.57,P<0.05)与膀胱的平均V30(t=-5.16,P<0.05)指标降低的最为显著。相比于人工计划中膀胱的平均Dmean,自动计划降低了(1.33±1.13) Gy。在自动计划中,OAR的剂量保护是以Avoidance平均剂量的轻微增加为代价。然而,自动计划对于Avoidance热点区域有着更好的控制,其平均D2%指标有所下降(t=-1.04,P<0.05)。
为了更好地展现自动计划对OAR剂量的控制,从20例直肠癌患者中选取2例具有代表性的病例展示自动计划与人工计划各个OAR的DVH曲线对比情况如图 3所示。在图 3中,自动计划对于股骨头和膀胱的剂量控制要优于人工计划,但对于Avoidance的DVH曲线,人工计划有着更好的结果;自动计划股骨头和Avoidance的DVH曲线相比于人工计划均有明显的改善趋势。对于膀胱的DVH曲线,自动计划与人工计划有着相近的剂量结果,自动计划在计划优化中并没有对膀胱有更好地剂量控制。
讨论
本研究在Eclipse计划系统平台上开发了一种直肠癌IMRT自动计划方法。该方法借助PB-AIO的部分思路,通过ESAPI接口模拟人工做计划的流程进行计划设计,利用迭代的方法不断评估和调整计划优化中的目标参数,以获得最终的剂量分布,从而实现IMRT计划设计的自动化[9-12]。Tol等[13]在Eclipse平台的基础上开发了一个交互界面,通过自动地移动屏幕上鼠标光标,检测屏幕上每个ROI的DVH线的位置,迭代地改变目标。与上述方法相比,本研究利用Eclipse自带的ESAPI脚本接口,在整个计划设计过程中更加方便的去读取、修改,回输所需的剂量学指标,通过迭代的方式生成IMRT自动计划,结果证明了该方法的有效性和实用性。与人工计划相比,自动计划在靶区剂量覆盖和OAR剂量控制方面均有所提升。如果临床物理师对自动计划的结果不满意,所获得的自动计划也可以作为进一步人工干预的起点。此外,OPTSA可以同时在后台执行多个IMRT自动计划,以缩减计划设计中的人力和人工耗时。
目前,除了PB-AIO方法外,基于KBP和MCO的自动计划设计也有报道。KBP方法通过先验知识建立患者数据库利用机器学习模型进行训练,可用于预测新患者的剂量分布和DVH曲线[14-16]。在Eclipse平台中,基于KBP方法的RapidPlan自动计划已经在不同病种下实现[17-18]。MCO方法是基于“帕累托最优解”的概念,试图在靶区覆盖和OAR保护之间寻求一个最佳的折中方案[19]。MCO自动化方法也有了商业化的模块。与Rapidplan自动计划相比,基于OPTSA算法的自动计划并不依赖于临床病例库和模型建立,其模拟临床物理师的计划设计流程并通过迭代的形式去逐步改善IMRT的计划质量。
在OPTSA初始阶段,初始优化目标参数集主要是参考资深物理师的临床经验给出。当对OAR剂量目标施加严格的约束时,可以实现较低的OAR剂量分布,但往往是以热点的出现和靶区CI、HI变差为代价。因此,为了扩展整体参数空间和平衡靶区与OAR的关系,本方法在参考临床经验的基础上给出了相对宽松的剂量目标限值,证明在满足临床处方剂量要求的基础上,算法可以尽可能地压低OAR的剂量。当然,由于探索空间会随着树深度的增加而呈指数增长,导致OPTSA只能在有限的时间和空间内寻找更优的计划,即获取局部最优解而非全局最优。当不同地区或不同医院做自动计划拓展时,可以参考本文的初始优化目标参数集或略作修改。在OPTSA迭代优化的过程中,优化结果评估函数(综合OAR得分)是根据OAR分类情况给出统一的OAR剂量权重并进行计算。由于直肠癌放疗计划病理结构并不复杂,且OAR数目较少,在本研究中使用统一的计划参数集合可以适用于所有病例。在接下来的工作中将尝试为自动计划添加剂量预测模块,通过剂量预测模型获取每个患者独立的初始优化目标参数,指导初始优化的同时缩短OPTSA迭代搜索的过程。对于复杂病种,会综合考虑各个结构之间的联系,并为优化结果评估函数作适应性设计。
在一些自动计划中也会存在表现不佳的方面,例如,自动计划对于股骨头和膀胱有更好的剂量控制往往是以Avoidance剂量略微增加为代价的。但是表现不佳的原因最主要是由于大多数术前的直肠癌患者存在较大的PTV体积,PTV与膀胱体积过大则会产生交叠现象[20]。在所有参与试验的直肠癌患者中,膀胱与PTV交叠体积占膀胱总体积的平均比重为10.32%。交叠现象的出现,使得靶区和OAR的优化目标发生冲突,为了满足靶区剂量要求,没有足够的距离来降低OAR的剂量,从而无法对OAR的剂量进行更好的控制。尽管如此,自动计划总体上显示出更好的OAR保护效果。在满足靶区剂量标准的同时,每个OAR的剂量可以控制在尽可能低的水平。当一些病例出现靶区和危及器官严重交叠现象时,基于OPTSA算法的自动计划在计划质量上无法比人工计划有着明显优越的表现,甚至需要进一步的人工干预。在未来的工作中,将更深一步探究靶区和危及器官交叠问题对于自动计划及剂量的影响。
基于脚本接口的OPTSA方法可以很容易地集成到任何TPS中进行应用。其应用灵活,并可以根据需要进行扩展,包括:修改各种处方剂量要求,拓展其他病种(如宫颈癌、鼻咽癌等),增减所需的OAR数量。该方法在处理临床治疗计划中涉及的各种情况时性能更稳定。因此,它可以作为一个有效的工具来减少不同物理师,甚至不同治疗系统之间计划质量一致性的差异。
借助ESAPI脚本接口,本研究在Eclipse平台下开发了一种高效的搜索算法,并实现了直肠癌的IMRT治疗计划设计的自动化。通过迭代优化,OPTSA可以逐步调整优化目标参数,以不断改善IMRT计划质量。结果表明,该方法能有效缩短计划设计中耗费的人力时间,提高计划质量,减少人为因素对计划质量的影响。
利益冲突 本项研究工作无任何形式的利益冲突。Varian医疗系统公司为本项工作提供了部分支持。作者独立取得实验数据,并对数据及其统计分析的准确性和可靠性负完全责任
作者贡献声明 王翰林负责数据采集和分析;刘嘉城参与代码编写和调试;姚凯宁和岳海振负责统计分析;张健提供人工计划;王若曦和张艺宝负责审阅论文;吴昊负责设计方法学,指导整个研究
[1] |
Ezzell GA, Galvin JM, Low D, et al. Guidance document on delivery, treatment planning, and clinical implementation of IMRT:Report of the IMRT subcommittee of the AAPM radiation therapy committee[J]. Med Phys, 2003, 30(8): 2089-2115. DOI:10.1118/1.1591194 |
[2] |
Nelms BE, Robinson G, Markham J, et al. Variation in external beam treatment plan quality:An inter-institutional study of planners and planning systems[J]. Pract Radiat Oncol, 2012, 2(4): 296-305. DOI:10.1016/j.prro.2011.11.012 |
[3] |
Hussein M, Heijmen B, Verellen D, et al. Automation in intensity modulated radiotherapy treatment planning-a review of recent innovations[J]. Br J Radiol, 2018, 91(1092): 20180270. DOI:10.1259/bjr.20180270 |
[4] |
Thieke C, Küfer KH, Monz M, et al. A new concept for interactive radiotherapy planning with multicriteria optimization:first clinical evaluation[J]. Radiother Oncol, 2007, 85(2): 292-298. DOI:10.1016/j.radonc.2007.06.020 |
[5] |
Van Esch A, Tillikainen L, Pyykkonen J, et al. Testing of the analytical anisotropic algorithm for photon dose calculation[J]. Med Phys, 2006, 33(11): 4130-4148. DOI:10.1118/1.2358333 |
[6] |
Zhang F, Zheng M. Dosimetric evaluation of conventional radiotherapy, 3-D conformal radiotherapy and direct machine parameter optimisation intensity-modulated radiotherapy for breast cancer after conservative surgery[J]. J Med Imaging Radiat Oncol, 2011, 55(6): 595-602. DOI:10.1111/j.1754-9485.2011.02313.x |
[7] |
van't Riet A, Mak AC, Moerland MA, et al. A conformation number to quantify the degree of conformality in brachytherapy and external beam irradiation:application to the prostate[J]. Int J Radiat Oncol Biol Phys, 1997, 37(3): 731-736. DOI:10.1016/s0360-3016(96)00601-3 |
[8] |
Lauve A, Morris M, Schmidt-Ullrich R, et al. Simultaneous integrated boost intensity-modulated radiotherapy for locally advanced head-and-neck squamous cell carcinomas:Ⅱ-clinical results[J]. Int J Radiat Oncol Biol Phys, 2004, 60(2): 374-387. DOI:10.1016/j.ijrobp.2004.03.010 |
[9] |
Song Y, Wang Q, Jiang X, et al. Fully automatic volumetric modulated arc therapy plan generation for rectal cancer[J]. Radiother Oncol, 2016, 119(3): 531-536. DOI:10.1016/j.radonc.2016.04.010 |
[10] |
Zhang X, Li X, Quan EM, et al. A methodology for automatic intensity-modulated radiation treatment planning for lung cancer[J]. Phys Med Biol, 2011, 56(13): 3873-3893. DOI:10.1088/0031-9155/56/13/009 |
[11] |
Xhaferllari I, Wong E, Bzdusek K, et al. Automated IMRT planning with regional optimization using planning scripts[J]. J Appl Clin Med Phys, 2013, 14(1): 4052. DOI:10.1120/jacmp.v14i1.4052 |
[12] |
Hazell I, Bzdusek K, Kumar P, et al. Automatic planning of head and neck treatment plans[J]. J Appl Clin Med Phys, 2016, 17(1): 272-282. DOI:10.1120/jacmp.v17i1.5901 |
[13] |
Tol JP, Dahele M, Peltola J, et al. Automatic interactive optimization for volumetric modulated arc therapy planning[J]. Radiat Oncol, 2015, 10: 75. DOI:10.1186/s13014-015-0388-6 |
[14] |
Schreibmann E, Fox T. Prior-knowledge treatment planning for volumetric arc therapy using feature-based database mining[J]. J Appl Clin Med Phys, 2014, 15(2): 4596. DOI:10.1120/jacmp.v15i2.4596 |
[15] |
Shiraishi S, Moore KL. Knowledge-based prediction of three-dimensional dose distributions for external beam radiotherapy[J]. Med Phys, 2016, 43(1): 378. DOI:10.1118/1.4938583 |
[16] |
Zhu X, Ge Y, Li T, et al. A planning quality evaluation tool for prostate adaptive IMRT based on machine learning[J]. Med Phys, 2011, 38(2): 719-726. DOI:10.1118/1.3539749 |
[17] |
Schubert C, Waletzko O, Weiss C, et al. Intercenter validation of a knowledge based model for automated planning of volumetric modulated arc therapy for prostate cancer. The experience of the German RapidPlan Consortium[J]. PLoS One, 2017, 12(5): e0178034. DOI:10.1371/journal.pone.0178034 |
[18] |
Hussein M, South CP, Barry MA, et al. Clinical validation and benchmarking of knowledge-based IMRT and VMAT treatment planning in pelvic anatomy[J]. Radiother Oncol, 2016, 120(3): 473-479. DOI:10.1016/j.radonc.2016.06.022 |
[19] |
Monz M, Küfer KH, Bortfeld TR, et al. Pareto navigation:algorithmic foundation of interactive multi-criteria IMRT planning[J]. Phys Med Biol, 2008, 53(4): 985-998. DOI:10.1088/0031-9155/53/4/011 |
[20] |
Wu B, Pang D, Simari P, et al. Using overlap volume histogram and IMRT plan data to guide and automate VMAT planning:a head-and-neck case study[J]. Med Phys, 2013, 40(2): 021714. DOI:10.1118/1.4788671 |