心血管模型比较与HRV相干性分析

本项目属于 Master ATSI 研究团队项目（2025-2026学年），目标是通过计算机仿真研究心肺系统的耦合关系，并分析心率变异性（HRV）的相干性特征。

生理学背景

心脏的跳动并非由心脏自身完全决定，而是受到自主神经系统（Autonomic Nervous System, ANS）的持续调控。ANS 分为两个功能相反的分支：交感神经和副交感神经。交感神经在人体面临压力或危险时被激活，触发所谓的"战斗或逃跑"反应，表现为心率加快、血压升高、呼吸加速等；副交感神经则在人体处于放松状态时占主导，促进"休息和消化"功能，使心率减慢、血压降低。

呼吸活动会周期性地影响这两个神经分支的活动强度。具体而言，吸气时交感神经活动增强，心率略微加快；呼气时副交感神经活动增强，心率略微减慢。这种心率随呼吸节律波动的现象在生理学上称为呼吸性窦性心律不齐（Respiratory Sinus Arrhythmia, RSA）。RSA 的存在说明心血管系统和呼吸系统之间存在天然的耦合关系。

心脏相干性呼吸

心脏相干性呼吸是一种特定的呼吸调节技术，要求以每分钟 6 次的频率进行缓慢而规律的呼吸，对应的呼吸频率为 0.1 Hz，即每个呼吸周期持续 10 秒。选择这个频率并非随意，而是因为人体心血管系统对 0.1 Hz 附近的输入具有共振特性，在此频率下呼吸能够最大限度地放大心率的波动幅度，从而最大化心率变异性（Heart Rate Variability, HRV）。研究表明，较高的 HRV 与良好的心血管健康、情绪调节能力以及压力适应能力相关联。

仿真管线架构

本项目的核心是构建一条完整的信号处理链，模拟从呼吸输入到 HRV 输出的全过程。整条管线可以用以下流程表示：

这条管线的物理含义如下：首先，呼吸波形 V_t(t) 代表潮气量随时间的变化，描述了呼吸的节律和深度。呼吸机械模型将这个机械信号转换为神经信号 MN(t)，代表从脑干（延髓和脑桥）发出的运动神经元活动。神经信号随后输入心血管模型，该模型模拟自主神经系统对心脏的调控作用，输出心率 HR 和血压 BP 的时间序列。最后，对心率信号进行 HRV 分析，提取各种时域和频域指标，评估心率变异性的特征和相干程度。

代码模块结构

整个项目的代码组织为 7 个 Cell，每个 Cell 对应仿真管线中的一个功能模块。Cell 1 负责导入所需的 Python 库并完成环境初始化。Cell 2 定义了一系列工具函数，包括时间轴生成、信号滤波、标准化和绑图等基础功能。Cell 3 实现呼吸波形的生成，支持正弦波、方波、矩形波等多种波形类型。Cell 4 包含两个呼吸机械模型，将潮气量转换为神经信号。Cell 5 定义心血管模型，将神经信号转换为心率和血压。Cell 6 实现 HRV 分析功能，包括 RR 间期提取、时域指标计算和频域相干性分析。Cell 7 整合同步性指标计算并提供交互式 GUI 界面。

运行环境为 Google Colab，仅需 CPU 即可完成全部计算，按照 Cell 1 到 Cell 7 的顺序依次执行。

Cell 1：依赖库导入与环境初始化

本 Cell 是整个项目的基础，导入后续所有模块所依赖的 Python 库。运行顺序上必须首先执行此 Cell，否则后续代码会因找不到相应模块而报错。

数值计算与数据处理库

NumPy 是 Python 科学计算的核心库，提供了高效的多维数组对象 np.ndarray 以及大量的数学运算函数。本项目中所有的时间序列数据（呼吸波形、心率、血压等）都以 NumPy 数组的形式存储和处理。NumPy 的底层用 C 语言实现，对数组的批量运算远比 Python 原生循环快得多。

Pandas 建立在 NumPy 之上，提供了 DataFrame 数据结构。DataFrame 类似于电子表格，每一列可以有不同的数据类型，并且支持按行名或列名索引。在本项目中，仿真过程产生的状态变量（如心率、血管张力、血压随时间的变化）被组织成 DataFrame，便于查看、分析和导出为 CSV 文件。

类型提示与数据容器

dataclass 是 Python 3.7 引入的装饰器，用于简化数据容器类的定义。普通的 Python 类需要手动编写 __init__ 方法来初始化属性，而使用 @dataclass 装饰器后，只需声明属性及其类型，Python 会自动生成构造函数、__repr__ 方法等。本项目定义了 SimRes 类来统一封装所有模型的仿真结果，包括时间轴、输出信号、内部状态和元信息。

typing 模块提供类型提示功能。例如 Dict[str, np.ndarray] 表示一个字典，其键是字符串类型，值是 NumPy 数组类型。类型提示不会影响代码的实际运行（Python 仍然是动态类型语言），但能让阅读代码的人快速理解函数参数和返回值的预期类型，也便于 IDE 提供智能补全和错误检查。

常微分方程求解器

solve_ivp 是 SciPy 提供的常微分方程（Ordinary Differential Equation, ODE）初值问题求解器。心血管模型的核心是一组描述生理状态动态变化的微分方程，其一般形式为：

其中 y 是状态向量（可能包含心率、血管张力、血压等多个变量），f 是描述系统动力学的函数。给定初始条件 y(t_0) = y_0 和求解的时间范围 [t_0, t_1]，solve_ivp 通过数值积分方法逐步推进，计算出 y(t) 在各个时刻的数值解。

solve_ivp 支持多种数值积分算法。RK45 是默认算法，基于经典的四五阶 Runge-Kutta 方法，通过比较四阶和五阶估计的差异来自适应调整步长，适用于大多数非刚性问题。Radau 是一种隐式 Runge-Kutta 方法，适合求解刚性问题（即系统中存在变化速度相差悬殊的分量）。BDF（Backward Differentiation Formula，向后差分公式）也是处理刚性问题的常用方法。

信号处理函数

welch 函数用于估计信号的功率谱密度（Power Spectral Density, PSD）。直接对信号做快速傅里叶变换（FFT）得到的频谱往往噪声较大，Welch 方法通过将信号分成多个重叠的段，对每段分别计算周期图，然后取平均，从而得到更平滑、更稳定的功率谱估计。在 HRV 分析中，Welch 方法用于计算 RR 间期序列的功率谱，进而分析不同频段的能量分布。

butter 函数用于设计 Butterworth 滤波器。Butterworth 滤波器的特点是在通带内具有最平坦的幅频响应，没有纹波，因此又称为最大平坦滤波器。设计滤波器时需要指定阶数和截止频率，函数返回滤波器的系数。

filtfilt 实现零相位滤波。普通的数字滤波器在滤除特定频率成分的同时会引入相位延迟，导致滤波后的信号在时间上发生偏移。filtfilt 的做法是先正向滤波一次，再将结果反转后反向滤波一次，两次滤波引入的相位延迟正好抵消，最终得到零相位延迟的结果。这种方法只能用于离线处理（需要完整的信号数据），不能用于实时滤波。

hilbert 函数计算信号的 Hilbert 变换，将实数信号转换为解析信号（复数形式）。解析信号的形式为 z(t) = x(t) + j\hat{x}(t)，其中 x(t) 是原始实信号，\hat{x}(t) 是其 Hilbert 变换。从解析信号可以方便地提取瞬时幅度 |z(t)| 和瞬时相位 \angle z(t)。在本项目中，Hilbert 变换用于计算两个信号之间的相位锁定值（Phase Locking Value, PLV），评估呼吸与心率之间的同步程度。

可视化与交互界面

Matplotlib 是 Python 最常用的绑图库，提供了类似 MATLAB 的绑图接口。本项目使用它来绑制呼吸波形、心率和血压的时间序列、功率谱密度图等。

ipywidgets 在 Jupyter Notebook 或 Google Colab 环境中提供交互式 GUI 组件，包括滑块（Slider）、按钮（Button）、下拉框（Dropdown）等。通过这些组件，用户可以在不修改代码的情况下调整仿真参数（如呼吸频率、仿真时长），并实时查看结果变化。gui_ok = True 是一个标志位，表示 GUI 组件可用，后续代码会根据此标志决定是否创建交互界面。

环境配置

np.random.seed(7) 固定了 NumPy 的随机数种子。虽然本项目的仿真过程是确定性的，不涉及随机数生成，但固定种子是一个良好的习惯，确保在任何涉及随机性的场合下结果都可复现。

plt.rcParams 用于配置 Matplotlib 的全局默认参数。设置图像尺寸为 11×7 英寸，并默认在所有子图中显示网格线，使得绘制的图表更易于阅读和分析。

Cell 2：工具函数

本 Cell 定义了一组基础工具函数，这些函数本身不执行任何仿真逻辑，而是作为底层支持被后续的呼吸模型、心血管模型和 HRV 分析模块反复调用。将这些通用功能抽取为独立函数，既避免了代码重复，也使得各个模块的主体逻辑更加清晰。

时间轴生成函数

gen_taxis 函数用于生成等间隔的时间轴数组。仿真的第一步是确定时间范围和时间分辨率，输入参数包括总时长（单位：秒）和采样频率（单位：Hz，即每秒采样点数），输出是一个从 0 到总时长的等间隔一维数组。

采样点数的计算方式为总时长乘以采样频率再加 1（包含起点和终点）。例如，若总时长为 10 秒，采样频率为 100 Hz，则生成 1001 个时间点，相邻点之间的间隔为 0.01 秒。函数内部使用 np.linspace 来生成这个等间隔序列，确保首尾两端的值精确等于 0 和总时长。

边缘平滑函数

smooth_edge 函数用于平滑信号中的跳变边缘。理想的方波或矩形波在高低电平之间瞬间跳变，但真实的生理信号不可能有如此剧烈的突变。为了让生成的呼吸波形更接近真实情况，需要对这些尖锐的边缘进行平滑处理。

函数采用双曲正切函数 \tanh 来实现平滑。对于输入信号 x 和平滑参数 e，输出为：

当 e 很小时，\tanh(x/e) 在 x=0 附近的过渡非常陡峭，接近阶跃函数；当 e 较大时，过渡区域变宽，曲线更加平缓。通过调节 e 的值，可以控制边缘平滑的程度。\tanh 函数的输出范围在 (-1, 1) 之间，天然具有饱和特性，不会因为输入值过大而产生数值溢出。

低通滤波函数

lpf 函数对信号进行低通滤波，保留低频成分而衰减高频成分。在信号处理中，低通滤波常用于去除高频噪声或提取信号的缓变趋势。

采样定理与 Nyquist 频率

理解低通滤波器的设计需要先掌握采样定理。根据 Nyquist-Shannon 采样定理，若要从离散采样点无失真地重建连续信号，采样频率 f_s 必须至少是信号最高频率分量的 2 倍。换言之，给定采样频率 f_s，能够正确表示的最高频率为：

这个频率称为 Nyquist 频率。任何高于 f_{nyq} 的频率分量在采样后会发生混叠（aliasing），被错误地映射到低频区域，导致信号失真。

在设计数字滤波器时，截止频率通常需要以 Nyquist 频率为基准进行归一化。scipy.signal.butter 函数要求截止频率参数在 0 到 1 之间，其中 1 对应 Nyquist 频率。因此，若实际截止频率为 f_c，采样频率为 f_s，传入 butter 的归一化截止频率应为 f_c / f_{nyq} = 2f_c / f_s。

函数内部首先调用 butter 设计指定阶数的 Butterworth 低通滤波器，获得滤波器系数 b 和 a，然后调用 filtfilt 进行零相位滤波。零相位滤波通过正向和反向各滤波一次来消除相位延迟，适合离线信号处理场景。

z-score 标准化函数

znorm 函数对信号进行 z-score 标准化，将其变换为零均值、单位方差的形式。标准化公式为：

其中 \mu 是信号 x 的均值，\sigma 是标准差。标准化后的信号 z 均值为 0，标准差为 1。

标准化的意义在于消除不同信号之间的尺度差异。例如，呼吸模型 A 输出的神经信号可能在 2 到 14 的范围内波动，而模型 B 的输出可能在 7 到 17 的范围内。如果心血管模型直接使用这些原始数值，模型参数就需要针对不同的输入范围分别调整。通过 z-score 标准化，无论输入信号的原始范围如何，标准化后都具有相同的统计特性，心血管模型的参数可以保持通用。

函数实现中，分母加上一个极小值 10^{-12} 是为了防止当信号为常数（标准差为零）时发生除零错误。

多信号绑图函数

plot_sigs 函数用于在一张图中绑制多个时间序列，每个信号占据一个垂直排列的子图。输入参数包括时间轴数组和一个信号字典，字典的键是信号名称（将显示为 y 轴标签），值是对应的信号数组。

函数根据字典中信号的数量创建相应数量的子图，设置 sharex=True 使所有子图共享同一个 x 轴，便于在时间上对齐比较不同信号的变化。只有最下方的子图显示 x 轴标签 "Time (s)"，避免重复标注。这种堆叠式的绑图方式非常适合展示仿真管线中各阶段信号的演变过程。

Cell 2：工具函数代码实现

本部分对 Cell 2 中各工具函数的代码实现进行逐行解析，重点说明每行代码的作用以及实现背后的技术考量。

时间轴生成函数 gen_taxis

def gen_taxis(dur: float, fs: float) -> np.ndarray:
    n_pts = int(np.floor(dur * fs)) + 1
    return np.linspace(0.0, dur, n_pts)

函数签名中的类型提示 dur: float 和 fs: float 表明两个参数都应为浮点数，-> np.ndarray 表明返回值是 NumPy 数组。参数 dur 代表仿真总时长（秒），fs 代表采样频率（Hz）。

第一行计算采样点数。dur * fs 得到总时长内包含的采样间隔数，np.floor 向下取整确保不超出时间范围，int() 将结果转为整数，最后加 1 是因为 n 个间隔对应 n+1 个端点（包含起点和终点）。例如，时长 10 秒、采样率 100 Hz 时，有 1000 个间隔，对应 1001 个采样点。

第二行使用 np.linspace(0.0, dur, n_pts) 生成从 0 到 dur 的等间隔数组。与 np.arange 不同，linspace 直接指定点数而非步长，能确保首尾两端的值精确等于指定值，避免浮点累积误差。

边缘平滑函数 smooth_edge

def smooth_edge(x: np.ndarray, e: float = 0.02) -> np.ndarray:
    return np.tanh(x / max(e, 1e-6))

参数 e 的默认值设为 0.02，这是一个经验值，在大多数场景下能产生适度的平滑效果。max(e, 1e-6) 确保分母不会为零或负数，即使调用者错误地传入了 0 或负值，函数也不会崩溃。

np.tanh 是双曲正切函数，其数学定义为：

该函数的输出范围严格在 (-1, 1) 之间，当 u \to +\infty 时输出趋近于 1，当 u \to -\infty 时输出趋近于 -1，当 u = 0 时输出为 0。将输入 x 除以 e 后，e 越小则 x/e 的变化越剧烈，\tanh 曲线越陡峭，接近阶跃；e 越大则曲线越平缓，过渡越圆滑。

低通滤波函数 lpf

def lpf(x: np.ndarray, fs: float, fc: float, n_ord: int = 3) -> np.ndarray:
    nyq = 0.5 * fs
    b, a = butter(n_ord, fc / nyq, btype="low")
    return filtfilt(b, a, x)

参数 fc 是截止频率（Hz），n_ord 是滤波器阶数，默认为 3。阶数越高，滤波器在截止频率处的过渡越陡峭，但同时计算量也越大，且可能引入数值不稳定性。3 阶是一个平衡的选择。

第一行计算 Nyquist 频率 f_{nyq} = f_s / 2。第二行调用 butter 设计 Butterworth 低通滤波器，参数 fc / nyq 是归一化截止频率（0 到 1 之间，1 对应 Nyquist 频率），btype="low" 指定滤波器类型为低通。函数返回滤波器的分子系数 b 和分母系数 a，这两组系数定义了一个无限脉冲响应（IIR）滤波器，其传递函数为：

第三行调用 filtfilt(b, a, x) 进行零相位滤波。普通的 lfilter 函数按时间顺序逐点计算，会引入与频率相关的相位延迟。filtfilt 先正向滤波得到中间结果，再将中间结果时间反转后反向滤波，两次滤波的相位延迟正好抵消。代价是滤波器的幅度响应被应用了两次，等效阶数翻倍（3 阶滤波器经 filtfilt 后等效于 6 阶）。

z-score 标准化函数 znorm

def znorm(x: np.ndarray) -> np.ndarray:
    x = np.asarray(x)
    return (x - np.mean(x)) / (np.std(x) + 1e-12)

第一行 np.asarray(x) 确保输入被转换为 NumPy 数组。如果 x 本身已经是数组，该操作不会复制数据；如果 x 是 Python 列表或其他可迭代对象，则会创建新的数组。这增强了函数的鲁棒性，使其能接受多种输入类型。

第二行实现 z-score 公式。np.mean(x) 计算均值 \mu，np.std(x) 计算标准差 \sigma（默认使用 n 作为分母，即总体标准差）。分母加上 10^{-12} 是一个常见的数值稳定性技巧：当信号为常数时标准差为零，直接除以零会产生 inf 或 nan，加上极小值后结果虽然很大但仍是有限数，后续计算不会崩溃。

多信号绑图函数 plot_sigs

def plot_sigs(t: np.ndarray, sigs: Dict[str, np.ndarray], title: str = ""):
    n = len(sigs)
    fig, axes = plt.subplots(n, 1, sharex=True)
    if n == 1:
        axes = [axes]
    for ax, (name, y) in zip(axes, sigs.items()):
        ax.plot(t, y, lw=1.6)
        ax.set_ylabel(name)
    axes[-1].set_xlabel("Time (s)")
    if title:
        fig.suptitle(title)
    plt.show()

参数 sigs 的类型提示 Dict[str, np.ndarray] 表明它是一个字典，键为字符串（信号名称），值为 NumPy 数组（信号数据）。Python 3.7 以后字典保持插入顺序，因此信号在图中的排列顺序与字典中的顺序一致。

plt.subplots(n, 1, sharex=True) 创建 n 行 1 列的子图网格，sharex=True 使所有子图共享 x 轴的缩放和平移操作。当 n=1 时，axes 返回的是单个 Axes 对象而非数组，为了后续代码能统一处理，用 if n == 1: axes = [axes] 将其包装成列表。

zip(axes, sigs.items()) 将子图对象和信号键值对一一配对。sigs.items() 返回 (name, y) 元组的迭代器。循环中，ax.plot(t, y, lw=1.6) 在当前子图绑制时间序列，lw=1.6 设置线宽为 1.6 磅，比默认值略粗，更易于观察。ax.set_ylabel(name) 将信号名称设为 y 轴标签。

axes[-1].set_xlabel("Time (s)") 仅在最后一个子图设置 x 轴标签，因为所有子图共享 x 轴，只需标注一次。fig.suptitle(title) 在整个图像顶部添加总标题（如果提供了 title 参数）。最后 plt.show() 显示图像。

Cell 3：呼吸波形生成

本 Cell 实现呼吸波形（潮气量 V_t）的生成功能，这是整条仿真管线的输入端。呼吸波形描述了呼吸过程中肺部气体容量随时间的变化，后续的呼吸机械模型将以此为驱动信号，产生相应的神经输出。

呼吸频率与心脏相干性

呼吸频率通常以每分钟呼吸次数（Breaths Per Minute, BPM）来度量。正常成人的静息呼吸频率约为 12-20 BPM，但在本项目中，我们关注的是一种特定的呼吸调节技术——心脏相干性呼吸，其推荐频率为 6 BPM。

6 BPM 意味着每分钟呼吸 6 次，换算成周期是每次呼吸持续 10 秒，对应的频率为：

选择 0.1 Hz 作为目标呼吸频率有其生理学依据。人体的压力感受器反射系统（baroreflex）对 0.1 Hz 附近的输入具有共振特性，在此频率下呼吸能够最大限度地激发呼吸性窦性心律不齐（RSA），使心率随呼吸周期产生最大幅度的波动，从而增强心率变异性（HRV）。

波形类型

本项目支持三种呼吸波形类型，分别适用于不同的仿真场景。

正弦波（coherence / sine）

正弦波是最平滑的周期波形，其数学表达式为：

其中 f 是呼吸频率（Hz），t 是时间（秒）。2\pi f t 称为相位，单位是弧度，在一个周期内从 0 增加到 2\pi。正弦波的特点是连续可微，没有任何突变或尖角，最接近自然放松状态下的呼吸模式。在代码中，coherence 和 sine 两个选项是等价的，都生成正弦波。

平滑方波（square）

方波在吸气和呼气之间快速切换，吸气时间和呼气时间各占周期的一半。理想方波在切换瞬间发生跳变，但这种无限陡峭的边缘在生理上不现实，在数值计算中也可能导致问题。因此，实际实现时先生成正弦波，再通过 smooth_edge 函数（基于 \tanh）进行压缩：

当平滑参数 e 较小时，\tanh 将正弦波的正半周期压向 +1，负半周期压向 -1，形成接近方波的形状，但边缘处保持平滑过渡。

矩形波（rect）

矩形波与方波类似，但允许通过占空比（duty cycle）参数控制吸气时间和呼气时间的比例。生成过程首先计算每个时刻在当前周期内的相对位置（0 到 1 之间），若相对位置小于占空比则输出高电平，否则输出低电平，形成理想矩形波。然后同样应用 smooth_edge 函数平滑边缘。例如，占空比为 0.4 表示前 40% 的周期时间处于吸气状态，后 60% 处于呼气状态。

波形参数

bpm 参数指定呼吸频率，单位是次/分钟。函数内部将其转换为频率 f = \text{bpm} / 60（Hz），用于计算相位。

amp 参数控制波形的幅度，对应于潮气量的大小。幅度越大，表示呼吸越深。

base 参数是基线偏移量。在生理上，肺部即使在呼气末期也保留一定的气体容量，称为功能残气量（Functional Residual Capacity, FRC）。基线参数可以模拟这种偏移。

duty 参数是占空比，仅对矩形波有效，取值范围为 0 到 1，表示高电平时间占整个周期的比例。

edge_s 参数控制方波和矩形波的边缘平滑程度，值越大过渡越平缓。

最终输出

无论选择哪种波形类型，内部生成的归一化波形 y(t) 都在 [-1, 1] 范围内。最终输出的潮气量信号由以下公式计算：

这样，通过调节 baseline 和 amplitude 两个参数，可以将归一化波形映射到任意所需的数值范围，灵活模拟不同深度和基线的呼吸模式。

Cell 3：呼吸波形生成代码实现

本部分对 gen_breath 函数的代码实现进行逐行解析，说明各种波形的生成逻辑以及参数处理方式。

函数签名与参数处理

def gen_breath(
    t: np.ndarray,
    kind: str = "coherence",
    bpm: float = 6.0,
    amp: float = 1.0,
    base: float = 0.0,
    duty: float = 0.5,
    edge_s: float = 0.05,
) -> np.ndarray:

函数接受时间轴数组 t 作为必需参数，其余均为带默认值的可选参数。默认配置对应心脏相干性呼吸：6 BPM 的正弦波，幅度为 1，无基线偏移。

kind = kind.lower().strip()

这行代码将波形类型字符串转为小写并去除首尾空格，增强输入容错性。用户输入 "Coherence"、" SINE " 等都能被正确识别。

f = bpm / 60.0

将呼吸频率从 BPM 转换为 Hz。6 BPM 除以 60 得到 0.1 Hz，对应周期 T = 1/f = 10 秒。

phi = 2 * np.pi * f * t

计算各时刻的相位。相位 \phi(t) 是角度的度量，单位为弧度，其计算公式为：

当 t 从 0 增加到 T = 1/f 时，\phi 从 0 增加到 2\pi，完成一个完整的周期。由于 t 是数组，phi 也是同样长度的数组，存储了每个时刻对应的相位值。

正弦波生成

if kind in ("coherence", "sine"):
    y = np.sin(phi)

当波形类型为 coherence 或 sine 时，直接对相位数组应用正弦函数。np.sin 是向量化操作，一次性计算整个数组的正弦值。正弦函数的输出范围是 [-1, 1]，在 \phi = \pi/2 时达到最大值 1，在 \phi = 3\pi/2 时达到最小值 -1。

平滑方波生成

elif kind == "square":
    y = smooth_edge(np.sin(phi), e=edge_s)

方波的生成分两步：首先计算正弦波 np.sin(phi)，然后通过 smooth_edge 函数进行压缩。回顾 smooth_edge 的定义：

当 e 较小（如 0.05）时，\tanh 函数的过渡区域很窄。正弦波在正半周期的大部分时间内，\sin(\phi)/e 是一个较大的正数，\tanh 输出接近 1；在负半周期的大部分时间内，输出接近 -1。只有在正弦波过零点附近的短暂时间内，输出才在 -1 和 1 之间快速过渡，形成类似方波但边缘平滑的波形。

调整 edge_s 参数可以控制过渡的陡峭程度。edge_s = 0.02 产生非常接近理想方波的波形，edge_s = 0.1 则产生较为圆滑的过渡。

矩形波生成

elif kind == "rect":
    frac = (t * f) % 1.0

矩形波的生成逻辑与方波不同，它直接基于时间位置而非相位。t * f 计算从时间起点到当前时刻经过了多少个完整周期（包含小数部分）。例如，若 f = 0.1 Hz，则 t = 25 秒时 t \times f = 2.5，表示已经过了 2.5 个周期。取模 1 后得到 0.5，表示当前处于第三个周期的中间位置。

frac 数组的每个元素都在 [0, 1) 范围内，表示该时刻在当前周期内的相对位置：0 对应周期起点，0.5 对应周期中点，接近 1 对应周期末尾。

y_raw = (frac < duty).astype(float)

这行代码生成理想矩形波。布尔表达式 frac < duty 对数组中每个元素判断是否小于占空比，返回布尔数组。.astype(float) 将布尔值转换为浮点数（True 变为 1.0，False 变为 0.0）。

结果是：当时刻处于周期的前 duty 比例部分时，y_raw 为 1；否则为 0。例如 duty = 0.4 时，每个周期的前 40% 时间输出 1，后 60% 时间输出 0。

y = 2.0 * y_raw - 1.0

将 0/1 信号线性变换为 -1/+1 信号。原本的 0 变为 2 \times 0 - 1 = -1，原本的 1 变为 2 \times 1 - 1 = 1。这使得波形关于零对称，与正弦波和方波的范围一致。

y = smooth_edge(y, e=edge_s)

对理想矩形波应用边缘平滑。smooth_edge 函数将 -1 和 +1 之间的跳变转化为平滑过渡。由于 \tanh(\pm 1/e) 在 e 较小时已经非常接近 \pm 1，平台部分基本不受影响，只有跳变边缘被平滑。

错误处理与最终输出

else:
    raise ValueError(f"Unknown breathing pattern kind: {kind}")

如果用户传入了不支持的波形类型，函数抛出 ValueError 异常，错误信息中包含用户输入的值，便于调试。

return base + amp * y

最终输出通过线性变换将归一化波形映射到所需范围。amp * y 将波形幅度从 [-1, 1] 缩放到 [-\text{amp}, \text{amp}]，加上 base 后整体向上或向下平移。最终输出范围为 [\text{base} - \text{amp}, \text{base} + \text{amp}]。

演示代码

t_d = gen_taxis(60, 50)

生成 60 秒时长、50 Hz 采样率的时间轴。根据 gen_taxis 的计算，采样点数为 \lfloor 60 \times 50 \rfloor + 1 = 3001，时间间隔为 0.02 秒。

v1 = gen_breath(t_d, kind="coherence", bpm=6, amp=1.0)
v2 = gen_breath(t_d, kind="square", bpm=6, amp=1.0, edge_s=0.08)

分别生成正弦波和平滑方波两种呼吸波形，呼吸频率均为 6 BPM，幅度均为 1。方波的边缘平滑参数设为 0.08，比默认的 0.05 稍大，过渡更圆滑。

plot_sigs(t_d, {"Vt_coherence": v1, "Vt_square": v2}, title="Breathing patterns demo")

调用绑图函数，将两种波形绑制在上下两个子图中进行对比。信号名称 Vt_coherence 和 Vt_square 将显示为 y 轴标签。在 60 秒的时间范围内，6 BPM 的呼吸会产生 6 个完整周期，图中可以清晰观察到两种波形的周期相同但形状不同：正弦波平滑连续，方波在高低电平之间快速切换但边缘有轻微圆角。

Cell 4：呼吸机械模型（Vt → MN）

本 Cell 实现呼吸机械模型，将呼吸波形信号 V_t 转换为神经信号 MN。这是仿真管线的第二阶段，承上启下地连接呼吸输入端与心血管模型。呼吸机械模型模拟的是脑干呼吸中枢接收呼吸相关感觉输入后，产生调控心血管系统的神经输出的过程。

神经解剖学背景

呼吸活动的节律性由脑干的呼吸中枢产生和调控，主要涉及两个结构：延髓（medulla oblongata）和脑桥（pons）。延髓位于脑干的最下部，紧邻脊髓，其中包含多个与呼吸相关的神经元群。背侧呼吸组主要包含吸气神经元，在吸气相放电；腹侧呼吸组包含吸气和呼气神经元，控制呼吸肌肉的节律性收缩。延髓是产生基本呼吸节律的核心区域，即使与脑的其他部分隔离，延髓仍能产生原始的呼吸节律。

脑桥位于延髓上方，其中的呼吸调节中枢（包括脑桥呼吸组）负责调节呼吸的平滑过渡。脑桥接收来自延髓的信号，并向下反馈调节，使得吸气到呼气的转换更加协调，避免出现突兀的切换。可以将延髓理解为呼吸节律的发生器，脑桥则是节律的调节器和平滑器。

这些脑干神经核团发出的信号通过自主神经系统（交感神经和副交感神经）传递到心脏，影响窦房结的放电频率，从而导致心率随呼吸周期波动，这就是呼吸性窦性心律不齐（RSA）的神经机制。本项目用 \text{MN}_{\text{medulla}} 和 \text{MN}_{\text{pons}} 两个信号分别近似表示来自延髓和脑桥的运动神经元活动，作为下游心血管模型的驱动输入。

SimRes 数据类

为了统一管理各个模型的仿真输出，项目定义了 SimRes 数据类作为标准容器。该类使用 Python 的 @dataclass 装饰器定义，包含四个字段。t 字段存储时间轴数组，与输入时间轴相同。signals 字段是一个字典，键为信号名称（如 "hr"、"bp"、"mn_med"），值为对应的信号数组。states 字段是可选的 DataFrame，用于存储模型的内部状态变量，便于调试和详细分析。meta 字段也是可选的字典，记录模型名称、所用求解器等元信息。

这种封装设计的好处是，无论呼吸模型还是心血管模型，它们的 simulate 方法都返回相同类型的 SimRes 对象，上层代码可以用统一的方式访问仿真结果，无需关心具体模型的内部实现差异。

模型 A：低通滤波 + tanh 饱和

模型 A 是一个简单的信号处理模型，不涉及微分方程求解，计算速度快，适合快速原型验证。其信号处理流程分为三步。

第一步是低通滤波。原始的呼吸波形 V_t 可能包含高频噪声或不平滑的成分，通过低通滤波器可以滤除这些高频分量，保留呼吸节律的主要信息。默认截止频率设为 0.7 Hz，足以保留 6 BPM（0.1 Hz）的呼吸信号及其若干次谐波。

第二步是 z-score 标准化。滤波后的信号被标准化为零均值、单位方差的形式，消除幅度和偏置的影响，得到驱动信号：

第三步是通过 \tanh 非线性映射生成神经信号。延髓信号的计算公式为：

其中 b 是基线偏置（bias），代表神经元的基础放电率；g 是增益（gain），控制神经信号的波动幅度。\tanh 函数的输出范围在 (-1, 1) 之间，因此 \text{MN}_{\text{medulla}} 的范围被限制在 (b - g, b + g) 之内。这种饱和特性模拟了真实神经元的输入-输出关系：当输入很强时，神经元的放电率不会无限增加，而是趋于饱和。

脑桥信号的计算引入了时间延迟：

其中 \Delta t = 0.1 秒表示脑桥信号相对于延髓信号的延迟。系数 0.7 使脑桥信号的幅度略小于延髓信号。这种设计反映了生理上脑桥接收延髓输出并进行调节的时序关系。代码中使用 np.roll 函数实现信号的时间平移。

模型 B：二阶动力系统 + sigmoid

模型 B 是一个基于微分方程的动态模型，将呼吸机械系统建模为受外力驱动的二阶阻尼振荡器。与模型 A 的瞬时响应不同，模型 B 具有惯性，对输入变化的响应需要经过一定的动态过程。

系统的状态用位置 x 和速度 \dot{x} 描述，满足二阶线性常微分方程：

这个方程可以改写为：

其中 \omega_n 是系统的自然频率（natural frequency），单位是弧度/秒，决定了系统自由振荡的快慢；\zeta 是阻尼比（damping ratio），无量纲，决定了系统振荡衰减的快慢。

当阻尼比 \zeta < 1 时，系统处于欠阻尼状态，对阶跃输入的响应会出现振荡并逐渐衰减到稳态值；当 \zeta = 1 时为临界阻尼，系统以最快速度回到稳态且不振荡；当 \zeta > 1 时为过阻尼，系统缓慢地趋近稳态。模型 B 默认的阻尼比为 0.6，属于欠阻尼，系统响应会有轻微的超调但能较快稳定。

方程右边的 \omega_n^2 V_t 是驱动项，表示呼吸信号作为外力驱动系统运动。当系统达到稳态（\ddot{x} = \dot{x} = 0）时，x = V_t，即系统输出跟随输入。但在动态过程中，由于惯性和阻尼的存在，x 的变化相对于 V_t 有滞后和平滑效果。

求解这个二阶 ODE 需要将其转化为一阶系统。定义状态向量 \mathbf{y} = [x, \dot{x}]^T，则：

这就是传递给 solve_ivp 的右端函数 ode_rhs 所计算的内容。

神经信号通过 sigmoid 函数从系统状态生成。对位置状态进行 z-score 标准化后，延髓信号计算为：

其中 \sigma(z) 是 sigmoid 函数：

sigmoid 函数的输出范围是 (0, 1)，与 \tanh 的 (-1, 1) 不同。因此使用 sigmoid 时，\text{MN}_{\text{medulla}} 的范围是 (b, b + g)，始终为正，更符合神经元放电率非负的生理事实。

脑桥信号的生成更加复杂，融合了驱动信号及其导数：

其中 \text{drive}' 是驱动信号的时间导数，通过 np.gradient 计算。导数项的引入使脑桥信号对呼吸变化率有一定敏感性，在呼吸相位转换（吸气转呼气或呼气转吸气）时产生额外响应，模拟脑桥在呼吸相位切换时的调节作用。

两个模型的特性对比

模型 A 的优势在于简单高效，没有微分方程需要数值求解，计算速度快，适合快速测试管线的整体功能。其响应是即时的，除了低通滤波引入的少量延迟外，输出几乎立即跟随输入变化。

模型 B 的优势在于更接近真实的生理动态。二阶系统具有惯性，对输入变化的响应是渐进的，而非瞬时跳变。这种平滑的动态过渡更符合真实生物系统的特性。但模型 B 需要调用 ODE 求解器，计算成本较高。

两个模型使用的非线性函数也有所不同。\tanh 的输出关于零对称，范围 (-1, 1)；sigmoid 的输出始终为正，范围 (0, 1)。在神经生理学语境下，sigmoid 更适合表示放电率这种非负量，而 \tanh 更适合表示相对于基线的偏差。

Cell 4：呼吸机械模型代码实现

本部分对 Cell 4 中两个呼吸机械模型的代码进行逐行解析，重点说明数据类的设计、两种模型的实现差异以及 ODE 求解的技术细节。

SimRes 数据类定义

@dataclass
class SimRes:
    t: np.ndarray
    signals: Dict[str, np.ndarray]
    states: Optional[pd.DataFrame] = None
    meta: Optional[Dict[str, Any]] = None

@dataclass 装饰器告诉 Python 自动为这个类生成 __init__ 方法。根据字段定义，自动生成的构造函数等价于：

def __init__(self, t, signals, states=None, meta=None):
    self.t = t
    self.signals = signals
    self.states = states
    self.meta = meta

前两个字段 t 和 signals 没有默认值，是必需参数；后两个字段有默认值 None，是可选参数。Optional[pd.DataFrame] 表示该字段可以是 pd.DataFrame 类型或 None。

使用 dataclass 而非普通类的好处是代码更简洁，且自动获得合理的字符串表示（__repr__），便于调试时查看对象内容。

RespModelA 类实现

构造函数

def __init__(self, fs: float, fc: float = 0.7, g: float = 6.0, b: float = 8.0):
    self.fs = fs
    self.fc = fc
    self.g = g
    self.b = b

构造函数接收四个参数并存储为实例属性。fs 是采样频率，必须与输入信号的采样频率一致，否则滤波和时间延迟计算会出错。fc 是低通滤波的截止频率，默认 0.7 Hz。g 和 b 分别是增益和偏置，默认值 6.0 和 8.0 使得输出范围约为 (2, 14)。

simulate 方法

def simulate(self, t: np.ndarray, vt: np.ndarray) -> SimRes:
    vt_f = lpf(vt, fs=self.fs, fc=self.fc, n_ord=3)

第一步调用低通滤波函数 lpf，传入信号 vt、采样频率 self.fs、截止频率 self.fc 和滤波器阶数 3。返回的 vt_f 是滤波后的信号，高于 0.7 Hz 的频率成分被衰减。

    drv = znorm(vt_f)

第二步对滤波后信号进行 z-score 标准化。标准化后 drv 的均值为 0，标准差为 1，大部分值落在 [-3, 3] 范围内（假设近似正态分布）。

    mn_med = self.b + self.g * np.tanh(drv)

第三步生成延髓信号。np.tanh(drv) 对整个数组逐元素计算双曲正切，输出范围 (-1, 1)。乘以增益 self.g 后范围变为 (-g, g)，加上偏置 self.b 后最终范围为 (b-g, b+g)。以默认参数为例，b=8, g=6，输出范围约为 (2, 14)。

    mn_pons = self.b + 0.7 * self.g * np.tanh(np.roll(drv, int(0.1 * self.fs)))

第四步生成脑桥信号。np.roll(drv, int(0.1 * self.fs)) 将数组循环平移指定的位置数。int(0.1 * self.fs) 计算 0.1 秒对应的采样点数：若 fs=50，则为 5 个点。正数表示向后平移，即信号延迟。平移后数组末尾的元素会被移到开头，形成循环，但由于仿真时间通常远长于 0.1 秒，这种边界效应的影响可以忽略。

系数 0.7 使脑桥信号的幅度为延髓信号的 70%，反映脑桥作为调节器而非主发生器的生理角色。

    return SimRes(
        t=t,
        signals={"vt": vt, "mn_med": mn_med, "mn_pons": mn_pons},
        meta={"model": "RespModelA"},
    )

最后构造并返回 SimRes 对象。signals 字典包含原始输入 vt 和两个输出信号。meta 记录模型名称，便于后续识别结果来源。

RespModelB 类实现

构造函数

def __init__(self, wn: float = 2*np.pi*0.25, zeta: float = 0.6, g: float = 10.0, b: float = 7.0):
    self.wn = float(wn)
    self.zeta = float(zeta)
    self.g = float(g)
    self.b = float(b)

默认自然频率 \omega_n = 2\pi \times 0.25 \approx 1.57 rad/s，对应的固有频率为 0.25 Hz，周期 4 秒。这意味着系统对周期约 4 秒的输入会有较强的响应。默认阻尼比 \zeta = 0.6 处于欠阻尼区域，系统响应会有轻微振荡但很快稳定。

float() 转换确保即使传入整数，存储的也是浮点数，避免后续计算中的整数除法问题。

simulate 方法中的 ODE 定义

def simulate(self, t: np.ndarray, vt: np.ndarray, solver: str = "RK45") -> SimRes:
    def ode_rhs(tt, y):
        x, xdot = y
        u = np.interp(tt, t, vt)
        xdd = (self.wn**2) * (u - x) - 2*self.zeta*self.wn*xdot
        return [xdot, xdd]

ode_rhs 是定义在 simulate 方法内部的嵌套函数，它可以访问外层函数的变量（如 t、vt、self），这种特性称为闭包。

参数 tt 是 ODE 求解器当前求解的时刻（标量），y 是当前状态向量。x, xdot = y 将状态向量解包为位置和速度两个分量。

np.interp(tt, t, vt) 在时刻 tt 对输入信号进行线性插值。ODE 求解器使用自适应步长，求解时刻 tt 不一定落在预定义的时间网格 t 上，因此需要插值获取该时刻的输入值。np.interp 的三个参数分别是：插值点、已知数据的 x 坐标、已知数据的 y 坐标。

xdd 是加速度 \ddot{x}，根据二阶系统方程计算：

函数返回 [\dot{x}, \ddot{x}]，即状态向量的时间导数。

ODE 求解

    y0 = [vt[0], 0.0]
    sol = solve_ivp(ode_rhs, (t[0], t[-1]), y0, t_eval=t, method=solver, rtol=1e-6, atol=1e-8)

初始条件设置位置 x(0) = V_t(0)（与输入初值相同），速度 \dot{x}(0) = 0（静止状态开始）。

solve_ivp 的参数说明：第一个参数是右端函数；第二个参数 (t[0], t[-1]) 是求解的时间区间，从时间轴的起点到终点；y0 是初始状态；t_eval=t 指定在哪些时刻输出解，这里设为输入的时间网格，确保输出与输入对齐；method=solver 指定数值方法；rtol 和 atol 分别是相对误差和绝对误差容限，控制求解精度。

返回的 sol 对象包含求解结果。sol.y 是形状为 (2, len(t)) 的数组，sol.y[0] 是位置 x(t) 的时间序列，sol.y[1] 是速度 \dot{x}(t) 的时间序列。

神经信号生成

    x = sol.y[0]
    drv = znorm(x)
    mn_med = self.b + self.g * (1 / (1 + np.exp(-drv)))

提取位置状态并标准化后，通过 sigmoid 函数生成延髓信号。1 / (1 + np.exp(-drv)) 是 sigmoid 函数的直接实现，等价于 scipy.special.expit(drv)。

    mn_pons = self.b + 0.8 * self.g * (1 / (1 + np.exp(-(np.gradient(drv) * 0.7 + drv))))

脑桥信号的计算更复杂。np.gradient(drv) 计算 drv 的数值导数，使用中心差分公式（边界处使用单侧差分）。导数项乘以 0.7 的权重后与原信号相加，再通过 sigmoid 变换。

这种设计使脑桥信号在 drv 变化剧烈时（导数大）产生额外响应，而在 drv 平稳时主要跟随 drv 本身。系数 0.8 使脑桥信号幅度略小于延髓信号。

    st = pd.DataFrame({"x": x, "xdot": sol.y[1]}, index=t)

将内部状态存储为 DataFrame，以时间 t 作为索引。这便于后续检查系统的动态行为，例如绑制相空间轨迹或分析瞬态响应。

演示代码

t_d = gen_taxis(60, 50)
vt = gen_breath(t_d, kind="coherence", bpm=6, amp=1.0)
resA = RespModelA(fs=50).simulate(t_d, vt)
resB = RespModelB().simulate(t_d, vt, solver="RK45")

创建 60 秒、50 Hz 的时间轴和 6 BPM 正弦呼吸波形后，分别用两个模型进行仿真。模型 A 需要指定采样频率 fs=50，模型 B 使用默认参数。两个模型的输出都是 SimRes 对象，通过 .signals["mn_med"] 访问延髓信号。

plot_sigs(t_d, {"Vt": vt, "MN_med_A": resA.signals["mn_med"], "MN_med_B": resB.signals["mn_med"]},
                title="Resp models output demo")

绑图展示输入波形和两个模型的延髓输出。在图中可以观察到，模型 A 的输出几乎与输入同步变化，仅有低通滤波引入的轻微平滑；模型 B 的输出则表现出明显的动态响应特性，相对于输入有滞后，且波形更加平滑，这是二阶系统惯性的体现。

Cell 5：心血管模型（MN → HR/BP）

本 Cell 实现心血管模型，将来自呼吸机械模型的神经信号（MN）转换为心率（HR）和血压（BP）。这是整条仿真管线的核心部分，模拟了自主神经系统对心血管系统的调控作用。

心血管调节的生理机制

心脏的跳动节律由位于右心房的窦房结（Sinoatrial Node, SA Node）自发产生，但这个固有节律受到自主神经系统的持续调控。自主神经系统的两个分支以相反的方式影响心率：副交感神经（主要通过迷走神经）释放乙酰胆碱，作用于窦房结使心率降低；交感神经释放去甲肾上腺素，使心率升高。在静息状态下，副交感神经占主导，因此静息心率（约 60-70 BPM）低于窦房结的固有频率（约 100 BPM）。

血管张力（Vascular Tone）描述血管壁平滑肌的收缩状态。交感神经兴奋使血管收缩，血管内径减小，外周阻力增加；交感神经抑制或副交感神经兴奋则使血管舒张，外周阻力降低。血管张力是决定血压的关键因素之一。

血压（Blood Pressure）由心输出量（Cardiac Output）和外周阻力（Peripheral Resistance）共同决定，可以用类似于欧姆定律的关系来理解：

心输出量等于心率乘以每搏输出量。在本项目的简化模型中，假设每搏输出量恒定，因此心输出量与心率成正比。外周阻力则与血管张力正相关。这意味着心率升高或血管收缩都会导致血压升高。

面向对象的模型架构

本项目采用面向对象设计，定义了心血管模型的抽象基类 CardioBase，所有具体的心血管模型（ToyCardio、Ott97、SH95）都继承自这个基类。基类定义了统一的仿真接口：

任何心血管模型的 simulate 方法都接受时间轴和两个神经信号作为输入，返回一个 SimRes 对象，其中 signals 字典至少包含 'hr' 和 'bp' 两个键。

这种设计的好处是实现了多态性。上层代码（如 GUI、HRV 分析模块）只需要与基类接口交互，无需关心具体使用的是哪个模型。当需要切换模型或添加新模型时，只要新模型实现了相同的接口，上层代码无需任何修改。

ToyCardio 简化模型

由于项目要求的 Ott97 和 SH95 模型的原始文献难以获取完整的方程和参数，本项目实现了一个工程上可用的简化模型 ToyCardio，用于验证整条仿真管线的功能。虽然简化，但该模型能产生合理的心率变异性和呼吸-心率耦合现象。

状态变量与一阶滞后动力学

模型包含三个状态变量：心率 HR、血管张力 Tone、血压 BP。每个状态变量都遵循一阶滞后动力学，即状态变量以指数形式趋近于其目标值，趋近的速度由时间常数 \tau 决定。

心率的动力学方程为：

这个方程可以改写为：

方程的含义是：心率的变化率与当前心率和目标心率之间的差值成正比。当 HR 小于目标值时，导数为正，HR 增加；当 HR 大于目标值时，导数为负，HR 减少。时间常数 \tau_{\text{hr}} 决定了这个趋近过程的快慢：\tau 越大，趋近越慢，系统惯性越大。

对于阶跃输入（目标值突然变化），一阶系统的解是指数衰减形式。若目标值在 t=0 时从 \text{HR}_0 跳变到 \text{HR}_{\text{target}}，则响应为：

经过一个时间常数 \tau_{\text{hr}} 后，HR 与目标值的差距缩小到初始差距的 e^{-1} \approx 37\%；经过 3\tau_{\text{hr}} 后缩小到约 5%，可以认为基本达到稳态。

血管张力和血压的动力学方程形式相同：

默认的时间常数设置为 \tau_{\text{hr}} = 3 秒、\tau_{\text{tone}} = 5 秒、\tau_{\text{bp}} = 4 秒。这些值反映了不同生理过程的响应速度：心率变化最快，血管张力变化最慢。

目标值的计算

各状态变量的目标值由神经信号通过非线性映射确定。神经信号首先经过 z-score 标准化，消除来自不同呼吸模型的尺度差异。

心率目标值的计算为：

其中 \text{HR}_0 是静息心率（默认 70 BPM），m_1 和 m_2 分别是标准化后的延髓和脑桥信号，k_{\text{hr}} 是耦合增益。延髓信号的权重 0.7 大于脑桥信号的权重 0.3，反映延髓在心率调节中的主导作用。\tanh 函数将输出限制在 \pm 8 BPM 范围内，因此心率目标值在 (62, 78) BPM 之间波动。

血管张力目标值的计算为：

这里脑桥信号的权重 0.8 大于延髓信号的权重 0.2，反映脑桥在血管调节中的主导作用。输出范围为 \pm 0.5，以零为中心波动。

血压目标值的计算结合了心率和血管张力：

第一项 \text{BP}_0 是静息血压（默认 100 mmHg）。第二项反映心率对血压的影响：心率偏离静息值时，血压相应变化。第三项反映血管张力对血压的影响：血管收缩（Tone 为正）使血压升高，血管舒张（Tone 为负）使血压降低。

占位类 Ott97 与 SH95

项目要求中提到的 Ott97 模型（Ottesen 1997）和 SH95 模型（Seidel & Herzel 1995）是文献中的经典心血管模型。Ott97 模型基于压力感受器反射的数学描述，SH95 模型关注心率变异性的非线性动力学。

由于原始文献的完整方程和参数难以获取，这两个类目前仅作为占位实现。它们继承自 CardioBase 基类，定义了与 ToyCardio 相同的接口，但 simulate 方法会抛出 NotImplementedError 异常，提示需要补充具体实现。

这种设计保留了扩展性：一旦获得模型的完整细节，只需在 simulate 方法中填入 ODE 右端函数和参数，整个仿真管线即可无缝切换到新模型。

max_step 参数的兼容性处理

在使用 solve_ivp 求解 ODE 时，max_step 参数用于限制积分步长的最大值。某些情况下需要限制步长以确保不错过信号的快速变化。然而，不同版本的 SciPy 对 max_step=None 的处理存在差异：某些版本会将 None 与数值进行比较（如 None <= 0），导致 TypeError。

解决方案是在调用 solve_ivp 前进行条件判断：如果 max_step 为 None，则不将该参数传入函数，让求解器使用默认的自适应步长策略；只有当 max_step 有具体数值时，才将其加入参数字典。这种处理方式确保代码在不同版本的 SciPy 下都能正常运行。

Cell 5：心血管模型代码实现

本部分对 Cell 5 中心血管模型的代码进行逐行解析，重点说明面向对象设计的实现方式、ODE 系统的构建以及参数的生理学含义。

辅助函数与数据类的重复定义

def znorm(x: np.ndarray) -> np.ndarray:
    x = np.asarray(x)
    return (x - np.mean(x)) / (np.std(x) + 1e-12)

@dataclass
class SimRes:
    t: np.ndarray
    signals: Dict[str, np.ndarray]
    states: Optional[pd.DataFrame] = None
    meta: Optional[Dict[str, Any]] = None

这两个定义与 Cell 2 和 Cell 4 中的完全相同。重复定义的目的是使 Cell 5 可以独立运行，便于调试和测试。在完整的 notebook 执行流程中，后定义的版本会覆盖先前的版本，但由于实现完全相同，不会产生任何功能差异。

CardioBase 基类定义

class CardioBase:
    name = "CardioBase"

    def simulate(
        self,
        t: np.ndarray,
        mn_med: np.ndarray,
        mn_pons: np.ndarray,
        solver: str = "RK45",
        rtol: float = 1e-6,
        atol: float = 1e-8,
        max_step: Optional[float] = None,
    ) -> SimRes:
        raise NotImplementedError

CardioBase 是所有心血管模型的抽象基类。类属性 name 存储模型名称，子类应覆盖此属性以标识自己。

simulate 方法定义了统一的仿真接口。参数包括时间轴 t、两个神经信号 mn_med 和 mn_pons，以及 ODE 求解器的配置参数。返回类型为 SimRes，约定 signals 字典必须包含 'hr' 和 'bp' 两个键。

基类的 simulate 方法直接抛出 NotImplementedError，这是 Python 中实现抽象方法的常用模式。任何直接调用基类 simulate 的尝试都会失败，强制要求使用者必须使用具体的子类实现。

ToyCardio 构造函数

class ToyCardio(CardioBase):
    name = "ToyCardio"

    def __init__(self, hr0: float = 70.0, bp0: float = 100.0):
        self.hr0 = float(hr0)
        self.bp0 = float(bp0)
        self.tau_hr = 3.0
        self.tau_tone = 5.0
        self.tau_bp = 4.0
        self.k_hr = 0.9
        self.k_tone = 0.6
        self.k_bp_hr = 0.25
        self.k_bp_tone = 12.0

ToyCardio 继承自 CardioBase，通过 class ToyCardio(CardioBase) 语法声明继承关系。子类自动获得父类的所有属性和方法，同时可以覆盖或扩展它们。

构造函数接受两个可选参数：静息心率 hr0 和静息血压 bp0。float() 转换确保即使传入整数也存储为浮点数。

时间常数 tau_hr、tau_tone、tau_bp 决定各状态变量的响应速度。以 tau_hr = 3.0 秒为例，当心率目标值发生阶跃变化时，实际心率需要约 3 秒才能完成 63% 的调整，约 9 秒（3 个时间常数）才能基本达到稳态。三个时间常数的相对大小（心率最快、张力最慢）反映了不同生理过程的特征时间尺度。

耦合增益 k_hr、k_tone、k_bp_hr、k_bp_tone 决定各变量之间的影响强度。这些参数经过调整，使模型在 6 BPM 呼吸输入下产生约 8 BPM 的心率波动幅度，与真实心脏相干性呼吸的效果相近。

ToyCardio simulate 方法：输入处理

def simulate(
    self,
    t: np.ndarray,
    mn_med: np.ndarray,
    mn_pons: np.ndarray,
    solver: str = "RK45",
    rtol: float = 1e-6,
    atol: float = 1e-8,
    max_step: Optional[float] = None,
) -> SimRes:
    t = np.asarray(t)
    mn_med = np.asarray(mn_med)
    mn_pons = np.asarray(mn_pons)

    if t.ndim != 1:
        raise ValueError("t must be 1D")
    if len(mn_med) != len(t) or len(mn_pons) != len(t):
        raise ValueError("mn_medulla/mn_pons must have same length as t")

方法开始时将所有输入转换为 NumPy 数组，增强对不同输入类型（如 Python 列表）的兼容性。随后进行输入校验：时间轴必须是一维数组，两个神经信号的长度必须与时间轴匹配。这种防御性编程可以在错误输入时立即给出清晰的错误信息，而非在后续计算中产生难以追踪的异常。

    m1 = znorm(mn_med)
    m2 = znorm(mn_pons)

对神经信号进行 z-score 标准化。不同的呼吸模型可能产生不同范围的神经信号（例如模型 A 输出范围约 2-14，模型 B 可能是 7-17），标准化后两者都变为零均值、单位方差，使得心血管模型的参数具有通用性，无需针对不同呼吸模型分别调整。

ToyCardio simulate 方法：ODE 右端函数

    def ode_rhs(tt, y):
        hr, tone, bp = y

        u1 = np.interp(tt, t, m1)
        u2 = np.interp(tt, t, m2)

ODE 右端函数 ode_rhs 定义为嵌套函数，可以访问外层的 t、m1、m2 等变量。参数 tt 是当前求解时刻（标量），y 是当前状态向量。

np.interp(tt, t, m1) 在时刻 tt 对信号 m1 进行线性插值。由于 ODE 求解器使用自适应步长，求解时刻 tt 不一定落在预定义的时间网格上，因此需要插值获取该时刻的输入值。

        hr_tgt = self.hr0 + 8.0 * np.tanh(self.k_hr * (0.7 * u1 + 0.3 * u2))

计算心率目标值。表达式的结构是：基础值 + 调制范围 × 非线性函数（增益 × 加权输入）。

加权组合 0.7 * u1 + 0.3 * u2 使延髓信号 u1 的影响大于脑桥信号 u2。乘以耦合增益 self.k_hr = 0.9 后，输入范围约为 \pm 0.9（假设标准化信号主要在 \pm 1 内）。\tanh(0.9) \approx 0.72，乘以调制范围 8.0 后，心率目标值在静息值 70 BPM 上下约 \pm 5.8 BPM 范围内波动。\tanh 函数的饱和特性确保即使输入异常大，输出也不会超出 \pm 8 BPM。

        tone_tgt = 0.5 * np.tanh(self.k_tone * (0.2 * u1 + 0.8 * u2))

计算血管张力目标值。这里脑桥信号的权重 0.8 大于延髓信号的权重 0.2，反映脑桥在外周血管调节中的主导作用。基础值为 0，表示张力围绕中性状态波动；调制范围 0.5 使输出在 \pm 0.5 范围内。

        dhr = (hr_tgt - hr) / self.tau_hr
        dtone = (tone_tgt - tone) / self.tau_tone

心率和张力的一阶滞后动力学。以心率为例，当目标值 hr_tgt 高于当前值 hr 时，差值为正，dhr 为正，心率增加；反之则减少。除以时间常数 tau_hr 控制变化速率：时间常数越大，变化越慢。

        bp_tgt = self.bp0 + self.k_bp_hr * (hr - self.hr0) + self.k_bp_tone * tone

计算血压目标值。第一项 self.bp0 是静息血压。第二项 self.k_bp_hr * (hr - self.hr0) 反映心率对血压的影响：当心率高于静息值时，该项为正，使血压目标值升高。第三项 self.k_bp_tone * tone 反映血管张力对血压的影响：张力为正（血管收缩）时血压升高，张力为负（血管舒张）时血压降低。

以默认参数为例，若心率升高 8 BPM（从 70 到 78），第二项贡献 0.25 \times 8 = 2 mmHg；若同时张力为 0.3，第三项贡献 12.0 \times 0.3 = 3.6 mmHg。血压目标值相对于静息值升高约 5.6 mmHg。

        dbp = (bp_tgt - bp) / self.tau_bp
        return [dhr, dtone, dbp]

血压同样采用一阶滞后动力学。函数返回三个状态变量导数组成的列表，供 ODE 求解器使用。

ToyCardio simulate 方法：ODE 求解

    y0 = [self.hr0, 0.0, self.bp0]

设置初始条件：心率和血压从静息值开始，血管张力从零（中性状态）开始。

    opts = dict(
        t_eval=t,
        method=solver,
        rtol=rtol,
        atol=atol,
    )
    if max_step is not None:
        opts["max_step"] = float(max_step)

    sol = solve_ivp(ode_rhs, (float(t[0]), float(t[-1])), y0, **opts)

这段代码展示了处理可选参数的技巧。首先创建包含必需参数的字典 opts，然后仅当 max_step 不为 None 时才将其加入字典。最后使用 **opts 语法将字典解包为关键字参数传递给 solve_ivp。

这种处理方式解决了 SciPy 版本兼容性问题：某些版本的 solve_ivp 在收到 max_step=None 时会尝试执行 None <= 0 的比较，导致 TypeError。通过条件判断避免传入 None 值，确保代码在各版本下都能正常运行。

float(t[0]) 和 float(t[-1]) 确保时间边界是 Python 原生浮点数而非 NumPy 标量，避免某些版本的 SciPy 对类型的严格检查。

    if not sol.success:
        raise RuntimeError(f"solve_ivp failed: {sol.message}")

检查求解是否成功。solve_ivp 返回的对象包含 success 布尔属性和 message 字符串属性。如果求解因误差超限、步长过小等原因失败，success 为 False，此时抛出错误并附上失败原因。

    hr, tone, bp = sol.y
    st = pd.DataFrame({"hr": hr, "tone": tone, "bp": bp}, index=t)

sol.y 是形状为 (3, len(t)) 的数组，每行对应一个状态变量的时间序列。通过序列解包提取三个变量。随后将状态变量组织为 DataFrame，以时间为索引，便于后续分析和可视化。

    return SimRes(
        t=t,
        signals={"hr": hr, "bp": bp},
        states=st,
        meta={
            "model": self.name,
            "solver": solver,
            "rtol": rtol,
            "atol": atol,
            "max_step": max_step,
        },
    )

构造并返回 SimRes 对象。signals 字典仅包含 hr 和 bp，符合基类接口约定；tone 作为内部状态存储在 states DataFrame 中，不直接暴露给外部。meta 字典记录了模型名称和所有求解器参数，便于追溯结果的生成条件。

占位类 Ott97 与 SH95

class Ott97(CardioBase):
    name = "Ott97 (TBD)"

    def __init__(self, params: Optional[Dict[str, Any]] = None):
        self.params = params or {}

    def simulate(self, t, mn_med, mn_pons, solver="RK45", rtol=1e-6, atol=1e-8, max_step=None) -> SimRes:
        raise NotImplementedError(
            "Ott97Model 需要你们提供 rhs() 方程与参数。\n"
            "也请采用与 ToyCardioModel 相同的写法：max_step 为 None 时不要传给 solve_ivp。"
        )

Ott97 和 SH95 两个类的结构相同，都是占位实现。它们继承自 CardioBase，定义了与 ToyCardio 相同的接口，但 simulate 方法直接抛出 NotImplementedError，错误信息中提示了实现要求。

构造函数接受一个可选的参数字典 params。表达式 params or {} 利用了 Python 的短路求值：如果 params 为 None（假值），则返回空字典 {}；否则返回 params 本身。这确保 self.params 始终是一个字典，后续代码可以安全地对其进行操作。

名称中的 (TBD) 表示 To Be Done（待完成），在 GUI 中显示时可以提醒用户该模型尚未实现。当用户在界面中选择这些模型并点击运行时，会看到 NotImplementedError 异常信息，而非程序崩溃。

Cell 6：HRV 分析（时域指标 + 相干性）

本 Cell 实现心率变异性（Heart Rate Variability, HRV）分析功能，这是整条仿真管线的最终输出环节。HRV 描述了心跳间隔的波动特征，反映自主神经系统对心脏的动态调控能力。通过分析 HRV，可以量化评估心脏相干性呼吸技术的效果。

从连续心率到 RR 间期序列

心血管模型输出的是连续的心率信号 HR(t)，单位是每分钟心跳次数（BPM）。然而，标准的 HRV 分析是基于 RR 间期进行的。RR 间期是心电图上相邻两个 R 波（心室去极化的峰值）之间的时间间隔，也称为 NN 间期（Normal-to-Normal interval，排除异常心搏后的正常间期）。

从瞬时心率到瞬时 RR 间期的转换关系是简单的倒数关系：

例如，心率为 60 BPM 时，每分钟跳 60 次，相邻心搏间隔为 1 秒；心率为 120 BPM 时，间隔为 0.5 秒。

然而，RR 间期序列与连续心率信号有本质区别。连续心率信号是等间隔采样的时间序列，例如每 0.02 秒（50 Hz）一个采样点；而 RR 间期序列是事件型数据，每次心跳产生一个数据点，采样时刻不等间隔。在静息心率 70 BPM 时，平均每 0.86 秒产生一个 RR 值；心率波动时，相邻 RR 值的时间间隔也在变化。

相位累积方法

本项目采用相位累积方法从连续心率信号模拟生成 RR 事件序列。其基本思想是：将心跳过程视为一个相位不断累积的振荡器，每当累积相位增加 1（完成一个完整周期），就发生一次心跳事件。

具体实现中，首先计算每个采样时刻的瞬时 RR 间期：

然后计算相位增量。在时间步长 \Delta t 内，相位增量等于该时间段占当前 RR 间期的比例：

对相位增量进行累积求和得到累积相位 \phi(t)。当 \phi(t) 的整数部分发生跳变（从 n 变为 n+1）时，表示完成了一个心跳周期，该时刻即为心跳发生的时刻。记录这些时刻及其对应的瞬时 RR 值，就得到了 RR 事件序列。

最后，对生成的 RR 序列进行生理合理性过滤，剔除 RR 小于 0.25 秒（对应心率大于 240 BPM）或大于 2.5 秒（对应心率小于 24 BPM）的异常值。

时域 HRV 指标

时域分析直接基于 RR 间期序列计算统计指标，不涉及频率变换，计算简单直观。

平均 RR 间期（Mean RR）

平均 RR 间期是所有 RR 值的算术平均，单位通常为毫秒（ms）：

平均 RR 间期的倒数乘以 60000 即为平均心率（BPM）。该指标反映心脏跳动的基本节奏。

SDNN

SDNN 是 RR 间期的标准差（Standard Deviation of NN intervals），反映 HRV 的总体水平：

SDNN 综合反映了所有影响心率变异的因素，包括长期的昼夜节律变化和短期的呼吸调制。SDNN 值越大，表示心率变异性越高，通常与更好的心血管健康状态相关。典型的健康成人 SDNN 约为 100-150 ms。

RMSSD

RMSSD 是相邻 RR 间期差值的均方根（Root Mean Square of Successive Differences）：

RMSSD 主要反映短期的逐拍变异性，与副交感神经（迷走神经）活动密切相关。由于副交感神经对心率的调节作用非常迅速（可在一个心跳周期内生效），而交感神经的作用相对缓慢，因此相邻心搏之间的快速变化主要反映副交感活动。RMSSD 是评估副交感神经张力的敏感指标。

pNN50

pNN50 是相邻 RR 间期差值绝对值大于 50 ms 的比例（Percentage of NN50）：

pNN50 与 RMSSD 类似，也主要反映副交感神经活动。50 ms 的阈值是根据经验确定的，能够有效区分正常的逐拍变异与较大的心率波动。pNN50 的优点是不受极端值的影响（因为只考虑是否超过阈值，不考虑超过多少），但缺点是当样本量较小或变异性较低时，该指标可能为零或接近零，分辨力下降。

频域分析与相干性

频域分析通过功率谱密度（Power Spectral Density, PSD）揭示 HRV 信号中不同频率成分的能量分布。

RR 序列的重采样

频域分析方法（如 FFT、Welch 方法）要求输入信号是等间隔采样的。然而，RR 事件序列天然是不等间隔的——心跳间隔本身就在变化。因此，在进行频域分析之前，需要将不等间隔的 RR 序列重采样为等间隔的时间序列。

重采样的方法是线性插值：在指定的等间隔时间网格上，对原始 RR 序列进行线性内插，得到等间隔采样的 RR 信号。常用的重采样频率为 4 Hz，即每 0.25 秒一个采样点。4 Hz 的选择是基于 HRV 感兴趣的频率范围（通常 0-0.5 Hz）和 Nyquist 定理的考虑：4 Hz 采样可以正确表示最高 2 Hz 的频率成分，远超 HRV 的有效带宽。

Welch 方法估计功率谱

Welch 方法是一种改进的周期图方法，通过分段平均减少功率谱估计的方差。具体步骤是：将信号分成多个重叠的段，对每段加窗后计算周期图（FFT 幅度平方），然后对所有段的周期图取平均。

与直接对整段信号做 FFT 相比，Welch 方法得到的功率谱更平滑、更稳定，代价是频率分辨率略有下降。在 HRV 分析中，信号长度有限且存在噪声，Welch 方法的平滑特性是有利的。

相干性分数的定义

相干性分数（Coherence Score）量化了 HRV 信号在目标呼吸频率附近的功率集中程度。其定义为目标频段功率与总功率的比值：

其中 P_{\text{total}} 是功率谱在全频段的积分（总功率），P_{\text{target}} 是目标频段内的功率积分。目标频段通常定义为目标频率 \pm 0.02 Hz，例如对于 6 BPM（0.1 Hz）的呼吸，目标频段为 0.08-0.12 Hz。

相干性分数的取值范围是 0 到 1。当相干性接近 1 时，说明 HRV 的能量高度集中在目标呼吸频率附近，心率波动与呼吸节律高度同步，这正是心脏相干性呼吸所追求的状态。相干性较低则表示心率波动较为杂乱，能量分散在多个频率上，呼吸与心率的耦合不够紧密。

功率的计算使用梯形积分法（np.trapz），对离散的功率谱密度值在频率轴上进行数值积分。

主导频率

除相干性外，还可以提取功率谱的主导频率（Dominant Frequency），即功率谱密度最大值对应的频率。理想情况下，在心脏相干性呼吸期间，主导频率应接近呼吸频率（如 0.1 Hz）。如果主导频率偏离目标频率较多，可能说明呼吸节律不稳定或心血管系统对呼吸调制的响应不佳。

全流程演示

Cell 6 末尾的代码将前面所有模块串联起来，形成完整的仿真管线：首先生成时间轴和呼吸波形，然后通过呼吸机械模型转换为神经信号，再由心血管模型产生心率和血压，最后从心率信号提取 RR 序列并计算 HRV 指标。这段演示代码验证了整条管线可以正常运行，各模块之间的接口正确对接。

运行演示代码需要先执行 Cell 1-5，确保所有依赖的函数和类已经定义。如果直接运行 Cell 6 而未执行前面的 Cell，会因为时间轴 t 和心血管仿真结果 cres 未定义而报错。

Cell 6：HRV 分析代码实现

本部分对 Cell 6 中 HRV 分析函数的代码进行逐行解析，重点说明相位累积方法的实现原理、时域指标的计算细节以及频域相干性分析的技术要点。

hr2rr 函数：心率到 RR 序列的转换

def hr2rr(t: np.ndarray, hr: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    hr_c = np.clip(hr, 30, 220)

函数接受时间轴 t 和心率信号 hr 作为输入，返回两个数组的元组：心跳发生时刻和对应的 RR 间期。

np.clip(hr, 30, 220) 将心率值限制在 30-220 BPM 范围内。低于 30 BPM 的值被截断为 30，高于 220 BPM 的值被截断为 220。这是一种防御性处理：极端的心率值（可能由仿真异常产生）会导致后续计算出现问题，例如 HR 接近 0 时 RR 趋向无穷大。30-220 BPM 涵盖了从严重心动过缓到极限运动的生理范围。

    irr = 60.0 / hr_c

计算瞬时 RR 间期数组。根据心率与 RR 间期的倒数关系，60（秒/分钟）除以心率（次/分钟）得到 RR 间期（秒/次）。结果 irr 与输入 hr 等长，每个元素表示对应时刻的瞬时 RR 值。

    dt = np.diff(t)
    dt = np.concatenate([dt, dt[-1:]])

计算时间步长数组。np.diff(t) 计算相邻时间点的差值，返回长度为 len(t) - 1 的数组。为了使 dt 与 t 等长，用最后一个步长值 dt[-1:] 进行填充。dt[-1:] 使用切片语法返回包含最后一个元素的数组（而非标量），便于直接拼接。

    phi = np.cumsum(dt / (irr + 1e-12))

计算累积相位。dt / (irr + 1e-12) 计算每个时间步内的相位增量：时间步长除以当前 RR 间期，得到该时间段内完成的心跳周期比例。分母加 10^{-12} 防止除零。np.cumsum 计算累积和，得到从起点到每个时刻的总相位。

相位的物理意义是：从仿真开始到当前时刻，累计完成了多少个心跳周期。例如 phi[i] = 5.3 表示到时刻 t[i] 为止已经完成了 5.3 个心跳周期，即发生了 5 次完整心跳，第 6 次心跳进行到 30% 的位置。

    beats = np.where(np.diff(np.floor(phi)) > 0)[0] + 1

这行代码是检测心跳发生时刻的核心。np.floor(phi) 对相位取整，得到每个时刻已完成的完整心跳数。np.diff 计算相邻元素的差值，当差值大于 0 时，表示完整心跳数增加了（从 n 变为 n+1），即该时刻发生了一次心跳。

np.where(...)[0] 返回满足条件的索引数组。由于 np.diff 使结果长度减 1 且索引对应的是差分的后一个元素，加 1 进行校正，使 beats 中的索引指向心跳实际发生的时刻。

    rr_t = t[beats]
    rr_val = irr[beats]

使用检测到的心跳索引提取心跳发生的时间戳 rr_t 和对应时刻的瞬时 RR 值 rr_val。注意这里使用瞬时 RR 值而非计算相邻心跳的实际时间差，在连续心率信号的情况下，两者近似相等。

    mask = (rr_val > 0.25) & (rr_val < 2.5)
    return rr_t[mask], rr_val[mask]

创建布尔掩码过滤生理不合理的 RR 值。RR < 0.25 秒对应心率 > 240 BPM，RR > 2.5 秒对应心率 < 24 BPM，这些都超出了正常生理范围。使用掩码索引返回过滤后的时间戳和 RR 值数组。

hrv_td 函数：时域 HRV 指标计算

def hrv_td(rr: np.ndarray) -> Dict[str, float]:
    rr_m = rr * 1000.0

将 RR 间期从秒转换为毫秒。HRV 领域的惯例是使用毫秒作为单位，这样 SDNN、RMSSD 等指标的数值范围更直观（通常在几十到几百毫秒）。

    df = np.diff(rr_m)

计算相邻 RR 间期的差值，结果长度为 len(rr_m) - 1。这个差值数组用于计算 RMSSD 和 pNN50 两个涉及逐拍变化的指标。

    return {
        "mean_rr_ms": float(np.mean(rr_m)) if len(rr_m) else np.nan,

计算平均 RR 间期。条件表达式 if len(rr_m) else np.nan 处理空数组的情况：如果没有有效的 RR 值，返回 NaN（Not a Number）而非引发错误。float() 确保返回 Python 原生浮点数而非 NumPy 标量，便于后续 JSON 序列化等操作。

        "sdnn_ms": float(np.std(rr_m, ddof=1)) if len(rr_m) > 1 else np.nan,

计算 SDNN。np.std 默认计算总体标准差（分母为 N），设置 ddof=1（delta degrees of freedom）后使用样本标准差公式（分母为 N-1）。条件 len(rr_m) > 1 确保至少有两个数据点，否则标准差无意义。

        "rmssd_ms": float(np.sqrt(np.mean(df**2))) if len(df) else np.nan,

计算 RMSSD。df**2 计算差值的平方，np.mean 取平均，np.sqrt 开方，完整实现了均方根的计算。公式展开为：

        "pnn50": float(np.mean(np.abs(df) > 50.0)) if len(df) else np.nan,

计算 pNN50。np.abs(df) > 50.0 返回布尔数组，每个元素表示对应的差值绝对值是否大于 50 ms。布尔数组可以直接参与数学运算，True 被视为 1，False 被视为 0，因此 np.mean 计算的就是 True 的比例，即满足条件的差值占总差值数的比例。

        "n_beats": int(len(rr_m)),
    }

记录心跳总数。int() 转换确保返回 Python 原生整数。心跳数是判断数据量是否充足的参考：一般认为 HRV 分析至少需要 5 分钟（约 300-400 次心跳）的数据才能得到可靠的结果。

rr2uni 函数：RR 序列重采样

def rr2uni(rr_t: np.ndarray, rr: np.ndarray, fs: float = 4.0) -> Tuple[np.ndarray, np.ndarray]:
    t0, t1 = float(rr_t[0]), float(rr_t[-1])

提取 RR 事件序列的时间范围：从第一个心跳到最后一个心跳。float() 转换将 NumPy 标量转为 Python 浮点数。

    tu = np.arange(t0, t1, 1.0/fs)

生成等间隔时间轴。np.arange(t0, t1, 1.0/fs) 从 t0 开始，以 1.0/fs 为步长，生成不超过 t1 的等差数列。默认 fs=4.0 Hz，步长为 0.25 秒。

    xu = np.interp(tu, rr_t, rr)
    return tu, xu

np.interp(tu, rr_t, rr) 进行线性插值。第一个参数 tu 是要插值的目标点，第二第三个参数 rr_t 和 rr 分别是已知数据的 x 坐标（时间戳）和 y 坐标（RR 值）。函数返回 tu 各点处的插值结果 xu，与 tu 等长。

calc_coh 函数：相干性分析

def calc_coh(rr_t: np.ndarray, rr: np.ndarray, f_tgt: float = 0.1) -> Dict[str, float]:
    if len(rr) < 20:
        return {"coherence": np.nan, "f_dom": np.nan, "pwr_tgt": np.nan, "pwr_tot": np.nan}

数据量检查。频域分析需要足够的数据点才能产生有意义的结果。20 次心跳（约 20 秒）是最低要求，低于此阈值直接返回 NaN 字典。

    tu, xu = rr2uni(rr_t, rr, fs=4.0)
    xu = znorm(xu)

调用重采样函数将不等间隔的 RR 序列转换为 4 Hz 等间隔信号，然后进行 z-score 标准化。标准化的作用是去除直流分量（均值）并归一化方差，使后续的功率谱分析专注于信号的波动特征而非绝对水平。

    f, pxx = welch(xu, fs=4.0, nperseg=min(512, len(xu)))

调用 scipy.signal.welch 计算功率谱密度。fs=4.0 指定采样频率，用于将频率轴从归一化频率转换为实际频率（Hz）。nperseg 是每段的长度，设为 512 或信号长度的较小值，避免信号过短时分段失败。

返回的 f 是频率轴数组，pxx 是对应的功率谱密度数组。频率范围从 0 到 Nyquist 频率（2 Hz），频率分辨率取决于 nperseg 和采样频率。

    pwr_tot = float(np.trapz(pxx, f))

计算总功率。np.trapz(pxx, f) 使用梯形法对功率谱密度进行数值积分。第一个参数是被积函数值（PSD），第二个参数是积分变量值（频率）。积分结果的单位是功率（PSD 的单位是功率/Hz，乘以 Hz 后得到功率）。

    band = (max(0.01, f_tgt - 0.02), f_tgt + 0.02)

定义目标频段。以目标频率为中心，向两侧各扩展 0.02 Hz。max(0.01, ...) 确保下界不小于 0.01 Hz，避免当目标频率很低时下界变为负数。

    m = (f >= band[0]) & (f <= band[1])
    pwr_tgt = float(np.trapz(pxx[m], f[m])) if np.any(m) else 0.0

创建布尔掩码选取目标频段内的频率点，然后仅对这些点进行积分得到目标频段功率。np.any(m) 检查掩码是否至少有一个 True，避免空数组导致积分失败。

    f_dom = float(f[np.argmax(pxx)])

提取主导频率。np.argmax(pxx) 返回 PSD 最大值的索引，f[...] 取出对应的频率值。

    coh = float(pwr_tgt / (pwr_tot + 1e-12))
    return {"coherence": coh, "f_dom": f_dom, "pwr_tgt": pwr_tgt, "pwr_tot": pwr_tot}

计算相干性分数：目标频段功率除以总功率。分母加 10^{-12} 防止总功率为零时除零。返回的字典包含相干性、主导频率、目标功率和总功率四个指标，便于后续分析和调试。

全流程演示代码解析

本段代码将前面定义的所有模块串联起来，形成完整的仿真管线：从呼吸波形输入到 HRV 指标输出。通过这段代码可以验证整条管线的各个环节能够正确对接和运行。

时间轴与呼吸波形生成

fs = 50
t = gen_taxis(180, fs)
vt = gen_breath(t, kind="coherence", bpm=6, amp=1.0)

首先设置采样频率为 50 Hz，这是一个适中的选择：足够高以捕捉心率的快速变化，又不会产生过多的计算负担。调用 gen_taxis 生成 180 秒（3 分钟）的时间轴，采样点数为 \lfloor 180 \times 50 \rfloor + 1 = 9001 个点。

3 分钟的仿真时长是经过考虑的：在 6 BPM 的呼吸频率下，180 秒包含 18 个完整的呼吸周期，足够观察呼吸与心率的耦合模式；同时约产生 200 次心跳（假设平均心率 70 BPM），为 HRV 分析提供了基本够用的数据量。

gen_breath 生成心脏相干性呼吸波形：6 BPM 的正弦波，幅度为 1。输出 vt 是与 t 等长的数组，数值在 [-1, 1] 范围内波动。

呼吸机械模型

res_r = RespModelA(fs=fs).simulate(t, vt)

创建 RespModelA 实例并调用 simulate 方法。构造函数参数 fs=fs 将采样频率传入模型，用于内部的低通滤波和时间延迟计算。simulate 方法接受时间轴和呼吸波形，返回 SimRes 对象 res_r。

res_r.signals 字典包含三个信号：原始输入 "vt"、延髓神经信号 "mn_med" 和脑桥神经信号 "mn_pons"。后两者将作为心血管模型的输入。

心血管模型

cres = ToyCardio().simulate(
    t,
    res_r.signals["mn_med"],
    res_r.signals["mn_pons"],
    solver="RK45",
    max_step=None,
)

创建 ToyCardio 实例（使用默认参数：静息心率 70 BPM，静息血压 100 mmHg）并调用 simulate 方法。输入包括时间轴和两个神经信号，后者通过 res_r.signals 字典访问。

solver="RK45" 指定使用 Runge-Kutta 4/5 阶方法求解 ODE，这是 SciPy 的默认算法，适用于大多数非刚性问题。max_step=None 表示不限制最大步长，让求解器自适应选择步长以满足误差要求。

返回的 cres 是 SimRes 对象，其 signals 字典包含心率 "hr" 和血压 "bp" 两个信号。

状态检查

print("OK: cres 已创建。signals =", list(cres.signals.keys()))
print("HR(min,max) =", float(np.min(cres.signals["hr"])), float(np.max(cres.signals["hr"])))

这两行打印语句用于验证心血管仿真是否成功完成。第一行确认 cres 对象存在且包含预期的信号键；第二行显示心率的最小值和最大值，可以直观判断心率波动是否合理。

在 6 BPM 呼吸驱动下，预期心率在静息值 70 BPM 附近波动，幅度约 \pm 5-8 BPM。如果看到心率范围例如 63-77 BPM，说明模型产生了预期的呼吸性心律不齐现象。

HRV 分析

rr_t, rr = hr2rr(t, cres.signals["hr"])

调用 hr2rr 函数将连续心率信号转换为 RR 间期事件序列。输入是时间轴 t 和心率数组 cres.signals["hr"]，输出是两个数组：心跳发生的时刻 rr_t 和对应的 RR 间期值 rr。

mets = {}
mets |= hrv_td(rr)
mets |= calc_coh(rr_t, rr, f_tgt=0.1)

创建空字典 mets 用于收集所有 HRV 指标。|= 是 Python 3.9 引入的字典合并运算符，等价于 mets.update(hrv_td(rr))。通过两次合并操作，将时域指标和频域相干性指标都添加到 mets 字典中。

hrv_td(rr) 返回包含 Mean RR、SDNN、RMSSD、pNN50 和心跳数的字典。calc_coh(rr_t, rr, f_tgt=0.1) 返回包含相干性分数、主导频率、目标功率和总功率的字典。f_tgt=0.1 指定目标频率为 0.1 Hz，与 6 BPM 的呼吸频率对应。

pd.DataFrame([mets])

将指标字典转换为 Pandas DataFrame 并显示。[mets] 将字典包装成单元素列表，使其成为 DataFrame 的一行数据，字典的键成为列名。在 Jupyter 环境中，这行代码会渲染出一个漂亮的表格，展示所有计算得到的 HRV 指标。

预期输出

运行这段代码后，预期看到的输出包括两行打印信息和一个指标表格。打印信息确认仿真成功，表格显示各项 HRV 指标的数值。

在心脏相干性呼吸条件下，预期的指标特征包括：相干性分数（coherence）较高，接近或超过 0.5，表示 HRV 能量集中在 0.1 Hz 附近；主导频率（f_dom）接近 0.1 Hz，与呼吸频率一致；SDNN 和 RMSSD 反映心率变异的幅度。

如果相干性很低或主导频率偏离 0.1 Hz，可能需要检查呼吸波形是否正确、模型参数是否合理、仿真时长是否足够等因素。

Cell 7：同步指标 + 完整管线 + GUI 界面

本 Cell 是整个项目的集成模块，将前面定义的所有组件整合为一个完整的应用系统。主要包含三个部分：用于评估信号同步性的指标函数、封装整条仿真流程的管线函数，以及基于 ipywidgets 的交互式图形界面。

相位锁定值（PLV）

相位锁定值（Phase Locking Value, PLV）是神经科学和生理信号分析中常用的同步性指标，用于量化两个振荡信号之间的相位耦合程度。与简单的相关系数不同，PLV 专注于相位关系而忽略幅度，能够检测出两个信号是否以稳定的相位差共同振荡。

瞬时相位的提取

计算 PLV 的第一步是从实数信号中提取瞬时相位。这通过 Hilbert 变换实现。对于实信号 x(t)，其 Hilbert 变换 \hat{x}(t) 定义为：

其中 P.V. 表示 Cauchy 主值积分。Hilbert 变换在频域的效果是将所有正频率分量的相位移动 -90°，负频率分量移动 +90°，而幅度保持不变。

将原始信号与其 Hilbert 变换组合，得到解析信号（Analytic Signal）：

解析信号是一个复数信号，其模 A(t) = |z(t)| 称为瞬时幅度，幅角 \phi(t) = \arg(z(t)) 称为瞬时相位。通过 np.angle 函数可以提取瞬时相位，再用 np.unwrap 进行相位解缠绕，消除相位从 \pi 跳变到 -\pi 的不连续性。

PLV 的计算

设两个信号的瞬时相位分别为 \phi_x(t) 和 \phi_y(t)，它们的相位差为 \Delta\phi(t) = \phi_x(t) - \phi_y(t)。PLV 定义为相位差的复指数的时间平均的模：

这个公式的几何意义非常直观。e^{j\Delta\phi(n)} 是单位圆上的一个点，其幅角等于相位差 \Delta\phi(n)。对所有时刻的这些点取平均，得到一个复数向量。如果相位差在整个时间段内保持恒定（完美相位锁定），所有点都指向同一方向，平均后的向量模长为 1；如果相位差随机变化（完全无同步），点均匀分布在单位圆上，相互抵消，平均向量的模长趋近于 0。

因此，PLV 的取值范围是 [0, 1]：

• PLV = 1 表示完美的相位锁定，两个信号的相位差始终恒定
• PLV = 0 表示完全没有相位同步，相位差随机分布
• 中间值表示部分同步，PLV 越高同步性越强

在本项目中，计算呼吸信号 V_t 与心率信号 HR 之间的 PLV，可以量化呼吸对心率调制的稳定性。心脏相干性呼吸的目标之一就是提高这种同步性。

互相关最佳延迟

互相关（Cross-Correlation）分析用于探测两个信号之间的时间延迟关系。在生理系统中，输入信号到输出响应之间通常存在一定的延迟：呼吸改变后，心率的相应变化需要经过神经传导和心脏响应等过程，不会立即发生。

延迟扫描

互相关分析的基本思路是：将一个信号相对于另一个信号在时间上进行平移，计算每个延迟值下两个信号的相关系数，然后找到相关系数最大的延迟值，即为最佳延迟。

设两个信号为 x(t) 和 y(t)，延迟为 \tau，延迟后的相关系数为：

在离散实现中，延迟 \tau 对应的采样点偏移为 k = \tau \times f_s。对于每个整数偏移 k，将信号 x 和 y 对齐后计算 Pearson 相关系数。由于偏移后信号的重叠部分长度会减少，当偏移过大时重叠太短，相关系数不可靠，因此需要设置最大搜索延迟（如 \pm 5 秒）。

结果解释

最佳延迟值的符号和大小都有生理意义。正延迟表示 x 滞后于 y，负延迟表示 x 领先于 y。在呼吸-心率分析中，通常预期心率变化略滞后于呼吸变化，最佳延迟为小的正值（如 0.5-2 秒）。

最佳延迟处的相关系数反映了两个信号的线性相关强度。相关系数接近 1 或 -1 表示强相关，接近 0 表示弱相关或无线性关系。

完整管线函数 run_sim

run_sim 函数将整条仿真流程封装为单一函数调用，接受所有可配置参数，返回完整的仿真结果。这种封装设计有两个主要好处：一是便于 GUI 回调函数调用，只需收集控件值并传入函数；二是便于批量实验，可以用循环或参数网格搜索的方式调用函数，比较不同参数组合的效果。

函数内部的执行流程严格按照仿真管线的顺序进行。首先根据时长和采样率生成时间轴，然后根据波形类型和参数生成呼吸波形。接着根据选择的呼吸模型（A 或 B）创建相应的模型实例并运行仿真。再根据选择的心血管模型（Toy、Ott97 或 SH95）创建实例并运行仿真。最后从心率信号提取 RR 序列，计算时域 HRV 指标、频域相干性和同步性指标。

函数返回一个字典，包含所有中间信号（呼吸、神经信号、心率、血压、RR 序列）、所有计算指标以及元信息（使用的模型和求解器）。这种丰富的返回值设计便于后续的深入分析和可视化。

结果展示函数 show_res

show_res 函数接受 run_sim 的返回字典，生成可视化图表并显示指标表格。图表采用堆叠子图形式，从上到下依次展示呼吸波形、延髓神经信号、脑桥神经信号、心率和血压，便于观察信号沿管线的演变过程以及各阶段之间的时序关系。

指标表格使用 Pandas DataFrame 格式化显示，将所有 HRV 指标和同步指标整合在一行中，列名清晰标识各指标含义。

GUI 界面

GUI 界面使用 ipywidgets 库构建，提供了无需修改代码即可调整参数的交互方式。整个界面分为三列布局，逻辑清晰。

左列：仿真与呼吸参数

左列包含基本仿真参数和呼吸波形参数。仿真时长滑块范围 30-600 秒，默认 180 秒，允许用户根据需要选择短时快速测试或长时详细分析。采样频率滑块范围 10-200 Hz，默认 50 Hz。

呼吸波形参数包括：波形类型下拉框（coherence/sine/square/rect）、呼吸频率滑块（2-12 BPM，默认 6）、幅度滑块、基线滑块、占空比滑块（仅对矩形波有效）、边缘平滑参数滑块。这些控件覆盖了 gen_breath 函数的所有参数。

中列：模型选择

中列包含两组切换按钮，分别用于选择呼吸模型（A 或 B）和心血管模型（Toy、Ott97 或 SH95）。切换按钮比下拉框更直观，选项数量少时操作更便捷。

由于 Ott97 和 SH95 尚未实现，选择这些模型并点击运行会看到 NotImplementedError 提示，而非程序崩溃。

右列：求解器参数与运行按钮

右列包含 ODE 求解器的配置参数：求解器算法下拉框（RK45/Radau/BDF）、相对误差容限下拉框、绝对误差容限下拉框、最大步长下拉框。这些选项允许高级用户根据模型特性调整数值精度。

Run 按钮触发仿真执行。按钮采用蓝色主样式（button_style="primary"），在界面中突出显示。点击后，回调函数收集所有控件的当前值，调用 run_sim 函数执行仿真，然后调用 show_res 显示结果。

输出区域

界面下方是输出区域（widgets.Output），用于显示图表和指标表格。每次点击 Run 时，先清除之前的输出（clear_output），再显示新的结果，避免输出累积导致页面过长。

环境检测

GUI 代码包裹在 if gui_ok: 条件块中。gui_ok 是在 Cell 1 中设置的标志位，表示 ipywidgets 是否成功导入。如果 ipywidgets 不可用（例如在某些不支持交互控件的环境中），条件为 False，跳过 GUI 创建代码，仅打印提示信息，避免程序报错。

Cell 7：同步指标 + 完整管线 + GUI 界面代码实现

本部分对 Cell 7 中同步指标函数、完整管线函数和 GUI 界面的代码进行逐行解析，重点说明 PLV 计算的实现细节、管线函数的参数传递逻辑以及 ipywidgets 界面的构建方式。

calc_plv 函数：相位锁定值计算

def calc_plv(x: np.ndarray, y: np.ndarray) -> float:
    xh = hilbert(znorm(x))
    yh = hilbert(znorm(y))

函数接受两个等长的时间序列，返回它们之间的 PLV 值。首先对两个信号进行 z-score 标准化，消除幅度和直流偏置的影响，使 PLV 仅反映相位关系。然后调用 scipy.signal.hilbert 计算 Hilbert 变换，将实信号转换为解析信号（复数形式）。

hilbert 函数的返回值是复数数组，其实部等于原始信号，虚部等于原始信号的 Hilbert 变换。解析信号的模是瞬时幅度，幅角是瞬时相位。

    ph_x = np.unwrap(np.angle(xh))
    ph_y = np.unwrap(np.angle(yh))

np.angle(xh) 提取复数数组的幅角（相位），返回值在 [-\pi, \pi] 范围内。当相位从接近 \pi 变化到接近 -\pi 时，会出现约 2\pi 的跳变，这是相位的周期性导致的，而非真实的相位突变。

np.unwrap 进行相位解缠绕，检测并消除这种 2\pi 跳变。算法检查相邻相位差，当差值超过阈值（默认 \pi）时，通过加减 2\pi 的整数倍使相位变化连续。解缠绕后的相位可以超出 [-\pi, \pi] 范围，但保持了真实的相位演变趋势。

    d = ph_x - ph_y
    return float(np.abs(np.mean(np.exp(1j * d))))

计算两个信号的瞬时相位差 d。np.exp(1j * d) 将相位差转换为单位圆上的复数点：相位差为 \theta 时，对应复数 e^{j\theta} = \cos\theta + j\sin\theta，模长为 1，幅角为 \theta。

np.mean 对所有时刻的复数点取平均。如果相位差 d 在整个时间段内保持恒定（假设为 \theta_0），则所有点都是 e^{j\theta_0}，平均值仍是 e^{j\theta_0}，模长为 1。如果相位差随机变化，点在单位圆上均匀分布，平均值趋向原点，模长趋向 0。

np.abs 取复数的模长，即 PLV 值。float() 转换确保返回 Python 原生浮点数。

xcorr_lag 函数：互相关延迟分析

def xcorr_lag(x: np.ndarray, y: np.ndarray, fs: float, max_lag_s: float = 5.0) -> Dict[str, float]:
    x = znorm(x)
    y = znorm(y)
    max_k = int(max_lag_s * fs)
    lags = np.arange(-max_k, max_k + 1)

函数首先对两个信号进行标准化，确保后续计算的相关系数在 [-1, 1] 范围内。max_k 是最大延迟对应的采样点数，例如 max_lag_s=5.0 秒、fs=50 Hz 时，max_k=250 个点。

lags 是延迟数组，从 -max_k 到 max_k（包含两端），共 2 \times \text{max\_k} + 1 个延迟值。负延迟表示将 x 向前移动（或等效地将 y 向后移动）。

    corr = []
    for k in lags:
        xs = x[max(0, k):len(x)+min(0, k)]
        ys = y[max(0, -k):len(y)+min(0, -k)]

循环遍历每个延迟值 k，计算该延迟下的相关系数。信号对齐的索引计算是这段代码的核心难点。

当 k > 0 时，x 相对于 y 向后移动 k 个点。为了对齐，取 x[k:] 和 y[:-k]（或等效地 y[:len(y)-k]）。代入公式验证：max(0, k) = k，len(x) + min(0, k) = len(x)，所以 xs = x[k:]；max(0, -k) = 0，len(y) + min(0, -k) = len(y) - k，所以 ys = y[:len(y)-k]，即 y[:-k]。

当 k < 0 时，x 相对于 y 向前移动 |k| 个点。代入公式：max(0, k) = 0，len(x) + min(0, k) = len(x) + k = len(x) - |k|，所以 xs = x[:len(x)-|k|]，即 x[:k]；max(0, -k) = |k|，len(y) + min(0, -k) = len(y)，所以 ys = y[|k|:]，即 y[-k:]。

当 k = 0 时，xs = x[0:len(x)] = x，ys = y[0:len(y)] = y，两个完整信号。

        if len(xs) < 10:
            corr.append(np.nan)
        else:
            corr.append(np.corrcoef(xs, ys)[0, 1])

当延迟较大时，重叠部分变短。如果重叠长度小于 10 个点，认为数据量不足以计算可靠的相关系数，记录 NaN。否则调用 np.corrcoef 计算 Pearson 相关系数。np.corrcoef 返回 2 \times 2 的相关矩阵，对角线元素为 1（自相关），非对角线元素为两个变量之间的相关系数，取 [0, 1] 位置的元素。

    corr = np.asarray(corr)
    idx = np.nanargmax(corr)
    return {"best_lag_s": float(lags[idx] / fs), "best_corr": float(corr[idx])}

将相关系数列表转换为数组，然后用 np.nanargmax 找到最大值的索引（忽略 NaN 值）。lags[idx] 是最佳延迟的采样点数，除以采样频率 fs 转换为秒。返回字典包含最佳延迟（秒）和对应的相关系数。

run_sim 函数：完整仿真管线

def run_sim(
    dur: float, fs: float, pat: str, bpm: float, amp: float, base: float,
    duty: float, edge_s: float, resp_k: str, cardio_k: str,
    solver: str, rtol: float, atol: float, max_step: Optional[float],
) -> Dict[str, Any]:

函数签名包含 14 个参数，涵盖了仿真管线中所有可配置的选项。参数名采用简短形式（如 dur 而非 duration_s），与前面代码重构的风格保持一致。返回类型为包含所有结果的字典。

    t = gen_taxis(dur, fs)
    vt = gen_breath(t, kind=pat, bpm=bpm, amp=amp, base=base, duty=duty, edge_s=edge_s)

步骤 1：生成时间轴和呼吸波形。参数直接传递给相应的函数，kind=pat 使用关键字参数形式，增强代码可读性。

    if resp_k == "A":
        res_r = RespModelA(fs=fs).simulate(t, vt)
    elif resp_k == "B":
        res_r = RespModelB().simulate(t, vt, solver=solver)
    else:
        raise ValueError("Unknown resp model")

步骤 2：根据选择创建并运行呼吸模型。模型 A 需要采样频率参数，模型 B 需要求解器参数。两者的 simulate 方法签名不同，但都返回 SimRes 对象。如果传入未知的模型标识，抛出 ValueError。

    if cardio_k == "Toy":
        card = ToyCardio()
    elif cardio_k == "Ott97":
        card = Ott97()
    elif cardio_k == "SH95":
        card = SH95()
    else:
        raise ValueError("Unknown cardio model")

步骤 3：根据选择创建心血管模型实例。三个模型都使用默认参数构造。Ott97 和 SH95 虽然可以创建实例，但调用 simulate 时会抛出 NotImplementedError。

    cres = card.simulate(
        t, res_r.signals["mn_med"], res_r.signals["mn_pons"],
        solver=solver, rtol=rtol, atol=atol, max_step=max_step
    )
    hr = cres.signals["hr"]
    bp = cres.signals["bp"]

步骤 4：运行心血管模型。输入包括时间轴和来自呼吸模型的两个神经信号，ODE 求解器参数也一并传入。提取心率和血压信号供后续分析使用。

    rr_t, rr = hr2rr(t, hr)
    mets = {}
    mets |= hrv_td(rr)
    mets |= calc_coh(rr_t, rr, f_tgt=bpm/60.0)

步骤 5：HRV 分析。将心率转换为 RR 序列，计算时域指标和频域相干性。相干性分析的目标频率设为 bpm/60.0，自动匹配当前的呼吸频率。例如呼吸 6 BPM 时，目标频率为 0.1 Hz。

    mets["plv(vt,hr)"] = calc_plv(vt, hr)
    mets |= xcorr_lag(vt, hr, fs=fs, max_lag_s=10)

步骤 6：同步指标。计算呼吸与心率之间的 PLV，直接赋值给字典的新键。计算互相关最佳延迟，最大搜索范围设为 \pm 10 秒，返回的字典通过 |= 合并到 mets 中。

    return {
        "t": t, "vt": vt,
        "mn_med": res_r.signals["mn_med"],
        "mn_pons": res_r.signals["mn_pons"],
        "hr": hr, "bp": bp,
        "rr_t": rr_t, "rr": rr,
        "metrics": mets,
        "meta": {"resp": resp_k, "cardio": cardio_k, "solver": solver}
    }

返回包含所有结果的字典。信号数据（时间轴、呼吸、神经信号、心率、血压、RR 序列）便于后续绑图；指标字典 mets 便于表格显示；元信息 meta 记录了使用的模型和求解器，便于追溯。

show_res 函数：结果展示

def show_res(out: Dict[str, Any]):
    plot_sigs(out["t"], {
        "Vt": out["vt"],
        "MN_med": out["mn_med"],
        "MN_pons": out["mn_pons"],
        "HR (bpm)": out["hr"],
        "BP (a.u.)": out["bp"],
    }, title=f"Simulation meta: {out['meta']}")
    display(pd.DataFrame([out["metrics"]]))

函数接受 run_sim 的返回字典，调用 plot_sigs 绑制五个信号的堆叠图，标题显示元信息。然后将指标字典转换为 DataFrame 并通过 display 显示。在 Jupyter 环境中，display 会渲染出格式化的表格。

GUI 控件定义

if gui_ok:
    w_dur = widgets.FloatSlider(description="Duration(s)", min=30, max=600, step=10, value=180)
    w_fs = widgets.IntSlider(description="fs(Hz)", min=10, max=200, step=10, value=50)

GUI 代码包裹在 if gui_ok: 条件中，确保 ipywidgets 可用时才执行。滑块控件通过 widgets.FloatSlider 或 widgets.IntSlider 创建，参数包括描述文字（显示在控件旁边）、最小值、最大值、步长和默认值。

FloatSlider 产生浮点数值，IntSlider 产生整数值。采样频率使用整数滑块，因为采样频率通常取整数值。

    w_pat = widgets.Dropdown(description="Pattern", options=["coherence", "sine", "square", "rect"], value="coherence")

下拉框通过 widgets.Dropdown 创建，options 是可选项列表，value 是默认选中项。下拉框适合选项数量有限且互斥的场景。

    w_resp = widgets.ToggleButtons(description="Resp", options=["A", "B"], value="A")
    w_cardio = widgets.ToggleButtons(description="Cardio", options=["Toy", "Ott97", "SH95"], value="Toy")

切换按钮通过 widgets.ToggleButtons 创建，选项以并排按钮形式显示，当前选中项高亮。切换按钮比下拉框更直观，适合选项少（2-5 个）的场景。

    w_rtol = widgets.Dropdown(description="rtol", options=[1e-3, 1e-5, 1e-6, 1e-7], value=1e-6)

误差容限使用下拉框而非滑块，因为这些参数通常取特定的数量级值（如 10^{-6}），滑块难以精确选择这类值。options 列表中直接使用科学记数法的浮点数。

    w_maxs = widgets.Dropdown(description="max_step", options=[None, 0.2, 0.1, 0.05, 0.02], value=None)

最大步长的选项包含 None，表示不限制步长。下拉框可以包含任意 Python 对象作为选项值，不限于字符串。

    btn_run = widgets.Button(description="Run", button_style="primary")
    out_box = widgets.Output()

按钮通过 widgets.Button 创建，button_style="primary" 设置蓝色主题样式。输出区域通过 widgets.Output 创建，它是一个特殊的容器控件，可以捕获并显示打印输出、图表等内容。

回调函数与事件绑定

    def on_click_run(_):
        with out_box:
            clear_output(wait=True)
            res = run_sim(
                dur=float(w_dur.value),
                fs=float(w_fs.value),
                # ... 其他参数 ...
                max_step=w_maxs.value,
            )
            show_res(res)

    btn_run.on_click(on_click_run)

回调函数 on_click_run 在按钮点击时被调用。参数 _ 是按钮事件对象，这里不使用，用下划线命名表示忽略。

with out_box: 是上下文管理器语法，在此块内的所有输出（print、图表等）都会被重定向到 out_box 控件中显示，而非直接输出到 notebook。

clear_output(wait=True) 清除 out_box 中之前的内容。wait=True 表示等待新内容准备好后再清除旧内容，避免闪烁。

回调函数从各控件的 .value 属性获取当前值，组装参数调用 run_sim，然后调用 show_res 显示结果。float() 和 str() 转换确保类型正确。w_maxs.value 可能是 None 或浮点数，直接传递无需转换。

btn_run.on_click(on_click_run) 将回调函数绑定到按钮的点击事件。每次点击按钮时，回调函数被调用一次。

布局与显示

    display(widgets.HBox([
        widgets.VBox([w_dur, w_fs, w_pat, w_bpm, w_amp, w_base, w_duty, w_edge]),
        widgets.VBox([w_resp, w_cardio]),
        widgets.VBox([w_solver, w_rtol, w_atol, w_maxs, btn_run]),
    ]))
    display(out_box)

VBox 将控件垂直排列，HBox 将多个 VBox 水平排列，形成三列布局。第一列包含 8 个仿真和呼吸参数控件，第二列包含 2 个模型选择控件，第三列包含 4 个求解器参数控件和运行按钮。

两次 display 调用分别显示控件面板和输出区域。输出区域位于控件下方，点击 Run 后结果会在此处显示。

else:
    print("ipywidgets 不可用：请先安装 ipywidgets。")

如果 gui_ok 为 False（ipywidgets 导入失败），跳过所有 GUI 代码，仅打印提示信息。这确保代码在不支持交互控件的环境中也能运行而不报错。