Алгоритм и программа идентификации параметров модели кровообращения Текст научной статьи по специальности «Математика»

УДК 612.172.6;577.4

АЛГОРИТМ И ПРОГРАММА ИДЕНТИФИКАЦИИ ПАРАМЕТРОВ МОДЕЛИ

КРОВООБРАЩЕНИЯ

Ю.Б. Нечаев, Б.Н. Воронков, И.С. Арзамасцев

Математическое моделирование сердечнососудистой системы и программная реализация оценки полученных параметров

Ключевые слова: модель, кровообращение, давление, исследование

В настоящее время инвалидизация и летальность при сердечнососудистых заболеваниях занимают одно из первых мест в шкале негативных факторов, влияющих на состояние общества, особенно в экономически развитых странах. Возможно, уровень этих негативных факторов был бы значительно ниже, если бы были известны закономерности процессов, которыми определяются тяжесть осложнений и течение болезней. Немаловажную роль в выявлении таких закономерностей играет математическое моделирование. Уже более полувека ведутся работы по исследованию сердечнососудистой системы и математическому моделированию кровообращения. Особенно актуальными оказались возможности математических моделей в кардиохирургической практике и их применение непосредственно в ходе экспериментальных исследований [1, 3]. Медицина наших дней характеризуется возрастанием вмешательства врача в глубоко интимные механизмы жизнедеятельности организма. В связи с этим возникает необходимость обеспечения наилучшего лечения для каждого отдельного (а не среднестатистического) больного. Эта задача может быть решена методами оперативной оценки параметров кровообращения, связанных в систему математической моделью.

Сердечнососудистая система рассматривается как совокупность п подсистем, описываемых в сосредоточенных параметрах. Относительно средних давлений в подсистемах описание модели имеет вид

[4]:

Р (г) = БЯТ Р(г) + Боа); (1)

У(1) = СР(1),

где Р(г) - п-мерный вектор состояния (давлений), Р(г) = [Р1 (г)]; 0(г) - п-мерный вектор тестовых воздействий, 0(г) = Ю^)], (0;(г) - объемная скорость восполнения крови в !-й элемент, например, из шприца); У(г) - т-мерный вектор наблюдений, У (г) = [У1 (г)]; Б, Я - матрицы п х п свойств элементов модели: Б = diag [е1 ,.. .,еп] (е1 - жесткость 1-го элемента), Я = [Р^] (Р^ - проводимость потока крови между элементами ) и 1); С - матрица т х п наблюдений с известными элементами, С = [С^].

В случае двухэлементной модели сердечнососудистая система рассматривается в виде двух упругих

Нечаев Юрий Борисович - ВГУ, д-р физ.-мат. наук, профессор, тел. 8 (960) 112-91-30

Воронков Борис Николаевич - ВГУ, канд. техн. наук, доцент, тел. 8 (473) 265-41-28

Арзамасцев Иван Сергеевич - ООО «ИГиТ», инженер, тел. 8 (950) 772-84-73

резервуаров - артериального (индекс 1) и венозного (индекс 2). Рассмотрим клиническую ситуацию, когда непрерывный контроль возможен только для артериального давления (венозное контролируется до воздействий). Тестовое воздействие (восполнение крови) приложено к артериальному резервуару (02 (г) = 0). Тогда полагаем в (1) п = 2 и задаем матрицы свойств подсистем Я и Б и матрицу наблюдений С в следующем виде:

Б = diag [еь е2];

Я = [Ы (Р12 = -Р11 = р);

р21 = р22 = (а + р)

где: р - общая периферическая проводимость; а - характеристика функции сердца. С = [1, 0].

Отметим, что а, р, е1, е2 являются обобщенными характеристиками сердечно-сосудистой системы. Так, а определяет результативную насосную функцию, отражая влияние частоты, фазовых соотношений, сократимости и диастолической активности сердца; р - общая проводимость сосудистого ложа, не отражающая соотношений сопротивлений различных тканей и органов. Аналогичное содержание обобщенных жесткостей артериальной и венозной систем имеют параметры моделей е1 и е2 .

Уравнение, описывающее реакцию артериального давления Р1 на тестовое воздействие 01, можем записать в виде

Р1 (г) + а1Р1© = Ьпу п© + Ьшуп © (2)

где: а1 = р е1 + р е2 + а е2 ;

Ьц = е1 ;

Ъю = (а + р) е1 е2;

УпО) =1 01(т) dт.

Определение параметров модели а, р, е1, е2 проводится в два этапа. На первом определяются коэффициенты уравнения (2) по контролю Р1(г) и 01(г) на интервале [0, Т], а на втором этапе рассчитываются собственно параметры.

Введем опорное уравнение со значениями коэффициентов а1*, Ъ11*, Ъ10*, являющихся априорными оценками а1, Ъ11, Ъ10:

X(г) = - а1* х(г) + Ъю* У1(г) + Ъ„* У 1© (3)

Различие решений уравнений (3) и (2) при одинаковых начальных условиях Р()(0) = х ()(0), 1 =

0,.1, обусловлено отличием соответствующих ко-

эффициентов. Предположим далее, что существуют функции N^1), М0(г), М1(г), которые при любых конечных отклонениях аь Ъп, Ъ10 от их опорных значений а1*, Ъ11*, Ъ10* приводят для г > 0 к выполнению неравенства

Р1О) = х(г) + Аа1 Nl(t) + аЪц м^) + +

АЪю М0(1), (4)

где: А а1= а1 - аД А Ъп= Ъп - Ъп*, А Ъю= Ъю -

Ъ10*.

Найдем выражения, из которых можно было бы определить N(1), М0(1), М1(1). Для этого, подставив Р1(1) и Р1(1) из (4) в (2), получим:

X (1) + Аа1 N 1(1) + АЪц М 1(1) +

+ АЪю М „(Ч) + (Аа1 + а1*) ( х (1) +

+ Аа1 N1(1) + АЪц М1(1) + АЪю М^(1)) =

(Ъц*+ АЪП ) у 11(1)+ (Ъю*+ А Ъю ) У11 (1). (5)

Выполняя преобразования с учетом уравнений (2) и (3) получим:

А а1 [ Р1© + N 1(1) + а1* N1(1) ] + А Ъп [М 1(1) +

+ а1* М1(1) -у п0)] +АЪ10[М 0(1) +

+ а1*М0(1)-уп(1)]=0 (6)

Поскольку мы предположили, что (4), а следовательно и (6) выполняются при любых конечных Аа1, АЪ11, АЪ10, то все выражения в квадратных скобках уравнения (6) должны быть равны нулю, т. е.

N 1(1) = - а1* N1(1) - Р1(1); "I

М 0 (1) = - а1*М0© + У11©; и (7)

М 1(1) = - а1* М1(1) + у пф J Система (7) является исходной для определения М0(1), М1(1), N1(t), называемых функциями чувствительности второго рода [6]. Начальные условия для М0(1), М1(1), N^1), нулевые, что следует из (4) и равенства начальных условий для Р1(1) и Х(1).

Чтобы получить несмещенные, с минимальной дисперсией оценки коэффициентов уравнения (2) по измеряемым с погрешностью величинам Р1(1) и 01(1) используем фильтр Калмана [5, 7]. Для этого представим (4) в дискретном виде:

г(к) = Н(к) Х(к) + v(k), (8)

где : к = 1к (по определению) ; 2(к) = Р1(к) - х(к);

Н (к) = [^(к),М1(к),М0(к)]; ЯТ = [А а1, А Ъп, АЪю]; v(k) - некоррелированная случайная последовательность с нулевым математическим ожиданием и дисперсией Я(к), характеризующей погрешность измерений. Поскольку параметры модели считаем неизменными на интервале, соизмеримом со временем затухания переходных процессов в системе кровообращения, то

Х(к+1) = Х(к) (9)

Уравнения (9) и (8) можем считать соответственно уравнениями состояния и наблюдения некоторой системы с вектором состояния X. Тогда для оценки X можем использовать алгоритм фильтра Калмана :

Л Л

Я (к+1) = Я (к) + W(k+1) [ г(к+1) - Н(к+1) Х(к) ]

Л

с начальным условием Я (0) = 0. Вектор коэффициентов усиления определяется из соотношений:

W(k+1)=S(k)HТ(k+1)[H(k+1)S(k)HТ(k+1)+R(k+1)]

8(к+1) = [I-W(k+1)H(k+1)]S(k)[I-W(k+1)H(k+1)]+ +W (к+1 )R(k+1)W (к+1),

где I = diag[1,...Д]. Начальное условие для матрицы ковариаций 8(0) = ^0.

По известной оценке вектора X на 1-м шаге определяются оценки коэффициентов уравнения (2) :

€1 ( 1 ) = а1* + А €1( 1 );

§п( 1 ) = Ъп* + А Ъгп( 1 );

§ю( 1 ) = Ъю* + А §ю( 1 ).

Требуется найти четыре параметра модели а, р, е1, е2, а имеется три независимых соотношения, связывающих их. Поэтому для однозначного определения параметров в качестве четвертого, недостающего соотношения используем связь (в рамках модели) артериального и венозного давления в установившихся режимах:

Р1 Р2-1 = (а + р ) р-1.

Тогда для расчета параметров получаем следующие формулы :

р = (а1 Ъп - Ъю) Ъп -2; а = (Р1 - Р2) Р2-1 р; е1 = Ъ11;

е2 = Ъ10 Ъ11 Р1 Р2р .

Процедура оценки продолжается до получения необходимой точности.

Погрешность идентификации (оценки параметров) характеризуем величиной критерия:

1 Ь

'■=Ь *

к=1

Р1(к) - Рт1(к)

где РТ1 - рассчитанная по модели кривая переходного процесса артериального давления, соответствующая оценкам параметров, полученным по Ь дискретным измерениям Р1 и 01 на интервале [ 0, Т].

На основе описанной математической модели разработана программа оценки параметров. Для проверки работоспособности программы необходимы данные, которые могут быть получены только в условиях медицинского эксперимента. Это значения артериального и венозного давлений до начала тестового воздействия, затем, контролируемые в период воздействия артериальное давление, объем и время воздействия. Эти данные должны фиксироваться с определенным шагом (например, 0,1с) в течение всего воздействия. Также необходимо знать априорные оценки параметров модели. Работоспособность программы проверялась с использованием данных [3], где оценивалось состояние сердечно-сосудистой системы собак в эксперименте. В качестве тестовых воздействий использовались кровопотеря и крово-восполнение в бедренную артерию или в нижнюю полую вену. Объем разового изменения циркулирующей крови 100-150 мл., скорость воздействия 8-12 мл./с. Во всех опытах параметры существенно индивидуальны. В работе программы использовались, описанные в [3] артериальное и венозное

давления, априорные оценки параметров модели, вид и общий объем тестового воздействия. Результаты работы программы сравнивались с результатами, приведенными в [3]. Оперативные данные по контролю артериального давления, объема крови и времени тестового воздействия в указанной работе не приводятся. Анализ литературных источников [2 - 4] показал, что во всех случаях кровопускание приводит к выраженному снижению артериального, венозного давлений и минутного объема крови, а восполнение к повышению. Поэтому в качестве промежуточных данных были взяты плавное повышение (или понижение) артериального давления.

Файл данных для работы программы должен иметь вид :

60 = п

60.00 = х

1.536 = аа1

5.405 = ЬЬ11

0.122 = ЬЬ10

0.80 = Р2

0.01 = 80

0.0001 = ЄрБІ

518.8 = УИ

1[І] Р[І] ОМ

1.00 60.00 0.00

2.00 60.00 0.00

3.00 60.00 0.00

4.00 60.00 0.00

5.00 60.00 0.00

5.20 60.25 2.00

5.40 60.40 4.00

5.60 60.65 6.00

14.80 71.52 98.00

15.00 71.61 100.00

16.00 71.60 0.00

17.00 71.40 0.00

18.00 71.20 0.00

19.00 71.05 0.00

20.00 71.05 0.00

Рис. 1. Результат работы программы

Файл результатов имеет вид:

Опорное уравнение имеет следующие коэфици-

енты:

х=60.00

аа1=1.5360

ЪЪ11=5.4050

ЪЪ10=0.1220

Венозное давление Р2=0.80

произведено 60измерений.

І : 1[І] Р[І] 0[І] : Y[i] :

1 : 1.00 60.00 0.00 -0.06

2 : 2.00 60.00 0.00 0.02 :

: 3 : 3.00 : 60.00 0.00 0.03 :

: 4 : 4.00 : 60.00 0.00 0.02 :

: 5 : 5.00 : 60.00 0.00 0.03 :

: 6 : 5.20 : 60.25 2.00 0.00 :

: 7 : 5.40 : 60.40 4.00 -0.11

: 8 : 5.60 : 60.65 6.00 0.08 :

: 9 : 5.80 : 60.90 8.00 0.02 :

59 : 19.00 : 71.05 : 0.00 : -0.05

60 : 20.00 : 71.05 : 0.00 : -0.07

Начальное условие матрицы ковариации 80=0.01 Дисперсия Я=0.0001

Оценки параметров

Характери- : Общая : Жесткость : Жесткость

стика : перефери-: артериаль- : венозного

: функции : ческая : ного : резервуара

: сердца : проводи- резервуара :

мость

: а : г1 : е1 : е2

: мл/мм.рт. : мл/мм. рт. : мл/мм. рт. : мл/мм. рт.

: ст. *с : ст. *с : ст. *с : ст. *с

==: 20.5113 0.0014 : : 0.2772 : 5.4393 :

шТ = 0.0294

Рис. 2. Зависимость оценок параметров от интервала наблюдения

Рассмотренный пример соответствует опыту 8, приведенному в [3].

Тестовое воздействие - восполнение 100 мл крови в артериальный резервуар. Результаты идентификации приведенные в [3] :

а=17,857 р=0,308 е1=4,739 е2=0,002.

Программа может быть использована не только как составная часть комплекса программ поддержки принятия решений врачом-кардиологом, но и при раз-

работке более детальных четырех, шести и четырнадцатиэлементных моделей гемодинамических систем.

Литература

1. Лищук В. А. Идентификация системы кровообращения / В. А. Лищук, В. Л. Столяр // Автоматика, 1979, №1. - С. 17 - 24.

2. Абрамян А. С. Математическая модель для исследования крововосполнения и кровопотри / А. С. Абрамян // Вопросы кибернетики : Регуляция и саморегуляция вегетативных функций. - М. : АН СССР, 1977. -С. 126 - 132.

3. Бураковский В. И. Анализ функции и состояния сердечно-сосудистой системы в эксперименте с помощью математической модели / В. И. Бураковский, В. А. Лищук, М. В. Соколов // Вестник АМН СССР, 1976. - № 10. - С. 57 - 68.

4. Лищук В. А. Математическая теория кровообращения / В. А. Лищук. - М. : Медицина, 1991. - 256 с.

5. Сейдж Э. П. Идентификация систем управления / Э. П. Сейдж, Дж. Л. Мелса. - М. : Наука, 1974. -248 с.

6. Сердюков В. А. Функции чувствительности второго рода / В. А. Сердюков, Л. Е. Сердюкова // Актуальные вопросы технической кибернетики. - М. : Наука, 1972. - С. 94 - 98.

7. Современные методы идентификации систем / Под ред. П. Эйкхоффа. - М. : Мир, 1983. - 400 с.

Воронежский государственный университет

Общество с ограниченной ответственностью «Инженерная геодезия и топография»

ALGORITHM AND PROGRAM OF IDENTIFYING THE PARAMETERS OF THE MODEL BLOOD CIRCULATION

Y.B. Nechaev, B.N. Voronkov, I.S. Arzamastsev

Mathematical modeling of the cardiovascular system and software implementation of the evaluation of the parameters

Key words: model, blood flow, pressure, research