随机微分方程与网络扩散模型：模拟阿尔茨海默病病理传播的不确定性

1. 项目概述：当数学遇见大脑，为阿尔茨海默病建模

作为一名长期在计算神经科学与生物统计交叉领域摸爬滚打的研究者，我常常思考一个问题：我们如何用冷冰冰的数学方程，去刻画像阿尔茨海默病（AD）这样复杂、多变且充满不确定性的生命过程？传统的确定性模型就像一张精确的列车时刻表，假设大脑这趟“疾病列车”会沿着固定轨道、准点到达每一个“病理车站”。但现实是，AD的进程充满了意外和个体差异，更像是在复杂地形中受风雨影响的徒步旅行。这正是我们这项工作的起点：将随机微分方程（SDE）的“不确定性”内核，与基于真实大脑连接数据的网络扩散模型相结合，试图构建一个更贴近现实、能“呼吸”的AD动态模拟器。

简单来说，这个项目的核心是用数学模拟蛋白质的“错误传播”。在AD患者的大脑中，两种错误折叠的蛋白——β淀粉样蛋白（Aβ）和tau蛋白——会像错误的谣言一样，沿着神经元之间的连接网络（即“连接组”）扩散、聚集，最终破坏认知功能。我们的模型以人类连接组计划（HCP）提供的真实大脑结构网络为“地图”，以Fisher-Kolmogorov方程描述蛋白质扩散和聚集的“基本交通规则”，再引入SDE来模拟生活中的各种“随机扰动”——比如昨晚没睡好、今天的饮食、长期的压力等无法精确测量的因素。最终，我们不仅想看疾病“平均”会怎么发展，更想量化这种发展的“可能性范围”，并通过贝叶斯推断从模拟数据中反推关键参数的分布。这对于理解为何不同患者疾病进展速度迥异，以及评估干预措施在个体层面的潜在效果，至关重要。

2. 核心思路与模型架构设计

2.1 为什么是“随机微分方程+网络扩散”？

在动手写一行代码之前，模型框架的选择决定了整个研究的视野和天花板。我们放弃了纯确定性模型，也并非完全转向黑箱机器学习，而是选择了“机理模型+随机过程”这条更具解释性的道路。

确定性模型的局限：经典的微分方程模型（如我们采用的Fisher-Kolmogorov方程）能优美地描述蛋白质扩散（∇· (D·∇c)）和局部聚集（αc[1-c]）这两个核心生物物理过程。但它有一个致命假设：只要初始条件和参数固定，结果唯一确定。这显然与临床观察相悖——同一年龄、相似基因型的患者，其认知衰退速度可能天差地别。

随机性的引入：因此，我们引入随机微分方程。在方程dc/dt = ...的右侧，我们增加了一个σc dW(t)项。这里的dW(t)就是著名的维纳过程（Wiener process），你可以把它想象成一系列随机的、无规律的“推搡力”。系数σ控制了推搡的力度。这个项直接模拟了所有未被模型显式包含的随机影响因素的总和，如微观层面的分子热运动、突触传递的随机性，以及宏观层面的生活方式波动等。这使得模型的每次模拟运行都像一次独立的“人生实验”，从而能生成一个可能结果的分布，而非一个单一轨迹。

网络作为传播骨架：蛋白质不会在均质的大脑中扩散。它们优先沿着轴突和突触连接传播，这已被多项病理学研究证实。因此，我们将大脑抽象为一个图G=(V, E)，其中129个脑区作为节点（V），脑区之间的白质纤维束连接强度作为边的权重（E）。这个网络结构通过拉普拉斯矩阵（Laplacian matrix）L被编码进扩散项。拉普拉斯矩阵是图论中描述扩散或平滑操作的核心工具，L作用于蛋白质浓度向量c，即-ΣL_ij c_j，精确刻画了蛋白质从高浓度节点向低浓度节点流动的趋势，流量由连接强度决定。

实操心得：模型复杂度的权衡在初期，我们曾考虑为Aβ和tau分别建模并引入相互作用项。但这会使得参数空间爆炸，且很多相互作用参数缺乏可靠的先验数据估计，反而可能引入更多噪声。因此，我们最终选择用一个综合浓度c来代表“病理负荷”，这是一个在保证模型可解性和计算可行性下的合理简化。关键在于，我们通过随机项σc dW(t)部分地囊括了这些未建模细节的集体效应。这种“聚合-随机化”策略在复杂系统建模中非常实用。

2.2 数据基石：人类连接组计划（HCP）与邻接矩阵构建

模型的血肉来自数学，但骨架必须来自真实数据。我们使用了人类连接组计划（HCP）公开的1064名健康成年人的弥散磁共振成像（dMRI）数据。处理流程如下：

1. 脑区分割与网络构建：使用标准的脑图谱（如Desikan-Killiany图谱）将每个个体的大脑皮层划分为129个感兴趣区域（ROI）。通过追踪dMRI数据中的白质纤维束，计算每两个脑区之间的纤维数量n_ij和平均纤维长度l_ij。

2. 加权邻接矩阵计算：连接强度并非简单的纤维计数。更短的纤维意味着更直接、更高效的连接。因此，我们采用了一个物理启发式的权重公式：A_ij = n_ij / (l_ij)^2.5。这里的指数2.5是一个经验参数，源于扩散张量成像研究，它强调了纤维长度对连接阻力的非线性影响（距离越长，有效连接强度衰减越快）。这个步骤将每个个体的大脑转化为一个129x129的加权邻接矩阵。

3. 群体平均模板生成：为了得到一个具有普适代表性的大脑结构网络模板，我们对1064个个体的邻接矩阵进行了元素级的平均。即，A_template[i,j] = mean(A_1[i,j], A_2[i,j], ..., A_1064[i,j])。这个平均矩阵（如图1所示）构成了我们所有模拟的静态结构基础。它代表了健康年轻人群体的大脑连接“平均蓝图”。

注意事项：使用群体模板的利与弊利：消除了个体差异的噪声，提供了一个稳定、可重复的计算基础，特别适合进行机制性、原理性的探索研究。弊：它抹杀了个体间连接组结构的特异性。而AD的起始和传播模式可能与个体特有的连接弱点（如默认模式网络）密切相关。因此，本模型更适用于揭示群体层面的统计规律。在向个性化预测迈进时，替换为个体特定的连接矩阵是必经之路。

2.3 模型方程的离散化与数值求解

有了网络L和方程理念，接下来就是将其转化为计算机可执行的形式。连续方程需要在离散的图节点上进行求解。

确定性部分（ODE）离散化： 原始的Fisher-Kolmogorov偏微分方程（PDE）∂c/∂t = ∇·(D∇c) + αc(1-c)在图上被离散化为一个常微分方程组（ODE）：dc_i/dt = - Σ_{j=1}^{N} L_ij * c_j + α * c_i * (1 - c_i)，对于i = 1, 2, ..., N。 这里，-Σ L_ij c_j项精确实现了在网络上的扩散，即节点i的浓度变化取决于它与所有邻居节点j的浓度差与连接强度的加权和。

随机部分（SDE）的加入： 将上述ODE转化为SDE：dc_i = [- Σ L_ij c_j + α c_i (1 - c_i)] dt + σ * c_i * dW_i(t)。 注意，我们为每个节点i都引入了一个独立的维纳过程W_i(t)。这意味着噪声是空间异质性的，不同脑区受到的随机扰动是相对独立的，这比使用一个全局共享噪声更符合生物学实际。

数值求解器选择： 对于确定性ODE部分，我们使用R语言中的deSolve包，它提供了高效可靠的数值积分器（如ode45类方法）。对于SDE，我们采用了Euler-Maruyama方法进行离散近似。这是求解SDE最经典的一阶方法，其迭代格式为：c_i^{n+1} = c_i^n + [-Σ L_ij c_j^n + α c_i^n (1 - c_i^n)] * Δt + σ * c_i^n * √(Δt) * Z_i^n。 其中，Z_i^n ~ N(0, 1)是标准正态分布的随机数，Δt是时间步长。我们通过大量测试，将Δt设为0.1（年），在数值稳定性和计算效率间取得了平衡。

3. 模拟实验设计与关键参数设置

3.1 初始条件与参数选择

一个模型的输出对其输入极为敏感，因此参数设置必须有据可依。

1. 病理种子位置：大量神经病理学证据表明，内嗅皮层是AD病理最早累及的区域之一，特别是tau蛋白的缠结。因此，我们设置模拟开始时，仅在内嗅皮层对应的节点上，初始病理浓度c = 0.1（即10%的蛋白发生错误折叠），其余所有节点c = 0。这模拟了疾病的“星星之火”。

2. 转化率常数 α：这个参数控制了蛋白质从正常形态转化为错误折叠形态的局部速率。我们参考了前驱研究，将其设为α = 0.5。这个值意味着在单个脑区内，病理转化过程具有一个中等的内在速率。在敏感性分析中，我们验证了在0.3-0.7范围内，疾病传播的时间尺度会相应缩放，但空间传播模式相对稳定。

3. 噪声强度 σ：这是本研究的核心变量之一。我们设计了两个对比场景：

   • 低噪声场景 (σ = 0.9)：模拟随机扰动较小的情形，可能对应遗传风险低、生活环境稳定的理想情况。
   • 高噪声场景 (σ = 2.0)：模拟随机扰动强烈的外在环境，可能对应多种风险因素叠加（如慢性疾病、长期睡眠障碍、高压力）的情况。 选择这两个值并非随意，而是通过预实验，使得σ=2时噪声的效应能够显著地、可视化地改变疾病的结局，从而便于对比分析。

4. 模拟时长与次数：我们将模拟时间设定为50年，覆盖了AD从临床前阶段到重度痴呆的典型时间跨度。对于每个噪声场景 (σ=0.9和σ=2.0)，我们分别运行1000次独立的随机模拟。这产生了两个各包含1000条可能疾病轨迹的集合，足以进行稳健的统计分析。

3.2 输出指标与分析方法

我们不仅关注“平均路径”，更关注“所有可能的路径”。

1. 节点与脑叶层面追踪：记录所有129个节点在每一个时间点的病理浓度c(t)。同时，根据脑图谱标签，将节点聚合到主要脑叶：额叶、颞叶、顶叶、枕叶，以及一个“其他”脑区组。计算每个脑叶内所有节点的平均浓度随时间的变化。

2. 确定性基线：首先运行无噪声的确定性模型 (σ=0)，得到一条“标准”疾病进展曲线，作为评估噪声影响的基准。

3. 随机模拟的可视化：

   • 时间序列区间图：绘制所有1000次模拟的脑叶平均浓度随时间变化的曲线，并用阴影区域表示其95%置信区间。这直观展示了疾病轨迹的“发散”程度。
   • 箱线图分析：在关键时间点（如第10、20、30年），绘制各脑叶病理浓度的箱线图，比较不同噪声强度下，浓度分布的均值、中位数、四分位距和异常值。

4. 贝叶斯推断：这是从模拟数据中“学习”的关键一步。我们聚焦于高噪声 (σ=2) 场景下，在特定时间点t（如第30年），某个脑叶（如额叶）的平均病理浓度μ_t。我们假设观测到的1000次模拟的浓度值服从正态分布N(μ_t, σ_t)。我们为均值μ_t设置一个无信息先验N(0, 10)，为标准差σ_t设置一个无信息先验Half-Cauchy(0, 2.5)。然后，使用马尔可夫链蒙特卡洛（MCMC）方法，通过近似贝叶斯计算（ABC）框架，从模拟数据中推断出μ_t和σ_t的后验分布。这告诉我们，基于我们的模型和模拟，在考虑了不确定性之后，我们对t=30年时额叶的平均病理浓度最合理的估计是什么，以及对这个估计的信心有多大。

4. 结果解析：噪声如何重塑疾病图景

4.1 确定性世界的“宿命论”路径

在没有噪声的确定性模型中，疾病展现出一条清晰、必然的路径（对应原文图3、4、5）。内嗅皮层作为“震中”，病理浓度迅速上升。疾病随后沿着连接网络有序传播。大约在第20年，全脑平均病理负荷超过50%。到第35年左右，几乎所有脑区的浓度都接近1（完全病变状态）。不同脑叶的易感性存在差异：颞叶（包含海马和内嗅皮层）受累最早最快，而额叶的进展最为迟缓。这与尸检研究和脑脊液生物标志物研究观察到的“颞叶-边缘系统优先”模式高度一致。这条路径是疾病的“平均力场”近似，但它描绘的是一幅过于确定、缺乏个体色彩的图景。

4.2 随机性引入后的“概率云”景观

一旦引入噪声，单一的疾病轨迹爆炸式地展开为一朵“概率云”。

低噪声 (σ=0.9) 场景：1000次模拟的轨迹彼此紧挨，95%置信区间较窄（如图6左）。虽然存在波动，但大多数轨迹仍在第40-45年左右汇聚到完全病变状态 (c=1)。噪声在此更像是对确定性路径的轻微“修饰”。

高噪声 (σ=2.0) 场景：情况截然不同（如图6右）。模拟轨迹从大约第15-20年开始剧烈分叉。关键的发现是：在模拟结束时（第50年），全脑平均病理浓度并未收敛于1，其均值大约在0.85-0.95之间，且存在一个显著的长尾分布。这意味着，在强随机扰动下，相当一部分模拟中，疾病并未达到理论上的“完全病变”状态。

脑叶异质性在噪声下的放大：箱线图（如图7）提供了更细致的快照。以第30年为例：

• 颞叶：即使在σ=2的高噪声下，其浓度分布依然较高且集中（中位数高，四分位距窄），说明其受累的“必然性”最强，受随机因素影响相对较小。
• 额叶：表现最为特殊。其浓度分布的中位数显著低于其他脑叶，且分布范围极广（箱体长，须线长），存在大量低浓度异常值。这表明额叶的病理进展不仅慢，而且其进程高度不可预测，极易受到随机因素的干扰而延迟。

核心发现解读：

1. “完美病变状态”的不可达性：高噪声下模型无法达到c=1，这具有深刻的生物学启示。它暗示在真实世界中，由于个体抵抗机制、环境缓冲等因素，AD可能很少以“全脑彻底沦陷”的极端形式存在，这与临床观察到的患者即使到晚期，某些功能仍可能部分保留的现象相符。
2. 疾病晚期的不确定性剧增：噪声的影响并非均匀。在疾病早期（前10-15年），所有模拟轨迹高度重叠，变异很小。这对应着临床前阶段，生物标志物虽已变化，但个体间差异不显。而中晚期（20年后）变异度急剧扩大，这精准地模拟了临床上的现象：从轻度认知障碍（MCI）向痴呆的转化时间，以及痴呆后的衰退速度，在患者间差异巨大，极难预测。
3. 额叶的“脆弱”与“韧性”：额叶进展最慢且最不稳定。慢，可能因为它距离病理播散的起点（内嗅皮层）网络距离较远，且其局部微环境对病理蛋白的“抵抗力”更强。不稳定，则可能因为它所负责的高级认知功能（执行功能、决策、社交）更易受到全身性随机因素（如全身炎症、代谢波动、睡眠质量日间差异）的影响。这为理解为何AD的行为和精神症状（BPSD）如此多变提供了计算视角。

4.3 贝叶斯推断：量化不确定性

对高噪声场景下t=30年时各脑叶平均浓度的贝叶斯推断，得到了收敛良好的后验分布（如图8）。以额叶为例，其后验均值μ约为0.55，95%最高密度区间（HPDI）为 [0.52, 0.58]。这个相对较窄的区间告诉我们：尽管单次模拟的结果千差万别（如图7中额叶的箱线图很宽），但当我们进行大量重复模拟后，其平均浓度的不确定性是可以被精确量化的。这个后验分布就是我们对“在如此高噪声环境下，疾病进展到30年时，额叶平均会是什么样子”这一问题的概率性回答。它比单一的均值报告包含了更丰富的信息。

5. 讨论、局限与未来方向

5.1 模型结果的临床与生物学关联

我们的计算发现与多项临床和生物学理论产生了有趣的共鸣：

• 网络退化假说：模型直观展示了病理如何沿结构连接网络传播，为这一假说提供了动力学模拟证据。
• 疾病异质性：模型生成的广泛结果分布，本身就是对AD巨大临床异质性的数学表征。噪声项σ可以视为个体一生中累积的“风险负荷”或“环境冲击”的抽象。
• 额叶功能保留：额叶进展慢且易受干扰的结果，或许部分解释了为何有些患者在记忆严重受损后，人格和部分执行功能仍能相对保留一段时间。

5.2 当前模型的局限性

没有完美的模型，只有不断改进的模型。我们必须清醒认识本工作的局限：

1. 蛋白种类的聚合：将Aβ和tau聚合为单一变量c，丢失了它们之间复杂的、有时序先后的相互作用（如Aβ驱动tau病理）。这是模型最大的生物学简化。
2. 静态的连接网络：我们使用了健康的、平均化的连接组。而实际上，AD病理本身会破坏突触和连接，导致网络结构随时间动态衰退。这是一个“病理影响网络，网络又影响病理传播”的耦合过程，当前模型未包含这一反馈。
3. 噪声的均质化假设：我们假设所有脑区受到噪声的强度 (σ) 相同。实际上，不同脑区对氧化应激、炎症等因子的敏感性可能不同。
4. 参数的经验性：转化率α、噪声强度σ等关键参数仍需更多来自纵向生物标志物研究的实证数据来校准和约束。

5.3 未来可拓展的方向

这项工作是起点，而非终点。基于此框架，至少有以下几个有潜力的拓展方向：

1. 双蛋白耦合模型：建立两个相互耦合的SDE系统，分别描述Aβ和tau的浓度，并引入Aβ促进tau磷酸化的项。这能模拟更经典的“淀粉样蛋白级联假说”动态。
2. 个性化预测：将群体平均连接矩阵替换为个体患者的弥散张量成像（DTI）连接组。同时，尝试将个体的遗传风险（如APOE ε4等位基因数量）、基线认知分数等作为先验信息，调整α或σ的初始分布，实现真正的个性化疾病轨迹风险预测。
3. 干预模拟：在方程中引入治疗项。例如，假设一种疗法能全局降低病理浓度（如增加清除率），或阻断特定脑区间的传播（如修改邻接矩阵A的特定元素），然后运行模拟，观察其对疾病轨迹概率分布的影响，从而在计算机上评估不同治疗策略的潜在效果。
4. 多模态数据融合：利用淀粉PET、tau PET、fMRI功能连接等多模态影像数据，不仅校准模型初始状态，还可以在模拟过程中引入动态约束，使模型输出与患者纵向的多模态数据变化相匹配。

最后一点个人体会：构建这样一个模型，最深刻的感受是“谦逊”。我们并非在预测某个特定患者大脑中蛋白质分子明天的确切位置，而是在探索疾病在千万种可能性中演进的概率地形图。随机微分方程的魅力在于，它诚实地承认了世界的复杂性和我们的无知（用σ dW(t)来概括），并以此为基础，生成一幅充满可能性的图景。这幅图景或许无法给出确切的答案，但它能更好地提出正确的问题：为什么有些脑区更脆弱？哪些随机因素影响最大？在哪个时间窗口进行干预，最有可能将疾病轨迹推向更好的方向？这才是计算模型在神经退行性疾病研究中真正的价值——一个用于整合知识、生成假设、量化不确定性的强大思维实验平台。