Сообщество Engee

Частотный анализ вегетативного контроля сердечного ритма

Автор
avatar-akmalovmakmalovm
Соавторы
avatar-dimabalakindimabalakin
Notebook

Частотный анализ вегетативного контроля сердечного ритма

Примечание

В данном примере термин "блуждающая активность" используется как синоним "вагусной активности" (от Nervus vagus - блуждающий нерв) и означает активность парасимпатической ветви вегетативной нервной системы.

Список используемых сокращений

АД - Артериальное давление

ДСА - Дыхательная синусовая аритмия

ЛЧМ - Линейно-частотно модулированный

СА - Синоатриальный

ФНЧ - Фильтр нижних частот

ЧСС - Частота сердечных сокращений

ABP - Arterial blood pressure

HR - Heart rate

RSA - Respiratory sinus arrhythmia

Дополнительные функции

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.3 Гц
    freq = freq[1:max_index]
    amplitude = amplitude[1:max_index]
    phase_deg = phase_deg[1:max_index]
    
    return amplitude, phase_deg, freq
end;

Введение

Краткосрочная регуляция частоты сердечных сокращений (ЧСС) и артериального давления (АД) в основном осуществляется за счет обратной связи через артериальные бароцепторы, оба параметра постоянно испытывают влияния дыхания. Дыхание воздействует на сердце и сосуды несколькими взаимосвязанными путями. Изменения давления внутри грудной клетки во время вдоха и выдоха непосредственно влияют на величину артериального давления; эти колебания АД, в свою очередь, через барорецепторные рефлексы изменяют частоту сердцебиения. Помимо этого, современные исследования показывают, что дыхательный центр в продолговатом мозге имеет прямые нейронные связи с вегетативными центрами, управляющими работой сердца. Совокупный результат этого комплексного влияния дыхания проявляется в виде дыхательной синусовой аритмии (ДСА) – естественного явления, при котором сердцебиение немного ускоряется на вдохе и замедляется на выдохе. Силу этого влияния дыхания на ЧСС (величину ДСА) можно количественно оценить с помощью специального математического анализа – частотной передаточной функции. Изменения характеристик этого анализа (например, сдвиг фазы или изменение амплитуды) служат индикатором изменений в работе вегетативной нервной системы, регулирующей сердечный ритм.

Реализация модели кровообращения ("Blood_circulation.engee") показана на рисунке ниже.

eb95fbf2_20ac_4fc9_a6b4_f834e7f5b117.png

Дыхание, измеряемое как изменение объема легких (V), предполагается оказывающим прямое влияние на вегетативные входы синоатриального узла (СА-узла): вдох приводит к уменьшению как блуждающей, так и симпатической эфферентной активности. Обратная связь от барорецепторов также непосредственно влияет на вегетативные входы сердца: повышение артериального давления (ABP) вызывает снижение симпатической активности и увеличение парасимпатической (блуждающей) активности.

Во время вдоха снижение блуждающей эфферентной активности действует на СА-узел, увеличивая ЧСС (HR). Передаточная функция, моделирующая динамику этого соотношения (блуждающая ветвь), представляет собой простой фильтр нижних частот (ФНЧ) с частотой среза (fp) и отрицательным коэффициентом усиления (-Kp). Напротив, реакция СА-узла на симпатическую стимуляцию значительно медленнее. Помимо задержки в 1-2 с, передаточная функция, описывающая динамику преобразования симпатической активности в ЧСС, имеет частоту среза (fs). В этом случае коэффициент усиления положителен. Изменения ЧСС предполагается влияющими на артериальное давление после задержки в 0.42 с. Для простоты передаточная функция, представляющая свойства артериального сосудистого русла, считается статической с коэффициентом усиления 0.01. Кроме того, поскольку преобразование АД (ABP) в выходной сигнал барорецепторов происходит с очень быстрой динамикой, предполагается, что барорефлекс может быть представлен статическим коэффициентом усиления равным 0.01, включенным последовательно с фиксированной задержкой 0.3 с. Также прямое механическое влияние дыхания на АД (ABP) моделируется как отрицательный дифференциатор, то есть вдох стремится снизить АД (ABP), а выдох - повысить его. Таким образом, модель симулирует дыхательную синусовую аритмию, допуская прямое вегетативное стимулирование ЧСС (HR). Кроме того, результирующие изменения ЧСС (HR) и прямое механическое влияние дыхания вызывают колебания АД (ABP), которые впоследствии влияют на ЧСС (HR) через барорефлексы.

Симуляция модели кровообращения

Инициализация параметров моделирования

Для определения частотной характеристики модели кровообращения человека будет использоваться генератор линейно-частотно модулируемого сигнала. В данном примере будут установлены параметры таким образом, чтобы начальная частота равнялась 0.005 Гц, а целевая частота 0.5 Гц. Поскольку амплитуда блока генератора ЛЧМ сигнала не изменяется, необходимо использовать коэффициент усиления 0.3, который ограничивает размах амплитуды до 0.6 литров.

Ts = 10e-3               # Размер шага, с
time_simulation = 300    # Время моделирования, с

# Параметры патерна дыхания (Chirp)
target_time = 300       # Целевое время, с
freq_initial = 0.005    # Начальная частота, Гц
freq_target = 0.5;      # Целевая частота, Гц

Временные задержки в соответствие с физиологическими параметрами.

time_delay_sym = 2      # Задержка СА-узла для симпатического влияния, с
time_delay_vas = 0.42   # Сосудистая задержка, с
time_delay_bar = 0.3;   # Задержка барорефлекса, с

Следующие номинальные значения параметров представляют здорового испытуемого в положении лежа на спине:

  • Коэффициент усиления передаточной функции СА-узла для блуждающего влияния, Kp = 6;

  • Коэффициент усиления передаточной функции СА-узла для симпатического влияния Ks = 18;

  • Частота среза передаточной функции СА-узла для блуждающего влияния, fp = 0.2 Гц;

  • Частота среза передаточной функции СА-узла для симпатического влияния, fs = 0.015 Гц;

  • Весовые коэффициенты для преобразования дыхательного драйва в эфферентную нервную активность составляют:

    • Для блуждающей ветви: Ap = 2.5,

    • Для симпатической ветви: As = 0.4.

Kp = 6      # Коэффициентусиления передаточной функции СА-узла для блуждающего влияния
Ks = 18     # Коэффициент усиления передаточной функции СА-узла для симпатического влияния
Ap = 2.5    # Весовой коэффицент для блуждающей ветви
As = 0.4    # Весовой коэффицент для симпатической ветви
fp = 0.2    # Частота среза передаточной функции СА-узла для блуждающего влияния, Гц
fs = 0.015; # Частота среза передаточной функции СА-узла для симпатического влияния, Гц

Запуск модели

Запустим расчет симуляции модели при помощи функции run_model.

out = run_model("Blood_circulation", @__DIR__);    # Запуск модели
Building...
Progress 0%
Progress 5%
Progress 11%
Progress 17%
Progress 23%
Progress 29%
Progress 34%
Progress 40%
Progress 48%
Progress 54%
Progress 59%
Progress 65%
Progress 71%
Progress 76%
Progress 85%
Progress 90%
Progress 96%
Progress 100%
Progress 100%

Извлечение полученных результатов.

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 секунд моделирования. Верхний график показывает паттерн дыхания, то есть генератор ЛЧМ сигнала. Средний график соответствует изменению ЧСС. Обратите внимание, что на низких частотах ЧСС колеблется почти синхронно с изменением объема легких; однако на более высоких частотах она начинает запаздывать относительно дыхания. Кроме того, амплитуда сигнала ЧСС уменьшается с ростом частоты, подчеркивая характер общей частотной характеристики фильтра. Нижний график показывает изменение артериального давления. Можно заметить, что по мере роста частоты дыхательно-индуцированные изменения АД (ABP) становятся больше. Это является следствием усиливающегося влияния прямых механических эффектов дыхания при увеличении частоты.

# Визуализация результатов
max_time = 200      # Время по оX, с

plot(
    # График объема (V)
    plot(time_model, V_normal;
         xlabel="Время, с", 
         ylabel="Объем, л",
         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="Давление, мм рт.ст.",
         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|, мин/л",
         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      # Коэффициентусиления передаточной функции СА-узла для блуждающего влияния
Ks = 9      # Коэффициент усиления передаточной функции СА-узла для симпатического влияния
Ap = 0.1    # Весовой коэффицент для блуждающей ветви
As = 4      # Весовой коэффицент для симпатической ветви
fp = 0.07   # Частота среза передаточной функции СА-узла для блуждающего влияния, Гц
fs = 0.015; # Частота среза передаточной функции СА-узла для симпатического влияния, Гц
Запуск модели
out_atropine = run_model("Blood_circulation", @__DIR__);    # Запуск модели
Building...
Progress 0%
Progress 5%
Progress 10%
Progress 15%
Progress 22%
Progress 27%
Progress 32%
Progress 37%
Progress 54%
Progress 60%
Progress 66%
Progress 72%
Progress 83%
Progress 89%
Progress 95%
Progress 100%
Progress 100%
V_atropine= out_atropine["V"].value;
HR_atropine = out_atropine["HR"].value;
ABP_atropine = out_atropine["ABP"].value;
Визуализация результатов моделирования первого сценария
# Визуализация результатов
max_time = 200      # Время по оX, с

plot(
    # График объема (V)
    plot(time_model, V_atropine;
         xlabel="Время, с", 
         ylabel="Объем, л",
         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="Давление, мм рт.ст.",
         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|, мин/л",
         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      # Коэффициентусиления передаточной функции СА-узла для блуждающего влияния
Ks = 1      # Коэффициент усиления передаточной функции СА-узла для симпатического влияния
Ap = 2.5    # Весовой коэффицент для блуждающей ветви
As = 0.1    # Весовой коэффицент для симпатической ветви
fp = 0.2    # Частота среза передаточной функции СА-узла для блуждающего влияния, Гц
fs = 0.015; # Частота среза передаточной функции СА-узла для симпатического влияния, Гц
Запуск модели
out_propranolol = run_model("Blood_circulation", @__DIR__);    # Запуск модели
Building...
Progress 0%
Progress 5%
Progress 11%
Progress 18%
Progress 23%
Progress 28%
Progress 34%
Progress 41%
Progress 50%
Progress 55%
Progress 61%
Progress 66%
Progress 71%
Progress 79%
Progress 86%
Progress 92%
Progress 98%
Progress 100%
Progress 100%
V_propranolol= out_propranolol["V"].value;
HR_propranolol = out_propranolol["HR"].value;
ABP_propranolol = out_propranolol["ABP"].value;
Визуализация результатов моделирования второго сценария
# Визуализация результатов
max_time = 200      # Время по оX, с

plot(
    # График объема (V)
    plot(time_model, V_propranolol;
         xlabel="Время, с", 
         ylabel="Объем, л",
         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="Давление, мм рт.ст.",
         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|, мин/л",
         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|, мин/л",
    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.03 Гц),

  • Ослабление отклика на высоких частотах (выше 0.1 Гц),

  • Увеличение фазового запаздывания (крутой наклон фазовой кривой).

Введение пропранолола (блокирующего симпатику) приводит к состоянию с доминированием парасимпатики:

  • Сильное ослабление отклика на низких частотах,

  • Незначительные изменения на частотах выше 0.05 Гц,

  • Малый фазовый сдвиг (0-0.3 Гц), что указывает на быстрый отклик ЧСС на дыхание.

Заключение

Данная модель кровообращения описывает краткосрочное влияние дыхания на регуляцию частоты сердечных сокращений и артериального давления. Она демонстрирует взаимодействие дыхательных процессов, вегетативной нервной системы, барорефлекса и прямых механических воздействий на формирование ЧСС и АД.

Анализ частотных характеристик модели позволяет наблюдать, как меняется динамика и общий вклада симпатического и парасимпатического отделов нервной системы в контроль ЧСС при различных физиологических состояниях и воздействиях различных лекарственных препаратов.

Список используемых источников

1 - Saul, J.P., R.D. Berger, P. Albrecht, S.P. Stein, M.H.Chen, and R.J. Cohen.Transferfunction analysis of the circulation: unique insights into cardiovascular regulation. Am. 1. Physiol. 261 (Heart Circ. Physiol. 30): HI231-HI245, 1991.