胶质母细胞瘤治疗反应预测：基于反应-扩散模型的个性化临床情景生成与验证

　　本研究提出一种基于反应-扩散方程（Fisher-Kolmogorov-Petrovsky-Piscunov方程）的数学建模方法，通过系统参数采样设计生成胶质母细胞瘤（GBM）治疗的个性化情景预测。模型利用临床MRI影像数据（T1增强和T2/FLAIR序列）进行初始化，在86%的影像间隔期内成功模拟出与实测肿瘤体积误差小于20%的结果，为复发GBM患者的治疗决策和临床咨询提供了创新工具。

　　

1 引言

胶质母细胞瘤（GBM）作为一种侵袭性原发性脑肿瘤，即使经过标准治疗仍常见复发，使得患者临床咨询面临巨大挑战。自20世纪初以来，数学建模已被用于理解肿瘤生长动力学和优化治疗策略。从最简单的指数增长模型到Gompertz模型和逻辑模型，这些模型通过不同方式模拟肿瘤细胞增殖的限制因素。

随着肿瘤学向个体化治疗方向发展，数学建模在优化和个性化治疗方面的潜力逐渐被探索，包括免疫治疗、化疗和放疗方案的优化。特别是在胶质瘤领域，Hormuth等研究者通过测量肿瘤细胞异质性预测未来生长，并利用影像和数学建模预测治疗反应。

新诊断的GBM患者通常接受最大安全切除、放疗和替莫唑胺化疗，但复发肿瘤尚无标准治疗方案。尽管过去20年GBM治疗未有突破性进展，但患者生存期的改善主要得益于支持性护理的进步。面对复发GBM患者的不良预后，患者、家属和临床医生之间的共同决策显得尤为重要。

本研究采用类似于财务规划师的思维方式，不试图模拟全球经济或预测股票市场，而是使用简单模型（如货币时间价值）和统计合理的先验回报率采样方法，基于患者当前状况生成一系列合理情景。同样，本研究的目标是确定一个简单的数学模型是否能够为复发GBM患者模拟未来2-3个月内个性化的、现实的可能结果范围，对应典型的随访影像间隔期。

GBM生物学复杂性极高，许多细节尚不清楚。影响肿瘤进展和治疗反应的因素包括基因组异常、新生血管形成、缺氧、免疫抑制、GBM干细胞和脑细胞外基质等。任何相关的速率方程最多只能是推测性的，并引入了无法在活体患者中测量的系数和初始条件。因此，我们专注于一个参数空间可抽样处理的模型，该模型可以从常规临床监测收集的影像数据中初始化并进行比较。

2 材料与方法

2.1 患者人群与数据预处理

磁共振成像（MRI）扫描来自Barrow神经学研究所（BNI）患者档案，经圣约瑟夫医院和医疗中心机构审查委员会（IRB）批准。数据包含来自75名独特患者的357次轴向监测扫描，患者年龄范围25-78岁（中位数62岁），这些患者先前均接受过最大安全手术切除，随后进行放疗和在大多数情况下使用替莫唑胺化疗。

影像处理流程包括以下步骤：首先，患者影像数据由MRI扫描仪生成的DICOM文件组成，识别信息被移除后，将轴向切片组装成单个三维图像并转换为NIfTI格式；其次，使用SPM-12将所有图像配准到一个共同时间点，并生成仅包含大脑的NIfTI文件（去除颅骨和眼睛），同时将大脑域进一步分割为脑脊液（CSF）、白质和灰质区域；最后，所有MRI扫描由神经影像专家手动分割，使用3D-Slicer图像处理平台"绘制"相关体素，并根据操作者的神经外科判断将其与解剖正常组织分离。肿瘤体素被分为三类：坏死核心（T1低信号）、增强肿瘤（T1对比增强）和肿瘤相关水肿（T2-FLAIR高信号）。

排除标准包括：SPM-12无法分割大脑域（可能由于扫描中的运动伪影）、T1对比后影像未见增强肿瘤、患者接受再次手术切除。最终获得46名独特患者的132个时间间隔的完整预处理扫描。

2.2 数学模型

建模工作聚焦于肿瘤总体扩张（或收缩）。肿瘤体积是总生存期的预后因素，可从作为数据源的监测影像中通过3D-Slicer平台轻松测量。

Fisher-Kolmogorov-Petrovsky-Piscunov（FKPP）方程是一个反应-扩散方程，最初提出用于模拟入侵物种的扩散，后来被建议用于描述GBM肿瘤细胞在整个脑实质中的扩散传播。该模型假设建模群体中的所有细胞具有相同的生长潜力，但GBM肿瘤的特征是存在缺氧或坏死"核心"。

作者先前提出了FKPP方程的修改版本，即所谓的"ASU Barrow"（ASUB）模型，将肿瘤细胞分为增殖（p）和静止（q）两个亚群：

?p/?t = ?·(D(x)?p) + g(p,q,t) - h(p,q,t)

?q/?t = h(p,q,t)

增殖细胞和静止细胞的密度分别由空间和时间依赖函数p = p(x,t)和q = q(x,t)给出。整体动力学很简单：增殖细胞以?·(D(x)?p)给出的速率扩散；它们以净人均速率g(p,q,t)生长，并以人均速率h(p,q,t)变为静止（即死亡或变为缺氧）。函数D(x)定义细胞浸润大脑的速率。静止细胞只是积累，反映高细胞性但没有净肿瘤生长。

在方程中，h(p,q,t)反映了治疗的净细胞杀伤效应；g(p,q,t)捕捉治疗抵抗和净增殖趋势。在本研究模拟的版本中，我们取：

g(p,q,t) = ρ(x)p[1-δ(p+q)]

h(p,q,t) = k(x)pδ(p+q)

其中ρ(x)和k(x)是空间依赖的最大生长和静止速率，δ是总细胞密度的单调递增函数。

2.3 模型模拟

SPM-12大脑分割定义了患者特定的计算域。在与颅骨和CSF的界面处施加无通量边界条件。模型可以使用Shampine等人的IRKC求解器高效积分。在典型的256×256×25患者大脑域上进行60天模拟大约需要18秒（现代个人计算机单CPU核心）。多个患者情景可以在多核机器上并行模拟。

生长和静止函数ρ(x)和k(x)在方程中是分段常数：如果x对应对比增强区域中的体素，则ρ(x)=ρ_c和k(x)=k_c；在水肿区域中ρ(x)=ρ_e和k(x)=k_e。扩散速率D(x)类似处理：当x对应白质束或肿瘤区域中的位置时，D(x)=D_w；在灰质中D(x)=D_w/2；在CSF中D(x)=0。

初始条件从各自的肿瘤分割中估算。Chang等人发现MR信号强度与GBM细胞密度之间存在近似线性关系，细胞性随T1加权对比后影像中的强度增加而增加，随T2/FLAIR信号强度而降低。在我们的模拟中，我们在被分割为增强肿瘤的区域中赋予一个初始种群，该种群随T1C图像强度线性增加 within区间I_c = [p^c₁, p^c₂]，并在被分割为水肿的区域中随T2/FLAIR信号在线性区间I_e = [p^e₁, p^e₂]内减少。

2.4 实验设计

我们不知道模型参数和初始条件的"真实"值，这些值可能在患者之间差异很大。我们的目标是确定产生临床代表性结果的值的可能范围。为了在如此大的参数空间中使计算易于处理，我们采用田口实验设计，这是一种拉丁超立方抽样。

对于模拟目的，表2中的每个参数被分配三个"水平"，对应"低"、"中"和"高"值。有3⁷ = 2187种可能组合，但我们仅通过优先考虑最敏感的参数k?₀和ρ?₀并根据正交数组抽样选择18种。正交数组由Kacker等人的表3中的第2-8列组成，包含18行。模拟使用两个最敏感参数的9种可能值对进行，对于每对，其余参数通过均匀抽样选择两个不同选项。目标是经济地估计模拟变异性作为参数的函数。

3 结果

我们总共评估了46名独特复发GBM患者的132个时间间隔的肿瘤进展。对于每个间隔，我们将初始扫描中的增强体素数量V₀与下一次扫描中的数量V₁进行比较，并将相对变化定义为：

Δ^obs = (V₁ - V₀)/V₀

Δ^obs的临床观察分布范围从-0.917（即减少91.7%）到扩大21.28倍；中位变化为0.0377。在75%的病例中，观察到的变化范围从42.2%缩小到约111%扩大。

使用表1中列出的同一组18个参数模拟每个患者和成像间隔。在每个模拟中，我们计算最终肿瘤细胞种群p + q位于区间I_c内的体素数量。定义模拟肿瘤体积的相对变化为：

Δ_i^sim = (V_i^sim - V₀)/V₀

Δ^sim的分布与观察到的肿瘤体积变化分布大致一致。表3显示了Δ^obs和Δ^sim落在选定百分比范围内的可能性（负百分比表示在下一次扫描时增强体素数量的相应减少）。

每个模拟产生18个肿瘤演化的实现（情景）。其中，令V^*是最接近下一次扫描时观察到的肿瘤体积V^obs的实现体积。我们定义：

R^best = (V^* - V^obs)/V^obs

在所有132个模拟时间间隔中，R^best的分布显示：在86%（65%）的病例中，相对误差小于20%（10%）。此外，18种参数选择的每一种都为某些患者成像间隔产生了"最佳"模拟肿瘤；表1中的组合没有一个是冗余的。

图2显示了一个代表性模拟结果。在这个例子中，实际肿瘤从1629个增强体素增长到2175个，即观察到的增长为33.5%。集成产生模拟肿瘤范围从1169到3495个体素；"最佳"参数集产生一个具有2167个体素的肿瘤，相对误差约为0.37%。

因为肿瘤的位置也很重要，我们定义了两个空间误差测量：令R表示实际患者扫描中分割为增强肿瘤的区域中的体素集合，S表示模拟肿瘤的体素集合。我们定义包含率C = |S∩R|/|R|和一致性A = |S∩R|/|S∪R|。

图5显示了本研究中进行的132次模拟中相对误差最大的模拟结果。在这种情况下，临床观察到的增强肿瘤体积显著缩小，从2245个体素减少到500个体素（Δ^obs = -0.773）。18个模拟肿瘤范围从2075到7246个增强体素。R^best的绝对值约为312%。肿瘤一致性和包含率测量范围从约0.02到0.18。

4 讨论

这项初步研究显示了简单数学模型在使用常规临床监测收集的影像数据生成个体化患者治疗反应情景方面的潜在效用。实验设计提供了患者肿瘤增长或缩小给定阈值的可能性的合理准确估计。在86%的研究时间间隔中，用于模拟每个肿瘤的18种参数中至少有一种选择产生的结果其体积与观察到的肿瘤体积相差在20%以内。

ASU模型相对于当代机器学习方法有三个优势：参数数量少（仅7个）；可以量化模型输出的不确定性；只需适度的计算资源。患者在我们的数据集中接受了各种治疗，包括放疗、化疗、抗血管生成药物和额外手术。除了手术（减少超出模型范围的肿瘤细胞群体），我们未尝试考虑治疗时间线。

未来的研究可能会根据选择的生物标志物对患者进行分层。例如，MGMT启动子甲基化状态和IDH突变可以影响GBM的治疗结果。修改的模型参数范围可能为选定的患者子集产生更临床相关的结果。

一旦给定患者的MR扫描上传到适当的数字存储库，可以应用第3.1节中概述的预处理工作流程。ASUB模型对分割边界的小变化不敏感；人工智能或其他半自动化方法可用于生成初步分割，只要在识别肿瘤区域没有严重错误。由此产生的预处理数据可以像任何其他影像文件一样下载，模拟程序本身可以打包为医生笔记本电脑上的"应用程序"。

配准错误是不可避免的，在某些情况下是显著的。配准错误的影响在切除腔附近尤其重要，肿瘤相关的大脑几何变化可能是一个影响因素。我们的模型模拟没有考虑质量效应，未来的工作可能受益于包含一个（但会涉及显著额外的计算复杂性）。

虽然对每位患者使用了T1和T2序列，但未遵循其他统一的成像协议。本研究使用的数据是在不同扫描仪上收集的。磁场强度不同，操作者定位患者的方式也不同。这些变化使图像比较复杂化，但它们也反映了"真实世界"的临床数据。

ASUB模型和FKPP模型假定肿瘤细胞有固定的承载能力。当以低于承载能力的细胞密度初始化时，两个模型都不能模拟减少的细胞群体。当增殖细胞在ASUB模型中变得足够多时，它们加入静止池。足够密集的静止细胞区域被假定不会在T1成像上增强（它们形成低信号区域）。因此，我们现有的模型参数化不能总是模拟增强体素数量的快速减少。

我们初始化肿瘤细胞群体和与患者扫描进行比较的方法简单，肿瘤分割仅基于相对像素强度。我们数据集中的一些患者接受了抗血管生成剂治疗，这可能会影响肿瘤在MRI上的外观（"假性反应"）。放疗可能导致"假性进展"，其中死亡细胞引起成像增强，难以与活跃生长的肿瘤细胞区分。灌注MRI的加入可能会改善我们对治疗效应的表征，并提供更好的模型初始化。

最后，通过比各向同性扩散更复杂的细胞运动模型可能会改善一致性和包含率测量。扩散张量成像（DTI）——无论是个体患者还是通过应用许多患者的平均DTI图——可能有助于胶质瘤进展的反应-扩散模型。类似的方法涉及粘性应力张量。此外，Fisher和ASU-Barrow模型的旅行波特性可能允许更好地初始化水肿区域中的肿瘤细胞群体。