МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЛИЯНИЯ УСТАНОВКИ КАВА-ФИЛЬТРА НА ГЕМОДИНАМИКУ КРОВЕНОСНОЙ СИСТЕМЫ
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЛИЯНИЯ УСТАНОВКИ КАВА-ФИЛЬТРА...
9 БИОМЕДИЦИНСКИЕ ТЕХНОЛОГИИ
УДК 51.76
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЛИЯНИЯ УСТАНОВКИ КАВА-ФИЛЬТРА НА ГЕМОДИНАМИКУ КРОВЕНОСНОЙ СИСТЕМЫ
Ю.А. Иванов, Т.К. Добросердова
Моделирование кровеносной системы является обширной областью для научной деятельности и включает в себя множество аспектов. В этом смысле вопрос о строении сосуда и его реакции на деформацию логично рассматривать отдельно. На основе биологических особенностях строения ткани реализуется волоконная модель эластичной стенки сосуда, идея которой была сформулирована в работах Чарльза Пескина [6]. Подробно рассматривается калибровка модели с помощью специальных входных данных – уравнений состояния сосуда. Отдельная часть работы посвящена интеграции волоконной модели эластичности и глобальной модели кровообращения. Взаимодействие моделей дает возможность исследования влияния установки модельного кава-фильтра на гемодинамику. Ключевые слова: тромбоэмболия, трансмуральное давление, уравнение состояния сосуда, кава-фильтр.
Введение
В настоящее время для описания динамики кровообращения создан ряд численных гидродинамических моделей, в основе которых лежит система уравнений в частных производных. Огромный интерес для науки представляет исследование влияния на гемодинамику различных факторов, в частности, некоторых болезней. Тромбоэмболия – смертельно опасное осложнение, при котором происходит закупорка сосудов тромбами. Борьба с ней входит в число главных проблем медицины. Одним из методов профилактики эмболии является установка в вены специальных имплантатов – кава-фильтров. Однако в этом случае возможно не только формирование областей замедленного кровотока, но и повреждение стенки сосуда или ее разрыв в случае неправильного подбора размеров конструкции или ее установки. Все это приводит к повторному тромбообразованию и эмболии легочной артерии не только тромбами, но и обломками кава-фильтра (в случае его разлома). Таким образом, крайне важно отследить влияние имплантанта на гемодинамику. В работе [1] с помощью гидродинамической модели глобальной циркуляции крови исследовано возмущение кровотока в случае попадания в кава-фильтр тромба посредством уменьшения эффективного сечения вены. Мы же исследуем влияние установленной конструкции, модифицировав одно из уравнений основной системы с помощью волоконной модели эластичной стенки сосуда. Данная технология позволяет рассчитывать силы, возникающие в стенке деформированного сосуда (в том числе при наличии кава-фильтра), фиксировать изменения его формы при повышении или понижении кровяного давления.
Модель глобального кровообращения
Предполагается, что кровь – вязкая несжимаемая жидкость с постоянной температурой, протекающая по сети эластичных трубок. Трубки связаны в четыре графа, задающих артериальную и венозную часть большого и малого кругов кровообращения. Введем следующие обозначения: t – время; x – координата по длине сосуда; S – площадь поперечного сечения сосуда; u – скорость кровотока, осредненная по поперечному сечению; ρ – плотность крови; p – трансмуральное давление, равное разности давлений в сосуде и того, которое оказывают на него окружающие ткани; c – скорость распространения малых возмущений; Sˆ – площадь поперечного сечения при p = 0; d – диаметр трубки; µ – коэффициент вязкости
крови. Для каждого сосуда будут выполняться следующие уравнения:
S t
S (u) x
0;
(1)
u t
u2 2
x
p
t,
x,
S,u;
(2)
p c2 f (S),
(3)
где
f
(S)
exp(SSˆ 1 1)
ln(SSˆ 1 ),
1, S
S Sˆ
Sˆ
.
Выражения
(1),
(2)
отражают
законы
сохранения
массы
и
количества
движения соответственно. Функция ψ задает вязкое трение. Уравнение (3) учитывает эластичные свойст-
ва стенок сосудов. В настоящей работе будет использоваться гидродинамическая модель, основанная на
системе (1)–(3) [2].
94 Научно-технический вестник Санкт-Петербургского государственного университета
информационных технологий, механики и оптики, 2010, № 4(68)
Ю.А. Иванов, Т.К. Добросердова
Решения задачи, полученные на каждом ребре каждого из четырех графов отдельно, необходимо сшить, задав граничные условия на концах сосудов. В случае, когда концевая вершина принадлежит точке стыковки вен с венами или артерий с артериями, используются закон Пуазейля и закон сохранения массы:
pk (S(t, x)) pi (t) k Rk Sk (t, x)uk (t, x),
m
k Sk (t, x)uk (t, x) 0,
k 1
где k – номер сосуда, pi – давление в точке стыковки; i 1 и x 0 , если сосуд выходит из этой точки, = 1 и x l в противном случае (l – длина сосуда); R – сопротивление, которое преодолевает жидкость, протекающая через эту область; m – количество сосудов, имеющих данный узел своим концом, суммирование проводится по всему множеству. Замыкается система уравнением совместности – условием, накладываемым характеристическими кривыми, идущими внутрь области интегрирования. Поскольку размеры артериол и венул очень малы, их моделирование невозможно и не необходимо, поэтому в узлах стыковки артерий с венами используется такая же система граничных условий, как и в предыдущем случае, с той лишь разницей, что сопротивление R берется значительно больше, заменяя поколения сосудов, которыми приходится пренебречь. В настоящей работе используется динамическая модель четырехкамерного сердца [2]. Система уравнений, лежащая в ее основе, дополненная условиями равенства давлений и потоков в соответствующих предсердиях (желудочках) и во входящих (выходящих) в сердце венах (артериях), а также уравнениями совместности в концевых точках этих же сосудов, определяет систему граничных условий в узлах сопряжения сосудистых деревьев с сердцем.
Рассмотрим более подробно уравнение состояния (3). Поскольку при постановке задачи о глобальном кровотоке предполагается, что сосуды – эластичные трубки, то каждому из них необходимо сопоставить особый параметр, отражающий эти свойства. Таким параметром является скорость распространения волн давления малой амплитуды по упругой трубке (в наших обозначениях – c). В зависимости от нее график уравнения состояния, представляющий собой S-образную кривую, имеет разный угол наклона. Точка перегиба в нем соответствует переходу от круглой формы сечения сосуда к эллипсообразной.
Функция f(S) выведена эмпирически [3], поэтому никак не отражает воздействий на сосуд. Попытаемся получить новое уравнение состояния, используя волоконную модель эластичной стенки.
Модель эластичной стенки сосуда
Волоконная модель эластичной стенки сосуда основывается на биологических особенностях ее строения [6]. Она представляет собой систему переплетенных волокон, обладающих упругими свойствами.
Рис. 1. Волоконная модель эластичной стенки
Существует три типа волокон: эластические, гладкомышечные и коллагеновые. Модули Юнга первых двух примерно равны 3105 Па и в данной модели будут объединены в один класс. Эластические волокна образуют достаточно густую сеть, особенно на внутренней оболочке стенки. Они создают эластическое напряжение, которое противодействует кровяному давлению, растягивающему сосуд. Коллагеновые волокна имеют модуль Юнга 1108 Па. Сеть, образованная ими в средней и наружной оболочках, оказывает сопротивление растяжению большее, чем сеть эластических волокон [4]. Коллагеновые волокна располагаются относительно свободно в стенке сосуда, иногда даже образуя складки, так что противодействие внутреннему давлению оказывают лишь тогда, когда сосуд уже растянут до определенного состояния. Модельная реализация описанной выше структуры представляет собой набор нескольких пар спиралевидных (геликовидных) волокон, закрученных вдоль стенки в противоположных направлениях (рис. 1). Каждая группа определяется своим набором параметров – модулем упругости, скоростью реакции на деформацию, степенью закрученности и т.д. Также в модель могут быть введены и другие типы волокон, например, кольцевидные, опоясывающие сосуд по всей длине, или прямолинейные (лучевидные).
Определим силы упругости, возникающие вследствие деформации волокна, заданного набором точек. Пусть s – переменная Лагранжа, определяющая расстояние вдоль волокна от данной точки до не-
которой начальной; X (s, t) – вектор-функция, описывающая положение физической точки с лагранжевой
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 4(68)
95
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЛИЯНИЯ УСТАНОВКИ КАВА-ФИЛЬТРА...
координатой s; t – время; Т – напряжение деформированного волокна. Предположим, что напряжение подчиняется нелинейному обобщению закона Гука:
T
X s
T*
X s
1
,
(4)
где T* – коэффициент упругости волокна. Тогда локальная сила упругости f , действующая со стороны
стенки, выражается следующим образом:
f
s
T
;
(5)
X s
X s
,
(6)
где – направление касательной в точке волокна.
Осуществим дискретизацию уравнений (4)–(6). Пусть Nb – количество расчетных точек на волокне; Xk – координата k-ой точки в трехмерном пространстве, где k = 1,...,Nb; S – расстояние между соседними точками волокна в недеформированном состоянии. Предполагаем, что S будет одинаково для точек,
принадлежащих однотипным волокнам. Далее в соответствии с формулами (4), (6) дискретизируем T и следующим образом:
Tk 12
X k 1 X k / S
X k 1 X k S
1 ;
(7)
k 12
X k 1 X k X k 1 X k
.
(8)
Дискретная сила упругости в k-ой точке определяется по формуле
fk
T k 12 k 12
T k
12
k
12
S
.
(9)
Реализуем модель эластичной стенки в виде двух наборов спиралевидных волокон с двумя различ-
ными коэффициентами упругости и параметрами, характеризующими степень закрученности волокон.
Первый набор, моделирующий связки эластина, реагирует на любую деформацию, второй, моделирую-
щий коллагеновые волокна – только если деформация превышает определенный порог. Тогда калибровка
модели сводится к отысканию пяти параметров: Te, Tc – коэффициенты упругости двух наборов волокон; Le, Lс – «шаги» спирали; Rc – параметр, определяющий степень деформации, до которой второй набор волокон (коллагена) не вносит вклада в реакцию упругой стенки сосуда.
Рассмотрим следующую задачу. Пусть в модель на вход подаются данные – уравнение состояния
некоторого сосуда. Оно может быть задано как дискретный набор данных (полученных в результате экс-
перимента), так и в виде некоторой функции, описывающей уравнение состояния реального сосуда. Зада-
ча – подобрать коэффициенты так, чтобы уравнение состояние нашей модели было наиболее близко в
смысле L2-нормы к целевой функции на заданном отрезке, т.е. на определенном диапазоне изменения сечения сосуда.Под уравнением состояния модели и под целевой функцией в случае дискретных входных
данных будем понимать кусочно-линейную интерполяцию. Предположим, что у нас есть только один
набор волокон, описываемый двумя параметрами T и L. Требуется отыскать их значения, чтобы наилучшим образом удовлетворить целевой функции P(S) на отрезке [S0, S1], и где P(S0)=0. Функция P(S) для реальных сосудов имеет экспоненциальный характер роста на данном промежутке. Заметим, что уравне-
ние состояния модели имеет линейный характер при малых значениях параметра l и усиливает выпуклость при его увеличении. Параметр T будет играть роль нормирующего коэффициента. Основываясь на
этих замечаниях и применяя метод золотого сечения, строим итерационный алгоритм, минимизирующий
L2-норму ошибки на заданном отрезке. Отыскание совокупности описанных ранее пяти параметров модели – Te, Tc, Le, Lс и Rc – основывается на тех же принципах. Если рассматривать растяжение сосуда при S
9 БИОМЕДИЦИНСКИЕ ТЕХНОЛОГИИ
УДК 51.76
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЛИЯНИЯ УСТАНОВКИ КАВА-ФИЛЬТРА НА ГЕМОДИНАМИКУ КРОВЕНОСНОЙ СИСТЕМЫ
Ю.А. Иванов, Т.К. Добросердова
Моделирование кровеносной системы является обширной областью для научной деятельности и включает в себя множество аспектов. В этом смысле вопрос о строении сосуда и его реакции на деформацию логично рассматривать отдельно. На основе биологических особенностях строения ткани реализуется волоконная модель эластичной стенки сосуда, идея которой была сформулирована в работах Чарльза Пескина [6]. Подробно рассматривается калибровка модели с помощью специальных входных данных – уравнений состояния сосуда. Отдельная часть работы посвящена интеграции волоконной модели эластичности и глобальной модели кровообращения. Взаимодействие моделей дает возможность исследования влияния установки модельного кава-фильтра на гемодинамику. Ключевые слова: тромбоэмболия, трансмуральное давление, уравнение состояния сосуда, кава-фильтр.
Введение
В настоящее время для описания динамики кровообращения создан ряд численных гидродинамических моделей, в основе которых лежит система уравнений в частных производных. Огромный интерес для науки представляет исследование влияния на гемодинамику различных факторов, в частности, некоторых болезней. Тромбоэмболия – смертельно опасное осложнение, при котором происходит закупорка сосудов тромбами. Борьба с ней входит в число главных проблем медицины. Одним из методов профилактики эмболии является установка в вены специальных имплантатов – кава-фильтров. Однако в этом случае возможно не только формирование областей замедленного кровотока, но и повреждение стенки сосуда или ее разрыв в случае неправильного подбора размеров конструкции или ее установки. Все это приводит к повторному тромбообразованию и эмболии легочной артерии не только тромбами, но и обломками кава-фильтра (в случае его разлома). Таким образом, крайне важно отследить влияние имплантанта на гемодинамику. В работе [1] с помощью гидродинамической модели глобальной циркуляции крови исследовано возмущение кровотока в случае попадания в кава-фильтр тромба посредством уменьшения эффективного сечения вены. Мы же исследуем влияние установленной конструкции, модифицировав одно из уравнений основной системы с помощью волоконной модели эластичной стенки сосуда. Данная технология позволяет рассчитывать силы, возникающие в стенке деформированного сосуда (в том числе при наличии кава-фильтра), фиксировать изменения его формы при повышении или понижении кровяного давления.
Модель глобального кровообращения
Предполагается, что кровь – вязкая несжимаемая жидкость с постоянной температурой, протекающая по сети эластичных трубок. Трубки связаны в четыре графа, задающих артериальную и венозную часть большого и малого кругов кровообращения. Введем следующие обозначения: t – время; x – координата по длине сосуда; S – площадь поперечного сечения сосуда; u – скорость кровотока, осредненная по поперечному сечению; ρ – плотность крови; p – трансмуральное давление, равное разности давлений в сосуде и того, которое оказывают на него окружающие ткани; c – скорость распространения малых возмущений; Sˆ – площадь поперечного сечения при p = 0; d – диаметр трубки; µ – коэффициент вязкости
крови. Для каждого сосуда будут выполняться следующие уравнения:
S t
S (u) x
0;
(1)
u t
u2 2
x
p
t,
x,
S,u;
(2)
p c2 f (S),
(3)
где
f
(S)
exp(SSˆ 1 1)
ln(SSˆ 1 ),
1, S
S Sˆ
Sˆ
.
Выражения
(1),
(2)
отражают
законы
сохранения
массы
и
количества
движения соответственно. Функция ψ задает вязкое трение. Уравнение (3) учитывает эластичные свойст-
ва стенок сосудов. В настоящей работе будет использоваться гидродинамическая модель, основанная на
системе (1)–(3) [2].
94 Научно-технический вестник Санкт-Петербургского государственного университета
информационных технологий, механики и оптики, 2010, № 4(68)
Ю.А. Иванов, Т.К. Добросердова
Решения задачи, полученные на каждом ребре каждого из четырех графов отдельно, необходимо сшить, задав граничные условия на концах сосудов. В случае, когда концевая вершина принадлежит точке стыковки вен с венами или артерий с артериями, используются закон Пуазейля и закон сохранения массы:
pk (S(t, x)) pi (t) k Rk Sk (t, x)uk (t, x),
m
k Sk (t, x)uk (t, x) 0,
k 1
где k – номер сосуда, pi – давление в точке стыковки; i 1 и x 0 , если сосуд выходит из этой точки, = 1 и x l в противном случае (l – длина сосуда); R – сопротивление, которое преодолевает жидкость, протекающая через эту область; m – количество сосудов, имеющих данный узел своим концом, суммирование проводится по всему множеству. Замыкается система уравнением совместности – условием, накладываемым характеристическими кривыми, идущими внутрь области интегрирования. Поскольку размеры артериол и венул очень малы, их моделирование невозможно и не необходимо, поэтому в узлах стыковки артерий с венами используется такая же система граничных условий, как и в предыдущем случае, с той лишь разницей, что сопротивление R берется значительно больше, заменяя поколения сосудов, которыми приходится пренебречь. В настоящей работе используется динамическая модель четырехкамерного сердца [2]. Система уравнений, лежащая в ее основе, дополненная условиями равенства давлений и потоков в соответствующих предсердиях (желудочках) и во входящих (выходящих) в сердце венах (артериях), а также уравнениями совместности в концевых точках этих же сосудов, определяет систему граничных условий в узлах сопряжения сосудистых деревьев с сердцем.
Рассмотрим более подробно уравнение состояния (3). Поскольку при постановке задачи о глобальном кровотоке предполагается, что сосуды – эластичные трубки, то каждому из них необходимо сопоставить особый параметр, отражающий эти свойства. Таким параметром является скорость распространения волн давления малой амплитуды по упругой трубке (в наших обозначениях – c). В зависимости от нее график уравнения состояния, представляющий собой S-образную кривую, имеет разный угол наклона. Точка перегиба в нем соответствует переходу от круглой формы сечения сосуда к эллипсообразной.
Функция f(S) выведена эмпирически [3], поэтому никак не отражает воздействий на сосуд. Попытаемся получить новое уравнение состояния, используя волоконную модель эластичной стенки.
Модель эластичной стенки сосуда
Волоконная модель эластичной стенки сосуда основывается на биологических особенностях ее строения [6]. Она представляет собой систему переплетенных волокон, обладающих упругими свойствами.
Рис. 1. Волоконная модель эластичной стенки
Существует три типа волокон: эластические, гладкомышечные и коллагеновые. Модули Юнга первых двух примерно равны 3105 Па и в данной модели будут объединены в один класс. Эластические волокна образуют достаточно густую сеть, особенно на внутренней оболочке стенки. Они создают эластическое напряжение, которое противодействует кровяному давлению, растягивающему сосуд. Коллагеновые волокна имеют модуль Юнга 1108 Па. Сеть, образованная ими в средней и наружной оболочках, оказывает сопротивление растяжению большее, чем сеть эластических волокон [4]. Коллагеновые волокна располагаются относительно свободно в стенке сосуда, иногда даже образуя складки, так что противодействие внутреннему давлению оказывают лишь тогда, когда сосуд уже растянут до определенного состояния. Модельная реализация описанной выше структуры представляет собой набор нескольких пар спиралевидных (геликовидных) волокон, закрученных вдоль стенки в противоположных направлениях (рис. 1). Каждая группа определяется своим набором параметров – модулем упругости, скоростью реакции на деформацию, степенью закрученности и т.д. Также в модель могут быть введены и другие типы волокон, например, кольцевидные, опоясывающие сосуд по всей длине, или прямолинейные (лучевидные).
Определим силы упругости, возникающие вследствие деформации волокна, заданного набором точек. Пусть s – переменная Лагранжа, определяющая расстояние вдоль волокна от данной точки до не-
которой начальной; X (s, t) – вектор-функция, описывающая положение физической точки с лагранжевой
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 4(68)
95
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЛИЯНИЯ УСТАНОВКИ КАВА-ФИЛЬТРА...
координатой s; t – время; Т – напряжение деформированного волокна. Предположим, что напряжение подчиняется нелинейному обобщению закона Гука:
T
X s
T*
X s
1
,
(4)
где T* – коэффициент упругости волокна. Тогда локальная сила упругости f , действующая со стороны
стенки, выражается следующим образом:
f
s
T
;
(5)
X s
X s
,
(6)
где – направление касательной в точке волокна.
Осуществим дискретизацию уравнений (4)–(6). Пусть Nb – количество расчетных точек на волокне; Xk – координата k-ой точки в трехмерном пространстве, где k = 1,...,Nb; S – расстояние между соседними точками волокна в недеформированном состоянии. Предполагаем, что S будет одинаково для точек,
принадлежащих однотипным волокнам. Далее в соответствии с формулами (4), (6) дискретизируем T и следующим образом:
Tk 12
X k 1 X k / S
X k 1 X k S
1 ;
(7)
k 12
X k 1 X k X k 1 X k
.
(8)
Дискретная сила упругости в k-ой точке определяется по формуле
fk
T k 12 k 12
T k
12
k
12
S
.
(9)
Реализуем модель эластичной стенки в виде двух наборов спиралевидных волокон с двумя различ-
ными коэффициентами упругости и параметрами, характеризующими степень закрученности волокон.
Первый набор, моделирующий связки эластина, реагирует на любую деформацию, второй, моделирую-
щий коллагеновые волокна – только если деформация превышает определенный порог. Тогда калибровка
модели сводится к отысканию пяти параметров: Te, Tc – коэффициенты упругости двух наборов волокон; Le, Lс – «шаги» спирали; Rc – параметр, определяющий степень деформации, до которой второй набор волокон (коллагена) не вносит вклада в реакцию упругой стенки сосуда.
Рассмотрим следующую задачу. Пусть в модель на вход подаются данные – уравнение состояния
некоторого сосуда. Оно может быть задано как дискретный набор данных (полученных в результате экс-
перимента), так и в виде некоторой функции, описывающей уравнение состояния реального сосуда. Зада-
ча – подобрать коэффициенты так, чтобы уравнение состояние нашей модели было наиболее близко в
смысле L2-нормы к целевой функции на заданном отрезке, т.е. на определенном диапазоне изменения сечения сосуда.Под уравнением состояния модели и под целевой функцией в случае дискретных входных
данных будем понимать кусочно-линейную интерполяцию. Предположим, что у нас есть только один
набор волокон, описываемый двумя параметрами T и L. Требуется отыскать их значения, чтобы наилучшим образом удовлетворить целевой функции P(S) на отрезке [S0, S1], и где P(S0)=0. Функция P(S) для реальных сосудов имеет экспоненциальный характер роста на данном промежутке. Заметим, что уравне-
ние состояния модели имеет линейный характер при малых значениях параметра l и усиливает выпуклость при его увеличении. Параметр T будет играть роль нормирующего коэффициента. Основываясь на
этих замечаниях и применяя метод золотого сечения, строим итерационный алгоритм, минимизирующий
L2-норму ошибки на заданном отрезке. Отыскание совокупности описанных ранее пяти параметров модели – Te, Tc, Le, Lс и Rc – основывается на тех же принципах. Если рассматривать растяжение сосуда при S