Моделирование электрокардиосигнала

В примере показано моделирование электрокардиосигнала на основе несимметричных гауссовых функций, включая формирование нормального сигнала и ряда патологических изменений.

Подготовка рабочего пространства и настройка визуализации

Функция engee.clear() выполняет очистку рабочего пространств:

engee.clear()

Warning: detected a stack overflow; program state may be corrupted, so further execution might be unreliable.

Для визуализации результатов доступны два режима отрисовки графиков:

Для статической, быстрой отрисовки - используйте gr().
Для интерактивной, динамической визуализации с возможностью масштабирования и просмотра значений - используйте plotlyjs().

Выберите одну из команд, закомментировав вторую. Оба бэкенда являются взаимоисключающими, и активирован будет последний вызванный.

gr()
# plotlyjs()

Plots.GRBackend()

Электрокардиосигнал

На рисунке ниже представлены два последовательных сердечных цикла и основные элементы нормальной электрокардиограммы.

На данном рисунке показаны основные элементы нормальной ЭКГ и их место в пределах одного сердечного цикла. Зубец P отражает распространение возбуждения по предсердиям. После него следует сегмент PQ, который соответствует проведению импульса через атриовентрикулярное соединение, а интервал PQ охватывает весь путь от начала возбуждения предсердий до начала возбуждения желудочков. Комплекс QRS отражает деполяризацию желудочков: зубец Q соответствует начальному возбуждению межжелудочковой перегородки, зубец R - распространению возбуждения по основной массе миокарда желудочков, а зубец S - завершающему этапу деполяризации в базальных отделах желудочков. После комплекса QRS располагается сегмент ST, который соответствует периоду полного охвата миокарда желудочков возбуждением и в норме располагается близко к изолинии. Далее следует зубец T, отражающий реполяризацию желудочков. Интервал QT охватывает участок от начала комплекса QRS до конца зубца T и характеризует суммарное время деполяризации и реполяризации желудочков, включающую деполяризацию и последующую реполяризацию.

Такое представление ЭКГ удобно для дальнейшего математического моделирования, поскольку позволяет рассматривать кардиоцикл как последовательность отдельных информативных фрагментов, каждый из которых может быть описан и параметризован независимо.

Описание модели

В данном примере рассматривается моделирование электрокардиосигнала на основе суммы несимметричных гауссовых функций. Такой подход позволяет представить один кардиоцикл как совокупность отдельных фрагментов. В используемой реализации сигнал формируется из компонент , , , , и причём две компоненты и используются для более гибкого описания комплекса . В основе модели лежит представление сигнала в виде суммы отдельных волн:

где параметр определяет амплитуду -й компоненты, параметр задаёт положение её экстремума по времени, а функция отвечает за ширину волны слева и справа от точки экстремума. Для того чтобы одна и та же функция могла описывать как симметричные, так и несимметричные фрагменты сигнала, ширина задаётся кусочно:

Если , соответствующая волна является симметричной, а если , то она становится несимметричной. Именно это свойство позволяет гибко задавать форму отдельных участков ЭКГ и получать более реалистичную морфологию сигнала.

Реализация данной модели строится на двух функциях. Первая функция формирует одну несимметричную гауссову компоненту сигнала. Если обозначить такую компоненту через , то она записывается в виде:

В коде эта идея реализуется функцией asym_gauss, которая для каждой точки временной оси проверяет, находится ли она левее или правее точки , и в зависимости от этого использует либо левую ширину , либо правую ширину . Таким образом, одна и та же функция может описывать как положительные, так и отрицательные волны, а также как симметричные, так и асимметричные фрагменты электрокардиосигнала.

Вторая функция предназначена ecg_complex для формирования всего кардиоцикла. Она суммирует несколько компонент, каждая из которых вычисляется той же самой функцией asym_gauss, но уже со своим набором параметров. В развёрнутом виде эта формула имеет вид:

function asym_gauss(t, A, mu, b1, b2)
    y = zeros(length(t))

    for i in 1:length(t)
        if t[i] <= mu
            y[i] = A * exp(-((t[i] - mu)^2) / (2 * b1^2))
        else
            y[i] = A * exp(-((t[i] - mu)^2) / (2 * b2^2))
        end
    end

    return y
end

function ecg_complex(t, A, mu, b1, b2)
    y = zeros(length(t))

    for k in 1:length(A)
        y += asym_gauss(t, A[k], mu[k], b1[k], b2[k])
    end

    return y
end

ecg_complex (generic function with 1 method)

Моделирование нормального электрокардиосигнала

Для моделирования сигнала зададим частоту дискретизации 1 кГц. При таком значении шаг по времени составляет 1 мс. Далее формируется временная ось одного кардиоцикла длительностью 1 секунда.

fs = 1e3
Ts = 1/fs
t = 0:Ts:1-Ts

0.0:0.001:0.999

В данной части примера задаются параметры нормального электрокардиосигнала и выполняется его формирование. Массивы , , и определяют амплитуды компонент, положения их экстремумов по времени, а также ширины слева и справа, задающие форму и асимметрию отдельных волн. После этого с помощью функции ecg_complex формируется один нормальный кардиоцикл, который затем отображается на графике. Полученный сигнал используется как базовый вариант для последующего сравнения с патологическими формами ЭКС.

A_norm  = [ 0.11, -0.11, 0.95, 0.02, -0.18, 0.00, 0.20 ]
mu_norm = [ 0.18,  0.476, 0.500, 0.510, 0.523, 0.60, 0.74 ]
b1_norm = [ 0.03,  0.010, 0.010, 0.006, 0.012, 0.040, 0.045 ]
b2_norm = [ 0.05,  0.010, 0.010, 0.007, 0.014, 0.040, 0.065 ]

y_norm = ecg_complex(t, A_norm, mu_norm, b1_norm, b2_norm)

p1 = plot(t, y_norm,
    xlabel = "Время, с",
    ylabel = "Амплитуда",
    title = "Нормальный цикл",
    label = "Норма"
)

После формирования нормального сигнала модель используется для имитации патологических изменений ЭКС. В данном примере рассматриваются патологический зубец , плоский, отрицательный, высокий и асимметричный зубец , а также депрессия и элевация сегмента . Каждый вариант получается изменением параметров соответствующих компонент без изменения общей структуры модели.

Моделирование патологий электрокардиосигнала

В данной части рассматривается моделирование патологических изменений электрокардиосигнала на основе параметрической модели кардиоцикла. Изменяя амплитуды, временные положения и ширины отдельных компонент сигнала, можно воспроизводить характерные отклонения формы ЭКГ без изменения общей структуры модели.

В рамках примера рассмотрены следующие варианты патологических изменений:

патологический зубец Q;
плоский зубец T;
отрицательный зубец T;
высокий зубец T;
асимметричный зубец T;
депрессия сегмента ST;
элевация сегмента ST.

Патологический зубец Q

A_q  = [ 0.11, -0.32, 0.82, 0.03, -0.16, 0.00, 0.22 ]
mu_q = [ 0.18,  0.468, 0.500, 0.512, 0.532, 0.60, 0.74 ]
b1_q = [ 0.03,  0.016, 0.010, 0.006, 0.012, 0.040, 0.045 ]
b2_q = [ 0.05,  0.018, 0.010, 0.007, 0.014, 0.040, 0.070 ]

y_q = ecg_complex(t, A_q, mu_q, b1_q, b2_q)

p2 = plot(t, y_q,
    xlabel = "Время, с",
    ylabel = "Амплитуда",
    title = "Патологический зубец Q",
    label = "Патология Q"
)

Плоский зубец T

A_flatT  = [ 0.11, -0.10, 0.95, 0.03, -0.18, 0.00, 0.06 ]
mu_flatT = [ 0.18,  0.470, 0.500, 0.515, 0.535, 0.60, 0.75 ]
b1_flatT = [ 0.03,  0.010, 0.010, 0.006, 0.012, 0.040, 0.040 ]
b2_flatT = [ 0.05,  0.010, 0.010, 0.007, 0.014, 0.040, 0.060 ]

y_flatT = ecg_complex(t, A_flatT, mu_flatT, b1_flatT, b2_flatT)

p3 = plot(t, y_flatT,
    xlabel = "Время, с",
    ylabel = "Амплитуда",
    title = "Плоский зубец T",
    label = "Плоский T"
)

Отрицательный зубец T

A_negT  = [ 0.11, -0.10, 0.95, 0.03, -0.18, 0.00, -0.18 ]
mu_negT = [ 0.18,  0.470, 0.500, 0.515, 0.535, 0.60, 0.75 ]
b1_negT = [ 0.03,  0.010, 0.010, 0.006, 0.012, 0.040, 0.042 ]
b2_negT = [ 0.05,  0.010, 0.010, 0.007, 0.014, 0.040, 0.065 ]

y_negT = ecg_complex(t, A_negT, mu_negT, b1_negT, b2_negT)

p4 = plot(t, y_negT,
    xlabel = "Время, с",
    ylabel = "Амплитуда",
    title = "Отрицательный зубец T",
    label = "Отрицательный T"
)

Высокий зубец T

A_highT  = [ 0.11, -0.10, 0.95, 0.03, -0.18, 0.00, 0.42 ]
mu_highT = [ 0.18,  0.470, 0.500, 0.515, 0.535, 0.60, 0.75 ]
b1_highT = [ 0.03,  0.010, 0.010, 0.006, 0.012, 0.040, 0.040 ]
b2_highT = [ 0.05,  0.010, 0.010, 0.007, 0.014, 0.040, 0.060 ]

y_highT = ecg_complex(t, A_highT, mu_highT, b1_highT, b2_highT)

p5 = plot(t, y_highT,
    xlabel = "Время, с",
    ylabel = "Амплитуда",
    title = "Высокий зубец T",
    label = "Высокий T"
)

Асимметричный зубец T

A_asymT  = [ 0.11, -0.10, 0.95, 0.03, -0.18, 0.00, 0.24 ]
mu_asymT = [ 0.18,  0.470, 0.500, 0.515, 0.535, 0.60, 0.74 ]
b1_asymT = [ 0.03,  0.010, 0.010, 0.006, 0.012, 0.040, 0.028 ]
b2_asymT = [ 0.05,  0.010, 0.010, 0.007, 0.014, 0.040, 0.085 ]

y_asymT = ecg_complex(t, A_asymT, mu_asymT, b1_asymT, b2_asymT)

p6 = plot(t, y_asymT,
    xlabel = "Время, с",
    ylabel = "Амплитуда",
    title = "Асимметричный зубец T",
    label = "Асимметричный T"
)

Депрессия сегмента ST

A_depST  = [ 0.11, -0.10, 0.95, 0.03, -0.22, -0.07, 0.18 ]
mu_depST = [ 0.18,  0.47, 0.50, 0.515, 0.535, 0.62, 0.76 ]
b1_depST = [ 0.03,  0.010, 0.010, 0.006, 0.012, 0.055, 0.045 ]
b2_depST = [ 0.05,  0.010, 0.010, 0.007, 0.014, 0.080, 0.070 ]

y_depST = ecg_complex(t, A_depST, mu_depST, b1_depST, b2_depST)

p7 = plot(t, y_depST,
    xlabel = "Время, с",
    ylabel = "Амплитуда",
    title = "Депрессия сегмента ST",
    label = "Депрессия ST"
)

Элевация сегмента ST

A_eleST  = [ 0.11, -0.10, 0.95, 0.03, -0.16, 0.10, 0.20 ]
mu_eleST = [ 0.18,  0.47, 0.50, 0.515, 0.535, 0.62, 0.75 ]
b1_eleST = [ 0.03,  0.010, 0.010, 0.006, 0.012, 0.055, 0.045 ]
b2_eleST = [ 0.05,  0.010, 0.010, 0.007, 0.014, 0.090, 0.070 ]

y_eleST = ecg_complex(t, A_eleST, mu_eleST, b1_eleST, b2_eleST)

p8 = plot(t, y_eleST,
    xlabel = "Время, с",
    ylabel = "Амплитуда",
    title = "Элевация сегмента ST",
    label = "Элевация ST"
)

Сравнение нормального и патологических ЭКС

В данной части выполняется сравнение нормального и патологических вариантов ЭКС. Для наглядности результаты представлены в двух вариантах: в первом случае все сигналы отображаются на одном общем графике, что позволяет оценить различия формы и взаимное смещение отдельных участков; во втором случае каждый сигнал выводится на отдельной области построения с использованием компоновки layout = (4, 2), что делает более удобным детальное рассмотрение каждой патологии.

plot(t, y_norm, 
    xlabel="Время, с", 
    ylabel="Амплитуда", 
    title="Сравнение ЭКГ-сигналов", 
    label="Норма",
    legend_position = :outertopright
)

plot!(t, y_q, label="Патологический Q")
plot!(t, y_flatT, label="Плоский T")
plot!(t, y_negT, label="Отрицательный T")
plot!(t, y_highT, label="Высокий T")
plot!(t, y_asymT, label="Асимметричный T")
plot!(t, y_depST, label="Депрессия ST")
plot!(t, y_eleST, label="Элевация ST")

Отобразим все сигналы раздельно, используя layout = (4, 2).

plot(    p1, p2, p3, p4, p5, p6, p7, p8,
    layout = (4, 2),
    size = (1200, 1400)
)

Расщепленный зубец R

Ранее при описании модели были введены две компоненты зубца - и . Однако в ранее рассмотренных примерах такой случай отдельно не анализировался. В связи с этим далее рассмотрим его отдельно и выполним моделирование электрокардиосигнала с расщеплённым зубцом .

A_splitR  = [ 0.11, -0.10, 0.70, 0.62, -0.16, 0.00, 0.22 ]
mu_splitR = [ 0.18,  0.47, 0.495, 0.520, 0.538, 0.61, 0.74 ]
b1_splitR = [ 0.03,  0.010, 0.008, 0.008, 0.012, 0.045, 0.045 ]
b2_splitR = [ 0.05,  0.010, 0.009, 0.009, 0.014, 0.045, 0.070 ]

y_splitR = ecg_complex(t, A_splitR, mu_splitR, b1_splitR, b2_splitR)

p9 = plot(t, y_splitR,
    xlabel = "Время, с",
    ylabel = "Амплитуда",
    title = "Расщеплённый зубец R",
    label = "Расщеплённый R"
)

Заключение

В данном примере было рассмотрено параметрическое моделирование электрокардиосигнала на основе суммы несимметричных гауссовых функций. Показано, что такой подход позволяет формировать один цикл электрокардиосигнала в виде набора отдельных информативных компонент и гибко управлять его формой за счёт изменения амплитуд, временных положений и ширин соответствующих волн.