КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ В СРЕДЕ MAPLE ИСКУССТВЕННЫХ ЭКГ С НАЛИЧИЕМ ПАМЯТИ

АННОТАЦИЯ

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

ABSTRACT

The paper proposes a software implementation of a numerical algorithm for solving the McSherry fractional mathematical model in the Maple environment. The paper presents a statement of the problem and a numerical algorithm for its solution. Next, the program code in the Maple environment is given, which allows visualizing the simulation results, plotting oscillograms and phase trajectories.

Ключевые слова: математическая модель, Maple, дробная производная, численный алгоритм.

Keywords: mathematical model, Maple, fractional derivative, numerical algorithm.

Введение

Электрокардиограмма (ЭКГ) представляет собой изменяющийся во времени сигнал, отражающий протекание ионного тока, который вызывает сокращение и последующее расслабление сердечных волокон. Поверхностную ЭКГ получают путем регистрации разности потенциалов между двумя электродами, помещенными на поверхность кожи. Один нормальный цикл ЭКГ представляет собой последовательную деполяризацию/реполяризацию предсердий и деполяризацию/реполяризацию желудочков, которые происходят при каждом сердечном сокращении. Их можно приблизительно связать с пиками и впадинами кривой ЭКГ, помеченными P, Q, R, S и T [11], как показано на рис. 1.

Рисунок 1. Морфология среднего PQRST-комплекса ЭКГ, записанной у здорового человека

Дробная динамическая модель и методика ее исследования

В работе [1] была предложена математическая модель для построения искусственных ЭКГ здорового человека, реализующих морфологию PQRST (рис.1). В работе [2] был проведено продолжение исследований модели на качественном уровне. Исследованы хаотические режимы с помощью максимальных показателей Ляпунова. В работе [3] была предложена модель, которая обобщает модель из статьи [1] в случае учета эффектов наследственности с помощью производных дробных порядков, которые изучаются в монографиях [4]-[7]. Математическую модель, в которой присутствуют дробные производные будем называть дробной математической моделью

Дробную математическую модель для построения искусственной ЭКГ здорового человека можно записать в виде следующей задачи Коши [3]:

                                                             (1)

когда , ,  (четырехквадрантный арктангенс действительных частей элементов x и y,  и  - угловая скорость траектории при ее движении по предельному циклу,  - время процесса. Блуждание базовой линии было введено путем связывания базового значения  в (1) с частотой дыхания f с использованием

                                                                                                             (2)

когда A=0.15мВ.

Здесь производные дробных порядков понимаются в смысле Герасимова-Капуто, где 0 < α,β,γ < 1.

,

,                                                                                                 (3)

.

Дробную математическую модель (1) будем называть дробной математической моделью Макшерри.

Заметим, что в частном случае α = β = γ = 1 дробная математическая модель Макшерри (1) переходит в обычную математическую модель Макшерри из работы [1].

Для решения задачи (1) воспользуемся теорией конечно-разностных схем. Для этого будем считать, что функции x(t),y(t),z(t) обладают необходимыми условиями гладкости. Рассмотрим равномерную сетку по временной координате. Разобьём отрезок [0,T] на N равных частей с шагом дискретизации τ = T/N. Тогда функции решения x(t),y(t),z(t) перейдут в сеточные функции x(t_k),y(t_k),z(t_k), где t_k = kτ, k = 1,..,N. Аппроксимации производных дробных (3) при n=1 даются следующими формулами:

                                                                                  (4)

где ,

.

С учетом аппроксимаций (4), приведенных выше, мы можем записать дискретный аналог задачи (1).

                                           (5)

Мы получили нелокальную явную конечно-разностную схему, которая обладает первым порядком точности и условно сходится. В этой работе мы не будем проводить исследования вопросов устойчивости и сходимости схемы (5). Скажем, что всегда можно выбрать шаг расчетной сетки достаточно малым, чтобы выполнялось условия сходимости и устойчивости. Предложенный численный алгоритм (5) был реализован на языке среды символьной математики Maple.

Код программы в Maple

Maple позволяет достаточно просто реализовать численные алгоритмы, подобные алгоритму (5) в виде циклов, а также позволяет выводить результаты моделирования на соответствующих графиках [8]. В программном коде были выбраны значения параметров модели из работы [2]. Рассмотрим код более детально.

>restart;

>with(MTM); # библиотека специальных функций

>with(plots); # библиотека визуализации результатов моделирования

>T:= 50; # время моделирования

>N:= 2000; # количество узлов расчетной сетки

>tau:= evalf(T/N); # шаг расчетной сетки

>omega:= 1; # угловая скорость

>a[1] := 60; a[2] := -250; a[3] := 1500; a[4] := -375; # параметры из [2]

>a[5] := 37.5; b[1] := 0.25; b[2] := 0.1; b[3] := 0.1;

>b[4] := 0.1; b[5] := 0.4;

>Theta[1] := evalf(-1/3*Pi); Theta[2] := evalf(-1/12*Pi);

>Theta[3] := 0; Theta[4] := evalf(1/12*Pi);

>Theta[5] := evalf(1/2*Pi); delta := 0.15; f := 0.25;

>alpha1 := 0.9; alpha2 := 0.8; alpha3 := 0.7; # значения порядков дробных производных

>A[1] := tau^(-alpha1)/GAMMA(2 - alpha1);

>A[2] := tau^(-alpha2)/GAMMA(2 - alpha2);

>A[3] := tau^(-alpha3)/GAMMA(2 - alpha3);

>xx[0] := 0.1;# начальные условия

>yy[0] := 0.1;

>zz[0] := 0.25;

>mp := (l, m) -> l - floor(l/m)*m;

>atan2 := (yyy, xxx) -> piecewise(0 < xxx, arctan(yyy/xxx), And(xxx < 0, 0 <= yyy), arctan(yyy/xxx) + Pi, And(yyy < 0, xxx < 0), arctan(yyy/xxx) - Pi, And(xxx = 0, 0 < yyy), 1/2*Pi, And(xxx = 0, yyy < 0), -1/2*Pi);

>xx[1] := ((1 - evalf(sqrt(xx[0]^2 + yy[0]^2)))*xx[0] - omega*yy[0] + A[1]*xx[0])/A[1];

>yy[1] := ((1 - evalf(sqrt(xx[0]^2 + yy[0]^2)))*yy[0] + omega*xx[0] + A[2]*yy[0])/A[2];

>R[0] := add(a[j]*mp(atan2(yy[0], xx[0]) - Theta[j], 2*evalf(Pi))*exp(-mp(atan2(yy[0], xx[0]) - Theta[j], 2*evalf(Pi))*mp(atan2(yy[0], xx[0]) - Theta[j], 2*evalf(Pi))/((2*b[j])*b[j])), j = 1 .. 5);

>zz[1] := (A[3]*zz[0] + R[0] - zz[0])/A[3];

> for k to N - 1 do

xx[k + 1] := ((1 - evalf(sqrt(xx[k]^2 + yy[k]^2)))*xx[k] - omega*yy[k] - A[1]*add(((n + 1)^(1 - alpha1) - n^(1 - alpha1))*(xx[k - n + 1] - xx[k - n]), n = 1 .. k - 1) + A[1]*xx[k])/A[1];

 yy[k + 1] := ((1 - evalf(sqrt(xx[k]^2 + yy[k]^2)))*yy[k] + omega*xx[k] - A[2]*add(((n + 1)^(1 - alpha2) - n^(1 - alpha2))*(yy[k - n + 1] - yy[k - n]), n = 1 .. k - 1) + A[2]*yy[k])/A[2];

R[k]:= add(a[j]*mp(atan2(yy[k], xx[k]) - Theta[j], 2*evalf(Pi))*exp(-mp(atan2(yy[k], xx[k]) - Theta[j], 2*evalf(Pi))*mp(atan2(yy[k], xx[k]) - Theta[j], 2*evalf(Pi))/((2*b[j])*b[j])), j = 1 .. 5);

zz[k + 1] := ((R[k] + (-zz[k] + delta*sin(2*evalf(Pi)*f*k*tau))) - A[3]*add(((n + 1)^(1 - alpha3) - n^(1 - alpha3))*(zz[k - n + 1] - zz[k - n]), n = 1 .. k - 1) + A[3]*zz[k])/A[3];

end do;

>X := seq([tau*j, xx[j]], j = 0 .. N);

>Y := seq([tau*j, yy[j]], j = 0 .. N);

>Z := seq([tau*j, zz[j]], j = 0 .. N);

>XYZ := seq([xx[j], yy[j], zz[j]], j = 0 .. N);

>pointplot3d([XYZ], style = line, labels = ["x(t)", "y(t)", "z(t)"], labelfont = ["HELVETICA", "ROMAN", "BOLD", 12], axesfont = ["HELVETICA", "ROMAN", "BOLD", 10]);

Рисунок 2. Фазовая траектория

>pointplot([X], style = line, labels = ["t", "x(t)"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", "ROMAN", "BOLD", 12], axesfont = ["HELVETICA", "ROMAN", "BOLD", 10]);

Рисунок 3. Осциллограмма x(t)

>pointplot([Y], style = line, labels = ["t", "y(t)"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", "ROMAN", "BOLD", 12], axesfont = ["HELVETICA", "ROMAN", "BOLD", 10]);

Рисунок 4. Осциллограмма y(t)

>pointplot([Z], style = line, labels = ["t", "z(t)"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", "ROMAN", "BOLD", 12], axesfont = ["HELVETICA", "ROMAN", "BOLD", 10]);

Рисунок 5. Осциллограмма z(t)

На рис. 2-5 приведены результаты моделирования по численной схеме (5). Мы видим, что фазовая траектория на рис.5 сохраняет морфологию PQRST как на рис. 1. Порядки дробных производных можно рассматривать как дополнительную параметризацию сигнала ЭКГ для более гибкого моделирования.

Заключение

В работе приведен код компьютерной программы на языке Maple, который дает возможность проводить визуализацию результатов построения искусственной ЭКГ здорового человека. Показано, что результаты моделирования сохраняют морфологию PQRST и поэтому могут быть использованы в медицинской диагностике. Результаты работы могут быть развиты по нескольким направлениям. Первое направление - это уточнение параметров модели (1) по экспериментальным данным, другое направление, связанное с развитием модели (1), например, с учетом переменных порядков дробной производной от времени [9], еще одно направление связано с изучением хаотических и регулярных режимов [10].

Работа выполнена в рамках гранта президента РФ "Развитие математических моделей дробной динамики с целью исследования колебательных процессов и процессов с насыщением № МД758.2022.1.1.

Список литературы:

1. McSharry P.E., Clifford G.D., Tarassenko L., Smith L.A. A dynamical model for generating synthetic electrocardiogram signals, IEEE transactions on biomedical engineering., 2003. vol. 50, no. 3, pp. 289-294.
2. Марценюк В.П., Сарабун Р.О. Исследование нелинейной динамики в модели МакШерри на основе экспонент Ляпунова // Вестник Воронежского государственного университета. Серия: Системный анализ и информационные технологии. – 2014. – №. 2. – С. 57-61.
3. Алимов Х.Т. Дробная математическая модель Макшерри / Х.Т. Алимов, Ф.Х. Дзамихова, Р.И. Паровик // Вестник КРАУНЦ. Физико-математические науки. – 2023. – Т. 42, № 1. – С. 164-179. – DOI 10.26117/2079-6641-2023-42-1-164-179. – EDN BSUNVQ.
4. Oldham K., Spanier J. The Fractional Calculus. Theory and Applications of Differentiation and Integration to Arbitrary Order. London: Academic Press, 1974. 240 pp.
5. Miller K., Ross B. An Introduction to the Fractional Calculus and Fractional Differntial Equations. New York: A Wiley-Interscience Publication, 1993. 384 pp.
6. Нахушев А.М. Дробное исчисление и его применение. Москва: Физматлит, 2003. 272 с.
7. Kilbas A.A., Srivastava H. M., Trujillo J. J. Theory and Applications of Fractional Differential Equations, vol. 204: Amsterdam, 2006. 523 pp.
8. Дьяконов В. Maple 9.5/10 в математике, физике и образовании. – Litres, 2022.
9. Ortigueira M.D., Valerio D., Machado J. T. Variable order fractional systems, Communications in Nonlinear Science and Numerical Simulation, 2019. vol. 71, pp. 231–243.
10. Parovik R.I. Studies of the Fractional Selkov Dynamical System for Describing the Self-Oscillatory Regime of Microseisms / R. I. Parovik // . – 2022. – Vol. 10, No. 22. – P. 4208. – DOI 10.3390/math10224208. – EDN MJIHSH.