自主心率控制的频率分析
注意事项
在这个例子中，术语"迷走神经活动"用作"迷走神经活动"（来自神经迷走神经-迷走神经）的同义词，并且意味着自主神经系统的副交感神经分支的活动。
使用的缩写列表
BP-血压
Dsa-呼吸窦性心律失常
LFM-线性频率调制
SA-Sinoatrial
低通低通滤波器
HR-心率
ABP-动脉血压
HR-心率
RSA-呼吸窦性心律失常
附加功能
using EngeeDSP
要运行模型，请使用附加函数run_model。
# 启动模型的功能
function run_model( name_model, path_to_folder ) # 定义用于运行模型的函数
Path = path_to_folder * "/" * name_model * ".engee"
if name_model in [m.name for m in engee.get_all_models()] # 检查将模型加载到内核的条件
model = engee.open( name_model ) # 打开模型
model_output = engee.run( model, verbose=true ); # 启动模型
engee.close( name_model, force=true ); # 关闭模型
else
model = engee.load( Path, force=true ) # 上传模型
model_output = engee.run( model, verbose=true ); # 启动模型
engee.close( name_model, force=true ); # 关闭模型
end
return model_output
end;
对于结果的频率分析，我们将使用freq_analysis函数。
function freq_analysis(V, HR, Fs)
n = length(V)
# 汉娜窗口的应用
window = EngeeDSP.Functions.hann(n)
V_win = V .* window
HR_win = HR .* window
# FFT计算
fft_V = EngeeDSP.rfft(V_win) / n
fft_HR = EngeeDSP.rfft(HR_win) / n
freq = EngeeDSP.rfftfreq(n, Fs)
# 互谱
power_V = abs2.(fft_V) * (n / sum(window.^2))
cross_power = fft_HR .* conj(fft_V) * (n / sum(window.^2))
# 传递函数
H = @. cross_power / (power_V + eps())
amplitude = abs.(H)
phase_deg = EngeeDSP.rad2deg.(EngeeDSP.angle.(H))
# 限制为0.3赫兹
max_index = findlast(f -> f <= 0.3, freq) # 查找频率<=0.3Hz的最后一个索引
freq = freq[1:max_index]
amplitude = amplitude[1:max_index]
phase_deg = phase_deg[1:max_index]
return amplitude, phase_deg, freq
end;
导言
心率（HR）和血压（BP）的短期调节主要是由于通过动脉压力感受器的反馈而进行的，这两个参数都不断受到呼吸的影响。 呼吸以几种相互关联的方式影响心脏和血管。 在吸气和呼气期间胸部内部压力的变化直接影响血压的量;血压的这些波动反过来通过压力感受器反射改变心率。 此外，现代研究表明，延髓的呼吸中心与控制心脏工作的自主中心有直接的神经连接。 这种复杂的呼吸效应的累积结果表现为呼吸窦性心律失常（TSA）的形式，这是一种自然现象，其中心跳在吸气时略微加速，呼气时减慢。 呼吸对心率（DSA值）的这种影响的强度可以使用特殊的数学分析–频率传递函数来量化。 这种分析的特征的变化（例如，相移或幅度的变化）作为自主神经系统工作变化的指标，其调节心率。
实施血液循环模型（"Blood_circulation.engee"）如下图所示。
.png)
呼吸，测量为肺体积（V）的变化，被认为对窦房结（CA节点）的自主输入有直接影响：吸入导致迷走神经和交感神经传入活动的减少。 来自压力感受器的反馈也直接影响心脏的自主输入：血压（ABP）的增加导致交感神经活动的减少和副交感神经（迷走神经）活动的增加。
在灵感期间，迷走神经传出活动的减少作用于CA节点，增加心率（HR）。 模拟该比率动态的传递函数（游荡分支）是一个简单的低通滤波器，具有截止频率（fp）和负增益（-Kp）。 相反，CA节点对交感神经刺激的反应要慢得多。 除了1-2秒的延迟之外，描述交感神经活动转化为心率的动力学的传递函数具有截止频率(fs)。 在这种情况下，增益为正。 为了简单起见，表示动脉血管床属性的传递函数被认为是静态的，增益为0.01。 此外，由于BP（ABP）转换为压力感受器的输出信号发生在非常快的动态中，因此假设baroreflex可以用0.01的静态增益表示，以0.3s的固定延迟顺序打开。 因此，该模型通过允许直接自主刺激心率（HR）来模拟呼吸窦性心律失常。 此外，由此产生的心率（HR）的变化和呼吸的直接机械作用引起血压（ABP）的波动，随后通过血压波动影响心率（HR）。
血液循环模型的模拟
仿真参数的初始化
线性频率调制信号发生器将用于确定人体循环模型的频率响应。 在此示例中，将设置参数，使初始频率为0.005Hz，目标频率为0.5Hz。 由于LFM信号发生器单元的振幅不改变，因此需要使用0.3的增益因子，其将振幅范围限制为0.6升。
Ts = 10e-3 # 步长，带
time_simulation = 300 # 模拟时间，与
# 呼吸模式参数(啁啾)
target_time = 300 # 目标时间，从
freq_initial = 0.005 # 初始频率，Hz
freq_target = 0.5; # 目标频率，Hz
按照生理参数的时间延迟。
time_delay_sym = 2 # 交感神经影响的CA节点延迟，具有
time_delay_vas = 0.42 # 血管延迟，具有
time_delay_bar = 0.3; # Baroreflex延迟,s
参数的以下标称值表示处于仰卧位的健康受试者:
*Ca节点的传递函数对于游移影响的增益，Kp=6;
*交感神经影响的CA节点传递函数的增益Ks=18;
*游移影响的CA节点传递函数的截止频率，fp=0.2Hz;
*交感神经影响的CA节点传递函数的截止频率，fs=0.015Hz;
*呼吸驱动转化为传出神经活动的加权因子为:
*对于流浪分支：Ap=2.5,
*对于交感神经分支：As=0.4。
Kp = 6 # 游移影响的CA节点传递函数放大系数
Ks = 18 # 交感神经影响的CA节点传递函数的放大系数
Ap = 2.5 # 流浪树枝的权重系数
As = 0.4 # 交感神经分支的权重因子
fp = 0.2 # 对于游荡影响的CA节点的传递函数的截止频率，Hz
fs = 0.015; # 交感神经影响的CA节点的传递函数的截止频率，Hz
启动模型
让我们开始使用run_model函数计算模型的仿真。
out = run_model("Blood_circulation", @__DIR__); # 启动模型
提取得到的结果。
V_normal= out["V"].value;
HR_normal = out["HR"].value;
ABP_normal = out["ABP"].value;
time_steps= length(V_normal)
time_model = range(0, step=Ts, length=time_steps);
模拟结果的可视化
下图显示了仿真结果，为了更清晰地显示了仿真的前200秒。 上图显示呼吸模式，即FM信号发生器。 平均图表对应于心率的变化。 注意，在低频率下，心率几乎与肺体积的变化同步波动;然而，在较高频率下，它开始相对于呼吸滞后。 此外，HR信号的幅度随着频率的增加而减小，强调了滤波器的整体频率响应的性质。 下图显示了血压的变化。 可以注意到，随着频率增加，呼吸诱导的血压变化（ABP）变得更大。 这是呼吸频率增加的直接机械效应影响的结果。
# 结果可视化
max_time = 200 # 时间根据哦，与
plot(
# 体积图(V)
plot(time_model, V_normal;
xlabel="时间，从",
ylabel="卷，l",
label="肺容积变化(V)",
legend=:topright,
linewidth=1.5,
color=:blue,
grid=true,
xlims=(0, max_time),
ylims=(-0.4, 0.4)),
# 心率图(HR)
plot(time_model, HR_normal;
xlabel="时间，从",
ylabel="心率，心跳/分钟",
label=" 心率（HR）的变化",
legend=:topright,
linewidth=1.5,
color=:red,
grid=true,
xlims=(0, max_time),
ylims=(-10, 10)),
# 血压图(ABP)
plot(time_model, ABP_normal;
xlabel="时间，从",
ylabel="压力，mmHg",
label="血压变化(ABP)",
legend=:topright,
linewidth=1.5,
color=:green,
grid=true,
xlims=(0, max_time),
ylims=(-1, 1)),
# 整个图表的常规设置
layout = (3, 1),
size = (1500, 1000),
margin = 40 * Plots.px,
plot_title = "生理指标"
)
频率响应
Fs = 1 / Ts
# 频率特性的计算
amplitude_normal, phase_deg_normal, freq_normal = freq_analysis(V_normal,HR_normal, Fs)
plot(
plot(freq_normal, amplitude_normal;
title = "幅度-频率响应",
xlabel = "频率，赫兹",
ylabel = "/H/，min/l",
legend = false,
label = false,
linewidth = 2.5,
color = :blue,
grid = true,
xlims = (0, 0.3)),
plot(freq_normal, phase_deg_normal;
title = "相位频率响应",
xlabel = "频率，赫兹",
ylabel = "相位，°",
legend = false,
label = false,
linewidth = 2.5,
color = :red,
grid = true,
ylims = (-180, 180),
xlims = (0, 0.3)),
layout = (2, 1),
size = (1500, 1000),
margin = 40*Plots.px,
plot_title = "传递函数分析：V→HR"
)
模仿药物管理
在这一部分中，将考虑考虑药物对所研究参数的影响的两种情况。
第一种情况
这种情况模拟一剂阿托品的给药。 这种药物导致完全副交感神经（迷走神经）封锁。 此外，模型参数被改变以模拟当患者处于站立位置时，当对心率的交感神经效应增加时的情况。 在这种情况下，心率控制主要由交感神经系统调节。
初始化第一个方案的参数
用于模拟阿托品剂量给药的模型参数的值如下所示:
Kp = 1 # 游移影响的CA节点传递函数放大系数
Ks = 9 # 交感神经影响的CA节点传递函数的放大系数
Ap = 0.1 # 流浪树枝的权重系数
As = 4 # 交感神经分支的权重因子
fp = 0.07 # 对于游荡影响的CA节点的传递函数的截止频率，Hz
fs = 0.015; # 交感神经影响的CA节点的传递函数的截止频率，Hz
启动模型
out_atropine = run_model("Blood_circulation", @__DIR__); # 启动模型
V_atropine= out_atropine["V"].value;
HR_atropine = out_atropine["HR"].value;
ABP_atropine = out_atropine["ABP"].value;
第一个场景的仿真结果的可视化
# 结果可视化
max_time = 200 # 时间根据哦，与
plot(
# 体积图(V)
plot(time_model, V_atropine;
xlabel="时间，从",
ylabel="卷，l",
label="肺容积变化(V)",
legend=:topright,
linewidth=1.5,
color=:blue,
grid=true,
xlims=(0, max_time),
ylims=(-0.4, 0.4)),
# 心率图(HR)
plot(time_model, HR_atropine;
xlabel="时间，从",
ylabel="心率，心跳/分钟",
label=" 心率（HR）的变化",
legend=:topright,
linewidth=1.5,
color=:red,
grid=true,
xlims=(0, max_time),
ylims=(-10, 10)),
# 血压图(ABP)
plot(time_model, ABP_atropine;
xlabel="时间，从",
ylabel="压力，mmHg",
label="血压变化(ABP)",
legend=:topright,
linewidth=1.5,
color=:green,
grid=true,
xlims=(0, max_time),
ylims=(-1, 1)),
# 整个图表的常规设置
layout = (3, 1),
size = (1500, 1000),
margin = 40 * Plots.px,
plot_title = "阿托品给药期间的生理参数"
)
# 频率特性的计算
amplitude_atropine, phase_deg_atropine, freq_atropine = freq_analysis(V_atropine,HR_atropine, Fs)
plot(
plot(freq_atropine, amplitude_atropine;
title = "幅度-频率响应",
xlabel = "频率，赫兹",
ylabel = "/H/，min/l",
legend = false,
label = false,
linewidth = 2.5,
color = :blue,
grid = true,
xlims = (0, 0.3)),
plot(freq_atropine, phase_deg_atropine;
title = "相位频率响应",
xlabel = "频率，赫兹",
ylabel = "相位，°",
legend = false,
label = false,
linewidth = 2.5,
color = :red,
grid = true,
ylims = (-180, 180),
xlims = (0, 0.3)),
layout = (2, 1),
size = (1200, 1000),
margin = 40*Plots.px,
plot_title = "传递函数分析：V→HR"
)
第二种情况
这种情况模拟了一剂普萘洛尔的给药。 这种药物导致β-肾上腺素能阻断。 此外，我们假设仰卧位，这使得徘徊调制成为主要的控制模式。
初始化第二个方案的参数
用于模拟普萘洛尔剂量给药的模型参数的值如下所示:
Kp = 6 # 游移影响的CA节点传递函数放大系数
Ks = 1 # 交感神经影响的CA节点传递函数的放大系数
Ap = 2.5 # 流浪树枝的权重系数
As = 0.1 # 交感神经分支的权重因子
fp = 0.2 # 对于游荡影响的CA节点的传递函数的截止频率，Hz
fs = 0.015; # 交感神经影响的CA节点的传递函数的截止频率，Hz
启动模型
out_propranolol = run_model("Blood_circulation", @__DIR__); # 启动模型
V_propranolol= out_propranolol["V"].value;
HR_propranolol = out_propranolol["HR"].value;
ABP_propranolol = out_propranolol["ABP"].value;
第二场景的仿真结果的可视化
# 结果可视化
max_time = 200 # 时间根据哦，与
plot(
# 体积图(V)
plot(time_model, V_propranolol;
xlabel="时间，从",
ylabel="卷，l",
label="肺容积变化(V)",
legend=:topright,
linewidth=1.5,
color=:blue,
grid=true,
xlims=(0, max_time),
ylims=(-0.4, 0.4)),
# 心率图(HR)
plot(time_model, HR_propranolol;
xlabel="时间，从",
ylabel="心率，心跳/分钟",
label=" 心率（HR）的变化",
legend=:topright,
linewidth=1.5,
color=:red,
grid=true,
xlims=(0, max_time),
ylims=(-10, 10)),
# 血压图(ABP)
plot(time_model, ABP_propranolol;
xlabel="时间，从",
ylabel="压力，mmHg",
label="血压变化(ABP)",
legend=:topright,
linewidth=1.5,
color=:green,
grid=true,
xlims=(0, max_time),
ylims=(-1, 1)),
# 整个图表的常规设置
layout = (3, 1),
size = (1500, 1000),
margin = 40 * Plots.px,
plot_title = "丙醇给药期间的生理参数"
)
# 频率特性的计算
amplitude_propranolol, phase_deg_propranolol, freq_propranolol = freq_analysis(V_propranolol, HR_propranolol, Fs)
plot(
plot(freq_propranolol, amplitude_propranolol;
title = "幅度-频率响应",
xlabel = "频率，赫兹",
ylabel = "/H/，min/l",
legend = false,
label = false,
linewidth = 2.5,
color = :blue,
grid = true,
xlims = (0, 0.3)),
plot(freq_propranolol, phase_deg_propranolol;
title = "相位频率响应",
xlabel = "频率，赫兹",
ylabel = "相位，°",
legend = false,
label = false,
linewidth = 2.5,
color = :red,
grid = true,
ylims = (-180, 180),
xlims = (0, 0.3)),
layout = (2, 1),
size = (1200, 1000),
margin = 40*Plots.px,
plot_title = "传递函数分析：V→HR"
)
对所得结果的分析
# 3种场景的仿真结果可视化
p_amp = plot(
title = "幅度-频率响应",
xlabel = "频率，赫兹",
ylabel = "/H/，min/l",
legend = :topright,
grid = true,
ylims = (0, 80),
xlims = (0, 0.3),
size = (1200, 1000)
)
p_phase = plot(
title = "相位频率响应",
xlabel = "频率，赫兹",
ylabel = "相位，°",
legend = :topright,
grid = true,
ylims = (-180, 180),
xlims = (0, 0.3)
)
# 正常情况
plot!(p_amp, freq_normal, amplitude_normal,
label = "正常", linewidth = 1.5, color = :blue)
plot!(p_phase, freq_normal, phase_deg_normal,
label = false, linewidth = 1.5, color = :blue)
# 丙醇注射液
plot!(p_amp, freq_propranolol, amplitude_propranolol,
label = "丙醇,丙醇", linewidth = 1.5, color = :green)
plot!(p_phase, freq_propranolol, phase_deg_propranolol,
label = false, linewidth = 1.5, color = :green)
# 阿托品注射液
plot!(p_amp, freq_atropine, amplitude_atropine,
label = "阿托品", linewidth = 1.5, color = :red)
plot!(p_phase, freq_atropine, phase_deg_atropine,
label = false, linewidth = 1.5, color = :red)
# 总时间表
plot_total = plot(p_amp, p_phase,
layout = (2, 1),
plot_title = "传递函数的比较：V→HR",
margin = 40*Plots.px,
size = (1200, 1000))
# 显示图形
display(plot_total)
阿托品（阻断副交感神经）的给药主要将心率控制转移到交感神经系统。 这反映在频率响应中:
*在低频（低于0.03hz）时增强响应,
*高频响应衰减（0.1Hz以上）,
*增加相位延迟（相位曲线的陡坡）。
普萘洛尔（阻断交感神经）的给药导致副交感神经占主导地位的状态:
*低频响应严重衰减,
*频率超过0.05赫兹时的轻微变化,
*小相移（0-0.3Hz），这表明心率对呼吸的快速反应。
结论
这种血液循环模型描述了呼吸对心率和血压调节的短期影响。 它展示了呼吸过程，自主神经系统，baroreflex和对心率和血压形成的直接机械影响的相互作用。
对模型频率特征的分析使我们能够观察交感神经和副交感神经系统对心率控制的动力学和总体贡献如何在各种生理条件下变化以及各种药物的影
使用的来源列表
1-Saul，J.P.，R.D.Berger，P.Albrecht，S.P.Stein，M.H.Chen和R.J.Cohen。循环的转移功能分析:对心血管调节的独特见解. 阿姆。 1. 生理醇。 261(心脏Circ. 生理醇。 30)：HI231-HI245，1991.