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

9

БИОМЕДИЦИНСКИЕ ТЕХНОЛОГИИ

УДК 51.76

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЛИЯНИЯ УСТАНОВКИ КАВА-ФИЛЬТРА НА ГЕМОДИНАМИКУ КРОВЕНОСНОЙ СИСТЕМЫ

Ю.А. Иванов, Т.К. Добросердова

Моделирование кровеносной системы является обширной областью для научной деятельности и включает в себя множество аспектов. В этом смысле вопрос о строении сосуда и его реакции на деформацию логично рассматривать отдельно. На основе биологических особенностях строения ткани реализуется волоконная модель эластичной стенки сосуда, идея которой была сформулирована в работах Чарльза Пескина [6]. Подробно рассматривается калибровка модели с помощью специальных входных данных - уравнений состояния сосуда. Отдельная часть работы посвящена интеграции волоконной модели эластичности и глобальной модели кровообращения. Взаимодействие моделей дает возможность исследования влияния установки модельного кава-фильтра на гемодинамику.

Ключевые слова: тромбоэмболия, трансмуральное давление, уравнение состояния сосуда, кава-фильтр.

Введение

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

Модель глобального кровообращения

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

^ + = 0; (1)

дt дх

аГ и-1

^-Ц-^^и); (2)

дt дх

р = рс2 / (5), (3)

Сехр(55-1 -1) -1, 5 > 5

где / (5) = <! . Выражения (1), (2) отражают законы сохранения массы и количества

[ 1П(55"'), 5 < 5

движения соответственно. Функция у задает вязкое трение. Уравнение (3) учитывает эластичные свойства стенок сосудов. В настоящей работе будет использоваться гидродинамическая модель, основанная на системе (1)-(3) [2].

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

Рк (50, х)) - р. О) = е^^ 0, х)ик 0, х\

т

Xе А& х)ик(t, х) = с1,

к=1

где к - номер сосуда, р{ - давление в точке стыковки; е. = -1 и х = 0, если сосуд выходит из этой точки, е = 1 и х = I в противном случае (I - длина сосуда); К - сопротивление, которое преодолевает жидкость, протекающая через эту область; т - количество сосудов, имеющих данный узел своим концом, суммирование проводится по всему множеству. Замыкается система уравнением совместности - условием, накладываемым характеристическими кривыми, идущими внутрь области интегрирования. Поскольку размеры артериол и венул очень малы, их моделирование невозможно и не необходимо, поэтому в узлах стыковки артерий с венами используется такая же система граничных условий, как и в предыдущем случае, с той лишь разницей, что сопротивление К берется значительно больше, заменяя поколения сосудов, которыми приходится пренебречь. В настоящей работе используется динамическая модель четырехкамерного сердца [2]. Система уравнений, лежащая в ее основе, дополненная условиями равенства давлений и потоков в соответствующих предсердиях (желудочках) и во входящих (выходящих) в сердце венах (артериях), а также уравнениями совместности в концевых точках этих же сосудов, определяет систему граничных условий в узлах сопряжения сосудистых деревьев с сердцем.

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

Функция /(5) выведена эмпирически [3], поэтому никак не отражает воздействий на сосуд. Попытаемся получить новое уравнение состояния, используя волоконную модель эластичной стенки.

Модель эластичной стенки сосуда

Волоконная модель эластичной стенки сосуда основывается на биологических особенностях ее строения [6]. Она представляет собой систему переплетенных волокон, обладающих упругими свойствами.

Рис. 1. Волоконная модель эластичной стенки

Существует три типа волокон: эластические, гладкомышечные и коллагеновые. Модули Юнга первых двух примерно равны 3-105 Па и в данной модели будут объединены в один класс. Эластические волокна образуют достаточно густую сеть, особенно на внутренней оболочке стенки. Они создают эластическое напряжение, которое противодействует кровяному давлению, растягивающему сосуд. Коллагеновые волокна имеют модуль Юнга 1-108 Па. Сеть, образованная ими в средней и наружной оболочках, оказывает сопротивление растяжению большее, чем сеть эластических волокон [4]. Коллагеновые волокна располагаются относительно свободно в стенке сосуда, иногда даже образуя складки, так что противодействие внутреннему давлению оказывают лишь тогда, когда сосуд уже растянут до определенного состояния. Модельная реализация описанной выше структуры представляет собой набор нескольких пар спиралевидных (гелико-видных) волокон, закрученных вдоль стенки в противоположных направлениях (рис. 1). Каждая группа определяется своим набором параметров - модулем упругости, скоростью реакции на деформацию, степенью закрученности и т.д. Также в модель могут быть введены и другие типы волокон, например, кольцевидные, опоясывающие сосуд по всей длине, или прямолинейные (лучевидные).

Определим силы упругости, возникающие вследствие деформации волокна, заданного набором точек. Пусть 5 - переменная Лагранжа, определяющая расстояние вдоль волокна от данной точки до некоторой начальной; X (5, I) - вектор-функция, описывающая положение физической точки с лагранжевой

координатой 5; t - время; Т - напряжение деформированного волокна. Предположим, что напряжение подчиняется нелинейному обобщению закона Гука:

(4)

где Т* - коэффициент упругости волокна. Тогда локальная сила упругости / , действующая со стороны стенки, выражается следующим образом:

/ 4 И; ®

Г д X > _Т * Г д X Л

ст 1

V дя V дя /

-_ д X дя

дя

' "дХ

дя

(6)

где - - направление касательной в точке волокна.

Осуществим дискретизацию уравнений (4)-(6). Пусть Ль - количество расчетных точек на волокне; Хк - координата к-ой точки в трехмерном пространстве, где к = 1,...,ЛЬ; Д5 - расстояние между соседними точками волокна в недеформированном состоянии. Предполагаем, что Д5 будет одинаково для точек, принадлежащих однотипным волокнам. Далее в соответствии с формулами (4), (6) дискретизируем Т и -следующим образом:

Тк+х=^1 -Хк| /Д)_1Хк+Д- Хк1 -1; (7)

X - X

- - к+1_±±. (8)

-к+>2 _ |Хк+1 - Хк!

Дискретная сила упругости в к-ой точке определяется по формуле

_ Тк+>2-к+12- Тк -12-к - >2

/к _ Д5 • ( )

Реализуем модель эластичной стенки в виде двух наборов спиралевидных волокон с двумя различными коэффициентами упругости и параметрами, характеризующими степень закрученности волокон. Первый набор, моделирующий связки эластина, реагирует на любую деформацию, второй, моделирующий коллагеновые волокна - только если деформация превышает определенный порог. Тогда калибровка модели сводится к отысканию пяти параметров: Те, Тс - коэффициенты упругости двух наборов волокон; ДЬг, ДЬс - «шаги» спирали; Кс - параметр, определяющий степень деформации, до которой второй набор волокон (коллагена) не вносит вклада в реакцию упругой стенки сосуда.

Рассмотрим следующую задачу. Пусть в модель на вход подаются данные - уравнение состояния некоторого сосуда. Оно может быть задано как дискретный набор данных (полученных в результате эксперимента), так и в виде некоторой функции, описывающей уравнение состояния реального сосуда. Задача - подобрать коэффициенты так, чтобы уравнение состояние нашей модели было наиболее близко в смысле ¿2-нормы к целевой функции на заданном отрезке, т.е. на определенном диапазоне изменения сечения сосуда.Под уравнением состояния модели и под целевой функцией в случае дискретных входных данных будем понимать кусочно-линейную интерполяцию. Предположим, что у нас есть только один набор волокон, описываемый двумя параметрами Т и Д£. Требуется отыскать их значения, чтобы наилучшим образом удовлетворить целевой функции Р(5) на отрезке [50, 51], и где Р(50)=0. Функция Р(5) для реальных сосудов имеет экспоненциальный характер роста на данном промежутке. Заметим, что уравнение состояния модели имеет линейный характер при малых значениях параметра Д1 и усиливает выпуклость при его увеличении. Параметр Т будет играть роль нормирующего коэффициента. Основываясь на этих замечаниях и применяя метод золотого сечения, строим итерационный алгоритм, минимизирующий £2-норму ошибки на заданном отрезке. Отыскание совокупности описанных ранее пяти параметров модели - Те, Тс, Д£г, ДЬс и Яс - основывается на тех же принципах. Если рассматривать растяжение сосуда при 5<250, исследование используемого алгоритма калибровки модели для различных уравнений состояния показало, что норма ошибки составляет не более 4% от нормы приближаемой функции Р(5).

Влияние кава-фильтра на гемодинамику

Рассмотрим процесс установки кава-фильтра в сосуд. Предположим, что имплантант представляет собой два соосно-соединенных металлических конуса, совмещенных вершинами и образованных чередующимися лучами с крючками или фиксаторами на концах, которыми он прикрепляется к стенкам вены (см., например, фильтр «Зонтик» компании «Комед» [5]). При этом фильтр не фиксирует сечение сосуда в этой точке - оно может меняться в зависимости от кровяного давления. В предположении, что это давле-

ние компенсируется реакцией стенки вне зоны крепления фильтра, с помощью модели эластичной стенки получим новую зависимость р(5, §), где - расстояние от данной точки до места имплантации (рис. 2).

Р, Па

Рис. 2. График модифицированного уравнения состояния: p(S, Ç), где Ç = 0 (кривая 1); 0,5 (2); 1 (3); 1,5 (4); 2 (5); 2,5 (6); 3 (7); 3,5 (8); 4 (9); 4,5 (10); 5 (11).

Пусть кава-фильтр стоит в середине сосуда, длина которого равна 10 см, диаметр 3 см, скорость распространения малых возмущений 350 см/с. Тогда на его концах уравнение состояния практически совпадает с первоначальной зависимостью (3). Чем ближе к точке установки имплантанта, тем больше угол наклона графика: как при повышении, так и при понижении трансмурального давления сосуд меньше подвержен деформации из-за близкого расположения упругого кава-фильтра, препятствующего изменению формы. Кроме того, поскольку имплантант растягивает вену, при равенстве внутреннего и внешнего давлений на стенку при движении от ее концов к середине площадь поперечного сечения становится больше.

dpmepuu Вены

Рис. 3. Упрощенная схема кровеносной системы

Используем модифицированное уравнение состояния в глобальной модели кровообращения. Проведем вычислительный эксперимент на упрощенной схеме кровеносной системы (рис. 3). Венозная часть большого круга будет представлена тремя ветвями, артериальная - пятью (таблица). Артериальные ветви 4, 5 соединяются с венами 2, 3 соответственно. В артериях 4, 5 берется нарочито большое сопротивление, чтобы сохранить правдоподобие модели - благодаря нему компенсируются поколения сосудов, которыми мы пренебрегли в упрощенной схеме.

тип номер l, см d, см R, Па-с/мл с, см/с

Арт. 1 10 2 115 500

Арт. 2 20 1,6 600 700

Арт. 3 20 1,6 600 700

Арт. 4 5 2 10000 1000

Арт. 5 5 2 10000 1000

Вен. 1 10 3 25 350

Вен. 2 20 2,4 110 400

Вен. 3 20 2,4 110 400

Таблица. Параметры модели

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

Рассмотрим параметры гемодинамики в несколько моментов времени: в начале сердечного цикла, а также через 0,1; 0,2; 0.3; 0,4 и 0,5 с после. При установке кава-фильтра характер течения крови существенно меняется только в вене с имплантатом. В остальных сосудах возможен небольшой сдвиг показателей на константу, причем, чем больше сосуд удален от места установки кава-фильтра, тем меньше ощущается его влияние. Это объясняется тем, что организм стремится сгладить внешние воздействия. В вене 1, куда произведена имплантация, наблюдается замедление кровотока (рис. 4). Это происходит вследствие расширения сосуда: благодаря упругости кава-фильтра вена все время находится в распертом состоянии (рис. 5), а поток жидкости сохраняется прежним. Давление крови практически линейно по длине вены и существенных изменений не претерпевает (рис. 6).

и, см/с

Рис. 4. Величина скорости течения крови вдоль сосуда через 0 (линия 0); 0,1 (1); 0,2 (2); 0,3 (3); 0,4 (4); 0,5 с (5) после начала сердечного цикла

Рис. 6. Величина давления по длине сосуда через 0 (линия 0); 0,1 (1); 0,2 (2); 0,3 (3); 0,4 (4); 0,5 с (5)

после начала сердечного цикла

Таким образом, кава-фильтр типа «Зонтик» вызывает замедление кровотока в области его установки. Сильное уменьшение скорости течения становится причиной нового тромбообразования, поэтому необходимо задавать имплантату такие параметры, чтобы он был «прозрачен» для кровотока. С помощью объединения моделей глобальной циркуляции крови и эластичной стенки сосуда становится возможным наблюдение за влиянием кава-фильтров с различными параметрами на гемодинамику, а также определение их оптимальных размеров.

Работа поддержана РФФИ, грант 08-01-00159, 09-01-00115, а также программой Президиума РАН 21-П «Фундаментальные науки в медицине».

Литература

1. Есикова Н.Б., Мухин С.И., Соснин Н.В., Фаворский А.П., Хруленко А.Б. Математическое моделирование течения крови с кава-фильтрами. Препринт. - М.: МАКС Пресс, 2004.

2. Холодов А.С., Симаков С.С. Численное исследование содержания кислорода в крови человека при низкочастотных воздействиях // Математическое моделирование. - 2008. - Т. 20. - № 4. - С. 87-102.

3. Каро К., Педли Т., Шротер Р., Сид У. Механика кровообращения / Под ред. С.А. Регирер, В.М. Хаю-тина. - М.: Мир, 1981.

4. Rosar M.E., Peskin C.S. Fluid Flow in Collapsible Elastic Tubes: A Three-Dimensional Numerical Model // New York J. Math. - 2001. - № 7. - Р. 281-302.

5. Савельев В.С., Гологорский В.А., Кириенко А.И. и др. Флебология / Под ред. В.С. Савельева - М.: Медицина, 2001.

6. Cheng Tu and Charles Peskin. Stability and Instability in the Computation of Flows with Moving Immersed Boundaries: A Comparsion of Three Methods // SIAM J. SCI. STAT. COMPUT. - November 1992. - Vol. 13. - № 6. - Р. 1361-1376.

Иванов Юрий Александрович - Московский государственный университет им. М.В. Ломоносо-

ва, студент, yura-vtm@yandex.ru Добросердова Татьяна Константиновна - Московский государственный университет им. М.В. Ломоносова, студент, DobroserdovaTK@gmail.com

УДК 519.711

ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ РАСПРОСТРАНЕНИЯ ТУБЕРКУЛЕЗА В.Н. Леоненко, К.К. Логинов

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

Введение

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

Индивидуум-ориентированные модели позволяют детально провести формализацию принципов передачи и развития туберкулеза и исследовать динамику заболеваемости туберкулезом от различных факторов. Одними из первых в этом направлении являются модели, представленные в работах [3, 4]. В работе [3] предложена стохастическая модель в форме высокоразмерной системы рекуррентных уравнений с дискретным временем, учитывающая зависимость развития заболевания от возраста индивидуумов. Модель, построенная в работе [4], имеет следующие отличительные свойства: используется непрерывное время, продолжительность пребывания индивидуумов в различных стадиях болезни задается с помощью случайных величин с произвольными законами распределения, вероятность контактов индивидуумов считается пропорциональной произведению численностей когорт.

В данной работе рассматривается стохастическая модель [5], созданная на основе дифференциальной модели динамики туберкулеза [2]. В индивидуум-ориентированной модификации модели, в отличие