一直做放疗计划优化的同事问我你们搞建模的人总提灵敏度分析这东西到底能给临床带来什么我当时给他打了个比方你想调好一台复杂的设备如果不知道哪个旋钮对最终输出影响最大就只能一个接一个去试试一次动一次测量。肿瘤生长模型的伴随灵敏度分析干的事情就是把这套“逐个试”变成一次数学上的反向积分。做时空放射治疗优化时模型里可能有几十个与肿瘤增殖、氧合、辐射杀伤有关的参数每个参数还随时间、空间变化如果还靠老办法逐个扰动来算梯度优化器根本跑不动。这篇文章我就把整个思路捋一遍模型长什么样伴随方程从哪来Matlab里怎么落地以及算出来的灵敏度在时空放疗优化里怎么用。适合正在做计算生物医学、放疗物理建模或者在Matlab里写PDE约束优化的人参考。1. 灵敏度分析的起点肿瘤生长模型里埋着多少个不确定旋钮1.1 一个可用的肿瘤生长模型参数远比你想的多做时空放射治疗优化的时候模型不能太简单太简单了无法刻画辐射响应也不能太复杂复杂到参数本身根本定不准。常用的折中方案是一类反应扩散方程把肿瘤细胞密度 (c(t,x)) 的演化写成增殖、扩散、辐射杀伤三个过程的叠加[ \frac{\partial c}{\partial t} \nabla \cdot (D \nabla c) r c \left(1 - \frac{c}{K}\right) - \alpha c - \beta c \int_0^t e^{-\lambda(t-s)} d(s) ds ]这里面 (D) 是扩散系数(r) 是增殖率 (K) 是环境容纳量 (\alpha) 和 (\beta) 描述辐射直接杀伤和间接杀伤 (d(s)) 是照射剂量率。临床实际里不同患者的 (r) 能差好几倍 (D) 可能在低氧坏死区和活跃增殖区也不一样更不用说放疗过程中肿瘤还会再增殖、乏氧区域会重新氧合。这些参数一旦设错后续优化出的剂量分布再漂亮也是纸上谈兵。所以灵敏度分析要回答的问题很直接目标函数对哪个参数最敏感如果某个患者某个区域的 (r) 偏移了20%最终治疗方案的TCP肿瘤控制概率会掉几个百分点在梯度优化框架里这其实就是求目标函数对每个模型参数的偏导数。1.2 逐参数扰动为什么在时空放疗优化里行不通最早我接触灵敏度分析时第一反应是有限差分把某个参数 (\theta_i) 扰动一个小量 (\varepsilon)重新跑一遍正问题比较目标函数差异。公式很简单[ \frac{\partial J}{\partial \theta_i} \approx \frac{J(\theta_i \varepsilon) - J(\theta_i)}{\varepsilon} ]但当参数维度上升到几十个、上百个之后这个方案就崩了。每一轮优化迭代都要重新计算正问题也就是重新求解整个时空域上的PDE如果还要反复扰动才能得到梯度相当于每走一步都要跑几十遍正问题。我做三维网格、时间步长取小时级的时候一个正问题就要几分钟几十个参数就意味着
