简介：本资源是一套面向神经工程与计算建模方向的MATLAB实践材料，聚焦丘脑深部脑刺激（DBS）的网络机制解析，适用于计算机、电子信息工程、应用数学等专业本科生开展课程设计、期末大作业及毕业设计。压缩包共98个文件，含58个.mat实验数据集（存储神经动力学仿真结果）、17个.m主程序与函数脚本（实现率编码网络模型、特征分析与路径优化）、10个.txt说明文档（含Readme与模型参数注解），辅以CSV数据、FIG/PNG可视化图、XLSX统计表及少量Python与Word辅助文件，整体体积11.18MB，结构模块化、层次清晰。已有81人学习下载。用户可直接运行全部代码，无需额外配置；参数化设计支持快速调整刺激强度、连接权重等关键变量；代码注释详尽、逻辑分层明确，并配套特征分析（如特征值对分析）、模型对比（尖峰网络vs速率网络）及可视化流程，显著降低神经建模入门门槛。
1. 这不是普通神经建模代码包：它用率编码网络（rate network model）+ 特征模态分析（eigenpairs）解构丘脑DBS的全脑传播路径
如果你正在做神经工程方向的课程设计，手头却只有单神经元放电图（spiking raster）和零散LFP频谱——那这个 ZIP 包里的 MATLAB 实现会直接改写你的工作流。它不模拟单个神经元的动作电位，而是构建一个参数化率编码网络模型（rate network model），把丘脑核团、皮层分区、基底节节点抽象为动态变量，再通过雅可比矩阵的特征模态（eigenpairs）定位刺激信号在全脑网络中的主导传播方向与衰减模式。实验数据部分包含真实大鼠/非人灵长类的多通道记录片段（.mat 格式），与模型输出严格对齐；compare_with_spiking_network_model子目录则提供轻量级脉冲网络对照组，用于验证率模型在低计算开销下对宏观动力学的保真度。适用场景非常明确：电子信息工程专业做“生物医学信号建模”课设、数学系做“动力系统稳定性分析”大作业、或神经科学方向本科生开展 DBS 机制初探——所有代码均基于 MATLAB 2014a 及以上版本编写，注释覆盖每个参数物理意义（如tau_e = 15; % ms, excitatory synaptic time constant），无需修改即可运行analysis_with_eigenpairs.m得到前 3 阶特征向量的空间分布热图。
2. 率编码网络建模原理与 MATLAB 实现：从微分方程到可调参状态空间
2.1 为什么选择率编码而非脉冲模型？——计算效率与可解释性的平衡点
在丘脑深部脑刺激（DBS）机制研究中，全尺度脉冲网络（spiking network model）虽能复现单细胞精度，但其计算复杂度随神经元数量呈超线性增长（O(N²) 突触更新）。本项目采用Wilson-Cowan 类型率编码模型，将每个脑区视为一个平均发放率变量 rᵢ(t)，其演化由以下常微分方程描述：
$$ \tau_i \frac{dr_i}{dt} = -r_i + F_i\left( \sum_j w_{ij} r_j + I_i^{\text{ext}} \right) $$
其中 $F_i$ 是 Sigmoid 型增益函数，$w_{ij}$ 为结构连接权重（来自公开的 Macaque connectome 数据集），$\tau_i$ 为时间常数。该形式将 10⁴ 级别神经元的仿真压缩至 10–50 个区域变量，且特征值分析可直接映射到功能连接梯度上。rate_network_model-main/目录下的build_network.m负责加载预置连接矩阵并施加 DBS 刺激项（I_i^{\text{ext}} = I_0 \cdot \delta_{i,\text{thal}} \cdot \sin(2\pi f t)），这是后续所有分析的起点。
提示：
build_network.m中I_0默认设为0.8（归一化强度），若需匹配特定实验刺激幅值（如 100 μA），需按比例缩放——查看data/experimental/README.md中的stimulation_amplitude_uA字段进行换算。
2.2 核心建模脚本结构解析与关键参数表
整个率网络模型由四个主脚本协同驱动，全部位于rate_network_model-main/下：
| 脚本名 | 功能 | 关键可调参数（含默认值） | 修改影响说明 |
|---|---|---|---|
build_network.m | 初始化连接矩阵、设定刺激靶点与参数 | region_list = {'Thal','PFC','M1','STN','GPi'};tau_e = 15; tau_i = 10; | 更改region_list会重定义状态变量维度；tau_*影响响应速度与振荡频率边界 |
simulate_dynamics.m | 数值求解 ODE（使用ode15s），输出时间序列r_all | tspan = [0 2000]; dt = 0.1;solver_opts = odeset('RelTol',1e-5); | tspan决定仿真总时长（ms）；dt过大会导致高频成分失真，建议 ≤0.5 ms |
compute_jacobian.m | 在稳态点计算雅可比矩阵 J，并返回特征值/向量 | eq_point = find_equilibrium(r_all);J = jacobian_numerical(...); | 若系统存在多稳态，eq_point需手动指定索引（见README.txt第7行说明） |
analysis_with_eigenpairs.m | 主分析入口：绘制特征向量空间分布、计算模态参与比（MPR） | n_eig = 3;threshold_mpr = 0.15; | n_eig控制提取前 N 阶模态；threshold_mpr用于筛选高参与度脑区（见 3.2 节） |
2.2.1 手动验证稳态点的必要性
由于 DBS 刺激可能诱导多稳态，compute_jacobian.m中的find_equilibrium函数并非全自动。需先运行simulate_dynamics.m获取r_all，再执行：
% 在 MATLAB 命令窗口中交互式定位稳态 figure; plot(r_all(1,:)); xlabel('Time (ms)'); ylabel('r_{Thal}'); grid on; % 观察曲线后 500 ms 是否进入平台期（如图中 1500–2000 ms 段） eq_idx = 15000; % 对应 t=1500 ms 的索引（因 dt=0.1 ms → 1500/0.1=15000） r_eq = r_all(:, eq_idx); % 提取该时刻各区域发放率此步骤确保雅可比矩阵在生理相关工作点处计算，避免特征分析结果漂移。
2.3 连接权重矩阵的加载与校准逻辑
build_network.m通过load('data/connectome_weights.mat')加载预存连接矩阵W_struct（尺寸：N×N，N=region_list 长度）。但真实 DBS 效应不仅取决于解剖连接，还受白质纤维各向异性影响。因此代码中嵌入了有效连接缩放因子：
% build_network.m 第 42 行起 W_effective = W_struct; if ~isempty(stim_target_region) && isfield(data_config,'dti_scaling') % 对刺激靶点（如 Thal）的传出连接应用 DTI 各向异性权重 idx_thal = find(strcmp(region_list, stim_target_region)); W_effective(idx_thal,:) = W_struct(idx_thal,:) .* data_config.dti_scaling; enddata_config.dti_scaling来自data/experimental/config_dti.mat，其值为[0.92, 1.05, 0.88, ...]—— 这些数值源自 Diffusion Tensor Imaging 实测的纤维密度归一化结果。若替换为自己的 DTI 数据，只需保证dti_scaling向量长度等于region_list长度，并重新保存为.mat文件即可。
3. 特征模态（eigenpairs）分析实战：从雅可比矩阵到功能通路可视化
3.1 特征值物理意义与稳定性判据
compute_jacobian.m输出的特征值 λₖ = αₖ + iβₖ 直接决定网络动力学行为：
- 实部 αₖ：反映模态衰减/增长速率。若所有 αₖ < 0，则稳态点渐近稳定；最大实部 αₘₐₓ > 0 表明存在自发振荡。
- 虚部 βₖ：对应振荡频率 fₖ = βₖ/(2π)（Hz）。DBS 常用 130 Hz，故需检查是否存在 βₖ ≈ 2π×130 的特征值。
- 特征向量 vₖ：其第 i 个分量 |vₖᵢ|² 表示第 i 个脑区对该模态的参与度（Participation Ratio），是空间分布可视化的基础。
运行analysis_with_eigenpairs.m后，控制台输出类似：
Eigenvalue analysis completed: - Dominant mode: λ₁ = -0.023 + 816.2i → f ≈ 130.0 Hz (matches DBS frequency) - MPR(λ₁) = 0.42 → Thal(0.31), PFC(0.28), M1(0.22) are top contributors这说明第一阶模态主导 130 Hz 振荡，且丘脑-前额叶-运动皮层构成核心环路。
3.2 绘制特征向量空间分布热图
analysis_with_eigenpairs.m调用plot_eigenvector_spatial.m生成热图。关键代码段如下：
% analysis_with_eigenpairs.m 第 89 行 figure('Position',[100,100,800,600]); subplot(1,2,1); bar(abs(eigvec(:,1)).^2, 'FaceColor', [0.2 0.6 0.8]); xticklabels(region_list); xtickangle(45); title('Participation Ratio of Mode 1 (\lambda_1)'); ylabel('PR_i = |v_{i1}|^2'); subplot(1,2,2); % 使用地理坐标映射（需提前加载 atlas_coords.mat） load('data/atlas_coords.mat'); % 包含 x,y,z 坐标及标签 scatter3(coords(:,1), coords(:,2), coords(:,3), 120, abs(eigvec(:,1)).^2, 'filled'); colorbar; title('Spatial Distribution of Mode 1');注意：
atlas_coords.mat中的坐标系为 MNI152 标准空间，单位 mm。若使用其他脑图谱（如 AAL3），需替换该文件并确保coords行数与region_list一致。
3.3 模态参与比（MPR）阈值筛选与通路提取
为定量识别“DBS 主导通路”，代码引入模态参与比（Modal Participation Ratio）：
$$ \text{MPR}k = \frac{ \left( \sum_i |v{ki}|^2 \right)^2 }{ \sum_i |v_{ki}|^4 } $$
MPR ∈ [1, N]，值越大表示能量越分散；MPR ≈ 1 表示仅 1–2 个脑区主导。analysis_with_eigenpairs.m中设置threshold_mpr = 0.15，即筛选满足|vₖᵢ|² > threshold_mpr × max(|vₖ|²)的脑区。实际执行逻辑为：
% analysis_with_eigenpairs.m 第 112 行 v1_abs2 = abs(eigvec(:,1)).^2; v1_max = max(v1_abs2); core_regions = region_list(v1_abs2 > threshold_mpr * v1_max); fprintf('Core regions for Mode 1: %s\n', strjoin(core_regions, ', ')); % 输出：Core regions for Mode 1: Thal, PFC, M1此结果可直接用于论文图 3 的“DBS 有效传播通路”示意图。
4. 与脉冲网络模型（spiking network model）的定量对比方法
4.1 对照实验设计：双模型同构输入与输出对齐策略
compare_with_spiking_network_model/目录提供轻量级 Izhikevich 脉冲网络实现，其拓扑结构、连接权重、外部输入完全复刻率模型。对比关键在于输出量纲统一：
- 率模型输出：各区域平均发放率
r_i(t)（Hz） - 脉冲模型输出：各区域单位时间脉冲计数
n_i(t)/Δt（Hz）
compare_models.m脚本自动完成对齐：
% compare_models.m 第 35 行 % 将脉冲计数转换为等效发放率（滑动窗 10 ms） spike_rate = zeros(N, T); for i = 1:N spike_times = spk_data{i}; % {cell array of spike times for region i} spike_rate(i,:) = histcounts(spike_times, 0:dt:T*dt, 'Normalization','pdf') * (1/dt); end % 此时 spike_rate(i,t) 单位为 Hz，与 r_all(i,t) 完全可比4.2 三种量化对比指标与 MATLAB 实现
对比不依赖主观观察，而采用三个鲁棒指标：
| 指标 | 计算公式 | MATLAB 代码片段 | 解读 |
|---|---|---|---|
| 时间域相关性 | corr(r_i, spike_rate_i, 'rows','complete') | corr_coef = corr(r_all(1,:), spike_rate(1,:)); | 值越接近 1，说明两模型在该区域动态轨迹一致性越高 |
| 功率谱密度（PSD）重叠度 | 1 - norm(psd_rate - psd_spike,'fro') / norm(psd_rate,'fro') | psd_rate = pwelch(r_all(1,:),[],[],[],fs); | 重叠度 >0.85 视为频谱特性高度一致 |
| 相位同步指数（PLV） | abs(mean(exp(1j*(phi_rate - phi_spike)))) | phi_rate = angle(hilbert(r_all(1,:))); | PLV ∈ [0,1]，>0.7 表示锁相关系强 |
运行compare_models.m后生成comparison_report.pdf，内含三组指标表格及 PSD 重叠图。典型结果：Thal-PFC 通路的 PLV 达 0.79，证实率模型在宏观同步性上具备足够保真度。
5. 排查 ZIP 解压与 MATLAB 运行常见故障：从 invalid zip archive 到特征值 NaN
5.1 ZIP 包完整性验证与解压异常处理
下载的“揭示丘脑深部脑刺激的网络机制….zip文件名含中文与省略号，易导致解压工具识别失败。必须重命名为英文无空格名称（如thalamus_dbs_model.zip）后再解压。若遇invalid zip archive: could not find eocd错误：
# Linux/macOS 终端检查 ZIP 结构 file thalamus_dbs_model.zip # 应返回：thalamus_dbs_model.zip: Zip archive data, at least v2.0 to extract # 若返回 "data" 或 "cannot open"，说明下载不完整，需重新获取 unzip -t thalamus_dbs_model.zip # 测试压缩包完整性Windows 用户请使用 7-Zip（非系统自带解压器），并在设置中勾选“使用 UTF-8 编码读取文件名”（防止中文路径乱码）。
5.2 MATLAB 运行时报错定位与修复方案
错误 1：Undefined function or variable 'eig'
原因：未安装 MATLABSymbolic Math Toolbox（eig函数依赖此工具箱）
修复：在 MATLAB 命令窗口执行ver查看已安装工具箱，若无Symbolic Math Toolbox，需通过Add-Ons → Get Add-Ons安装。
错误 2：Error in compute_jacobian (line 67): J = jacobian_numerical(...)报NaN
原因：r_eq中存在Inf或NaN，导致数值微分失效
诊断：
% 运行 simulate_dynamics.m 后立即执行 disp(['r_eq contains Inf: ', num2str(any(isinf(r_eq)))]); disp(['r_eq contains NaN: ', num2str(any(isnan(r_eq)))]); % 若为 true，检查 build_network.m 中的增益函数 F_i 是否饱和（如 Sigmoid 输入过大）修复：在build_network.m中限制输入范围：
% 替换原 Sigmoid 行（约第 85 行） % x = sum(w_ij * r_j) + I_ext; x = max(-20, min(20, sum(w_ij * r_j) + I_ext)); % 截断至 [-20,20] r_i = 1 ./ (1 + exp(-(x - theta)/sigma));错误 3：analysis_with_eigenpairs.m绘图空白或坐标轴错乱
原因：data/atlas_coords.mat缺失或格式错误
验证：
load('data/atlas_coords.mat'); whos coords % 应显示 size: N×3, class: double assert(size(coords,2)==3, 'atlas_coords must have 3 columns (x,y,z)');若缺失，从data/experimental/中复制atlas_coords_template.mat并按实际脑区顺序填充坐标。
5.3 快速验证模型是否正常工作的三步法
无需运行全部脚本，用以下命令链 60 秒内确认核心功能：
% 步骤 1：检查数据加载 cd rate_network_model-main; load('../data/experimental/stim_data.mat'); disp(['Stimulation data loaded: ', num2str(size(stim_data,1)), ' time points']); % 步骤 2：生成单步雅可比矩阵（跳过ODE求解） r_test = ones(length(region_list),1) * 0.5; % 人工设定测试点 J_test = compute_jacobian(r_test, @build_network); disp(['Jacobian size: ', num2str(size(J_test,1)), 'x', num2str(size(J_test,2))]); % 步骤 3：提取并显示首阶特征向量 [eigvec_test, eigval_test] = eig(J_test); disp(['Dominant eigenvalue: ', num2str(eigval_test(1,1))]); disp(['Top contributor: ', region_list{find(abs(eigvec_test(:,1))==max(abs(eigvec_test(:,1))))}]);若输出显示Dominant eigenvalue为复数且Top contributor为Thal，则模型环境已就绪，可进入完整流程。
本文还有配套的精品资源，点击获取
