Математическая модель роста тромбоцитарного тромба и приложения к моделированию хронических болезней почек

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

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

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

4. Алгоритм, позволяющий автоматически изменять параметры модели, записанной в терминах обыкновенных дифференциальных уравнений, так чтобы решение системы ОДУ было максимально близко к результатам экспериментальных данных.

Теоретическая и практическая ценность

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

2. Исследована устойчивость некоторых автоволновых режимов (останавливающийся импульс, спиральная волна) в математической модели свертывания крови с учетом переключения активности тромбина.

3. Разработан алгоритм, позволяющий автоматически находить оптимальные коэффициенты системы ОДУ. Его можно использовать при расчетах различных задач медицины и биологии, где зачастую различные параметры известны с точностью до порядка. По результатам расчетов с использованием алгоритма была существенно изменена математическая модель адсорбции альбумин-билирубинового комплекса на поверхности УПП. Кроме того, алгоритм успешно использовался для нахождения оптимальных коэффициентов в задачах горения.

Методы исследований

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

Моделирование производилось в среде программирования Microsoft Visual Studio 2003, а также с использованием операционной системы Linux (расчеты на компьютерном кластере).

Апробация и публикации

Работа изложена на 108 страницах, состоит из введения, четырех глав, заключения и списка использованных источников, включающего 86 наименований.

2. СОДЕРЖАНИЕ РАБОТЫ

Во ВВЕДЕНИИ обосновывается актуальность темы диссертации, формулируется ее цель, представляются результаты, выносимые на защиту, а также определяется научная новизна, теоретическая и практическая значимость полученных результатов.

В ГЛАВЕ 1 представлен краткий обзор публикаций по математическому моделированию тромбоцитарного и плазменного путей тромбообразования.

В ГЛАВЕ 2 описана разработанная с участием автора математическая модель образования тромбоцитарного тромба.

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

Будем считать, что скорость течения достаточно велика, чтобы переносом тромбоцита вдоль оси сосуда можно было пренебречь по сравнению с переносом за время между столкновениями тромбоцитов.

Частота столкновений тромбоцитов в кровотоке н зависит от размера тромбоцитов, их концентрации в кровотоке, скорости течения и удаления от оси сосуда:

В этом интеграле относительная скорость Дv(r,x) может быть вычислена, если полагать, что скорость течения крови в сосуде описывается решением уравнений механики сплошной среды, в рамках данной модели -- уравнений Навье-Стокса.

В модели использованы упрощенные представления об эволюции состояния тромбоцитов. Предполагается, что следует различать тромбоциты по следующим признакам: пассивные-активные и полные-пустые. При этом для перехода из пассивного в активное состояние требуется (поглощается тромбоцитом) активатор (тромбин), а при переходе из полного в пустое состояние активатор выбрасывается тромбоцитом в кровоток.

Кинетика переходов в приведенной схеме описывается системой четырех уравнений, в которых приняты следующие обозначения:

w -- концентрация активатора;

cp -- концентрация пассивных тромбоцитов;

cf -- концентрация активных тромбоцитов;

c -- концентрация всех активных тромбоцитов;

k1, k2, kw -- константы скоростей мономолекулярных реакций;

k, w0, m -- константы функции, определяющей характер реакции перехода пассивных тромбоцитов в активные при взаимодействии с активатором (k -- коэффициент реакции, w0 -- константа Михаэлиса, m -- степень):

Предполагается, что изменения концентрации описываются уравнениями

где через обозначена проекция вектора градиента концентрации на направление вектора скорости течения жидкости, -- проекция вектора градиента на направление, нормальное к направлению скорости, мало, а , где н определяется выражением (1).

Значения параметров: kw = 2,0, k1 = 15,0, k2 = 0,5, k = 20,0, m = 2.0.

Для уравнений переноса тромбоцитов в потоке задаются граничные условия. На входе в сосуд заданы значения всех концентраций как функции времени, на выходе -- свободные (неотражающие) граничные условия. На оси сосуда ставятся условия симметрии, на стенке сосуда диффузионный поток активатора (тромбина) по нормали к стенке равен нулю , поток тромбоцитов на активированный участок стенки (см. ниже) определяется соотношением где -- значение концентрации тромбоцитов на расстоянии от стенки сосуда.

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

Рис. 1. Примерный вид сосуда и сечения расчетной области.

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

· течение вязкой жидкости в сосуде переменного сечения;

· моделирование тромбообразования в случае почечной патологии, когда все тромбоциты в потоке крови можно считать активированными;

· моделирование тромбообразования для полной системы уравнений с учетом распределения тромбоцитов по состояниям.

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

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

Описание методики расчета гемодинамического течения проводится с использованием системы уравнений Навье-Стокса, записанной в цилиндрической системе координат. Течение предполагается осесимметричным, поэтому задача рассматривается в r-z геометрии. В основу метода положена модель искусственной сжимаемости, в которой эллиптическое уравнение неразрывности заменяется гиперболическим уравнением.

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

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

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

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

Рис. 2. На рисунке изображено распределение (в изолиниях) радиальной компоненты скорости (Vr) для случая сформировавшегося возвратного течения за препятствием. Re~100. Стрелками показано векторное поле скоростей.

Рис. 3. На рисунке изображено распределение (в изолиниях) осевой компоненты скорости (Vz) для случая сформировавшегося возвратного течения за препятствием. Re~100. Стрелками показано векторное поле скоростей.

Кроме того, был выполнен расчет течения в стенозированном сосуде при больших числах Рейнольдса. За препятствием также образуется зона возвратного течения (рис. 2 и 3).

Для решения системы (2-5) применяется метод разделения по физическим процессам. На первом шаге считается реакционная часть методом Гира, на втором шаге считается диффузионная часть вариационно-разностным методом.

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

где с -- концентрация всех активных тромбоцитов.

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

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

Рис. 4. Распределение концентрации активных тромбоцитов (трехмерная поверхность), граница сосуда и граница активной зоны (жирные линии) и поле скоростей потока жидкости (стрелки). Число Рейнольдса Re~1. Рис. 5. Распределение концентрации активных тромбоцитов (трехмерная поверхность), граница сосуда и граница активной зоны (жирные линии) и поле скоростей потока жидкости (стрелки). Число Рейнольдса Re~100.

ГЛАВА 3 посвящена исследованию устойчивости некоторых режимов в математической модели свертывания крови с учетом переключения активности тромбина. Для описания модели используется система уравнений, записанных в безразмерном виде, приведенная ниже. Первое уравнение системы описывает изменения концентрации активатора системы свертывания -- тромбина (переменная u). Тромбин является автокаталитической переменной. Второе уравнение описывает изменение концентрации ингибитора свертывания -- протеина C (переменная v). Третья переменная модели (w) соответствует XI фактору свертывания.

,

,

Для исследования устойчивости используется аппарат сопряженных уравнений. Решение сопряженной системы позволяет выделить класс возмущений, к которым наиболее чувствительна основная задача. Поведение нормы решения сопряженной задачи позволяет судить об устойчивости решения основной задачи к данному классу возмущений (в норме L2). В работе исследована устойчивость останавливающихся импульсов тромбина, а также спиральных волн.

Для нелинейной задачи сопряженная система может быть построена неединственным образом. Для исследования устойчивости получающихся стационарных решений задачи (6-8) сопряженная система нелинейных уравнений строилась с помощью формулы Тейлора.

Сопряженная система имеет следующий вид:

где функции , и могут быть представлены в виде интегралов

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

На рис. 6 показан класс возмущений v1, к которым наиболее чувствительны решения типа спиральной волны. Оно имеет максимум в окрестности центра спиральной волны. Норма решения сопряженной задачи убывает со временем, что вполне согласуется с устойчивостью автоволн при данных значениях параметров. Отметим, что такие возмущения не разрушают спиральные волны в системе, а переводят одну спиральную волну в другую, с близким расположением центра волны.

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

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

Рис. 6. Характерное распределение первой переменной сопряженной задачи для решения типа спиральной волны.

В ГЛАВЕ 4 рассматриваются некоторые проблемы моделирования современных установок диализа. В частности, рассматривается математическая модель взаимодействия угольного пирополимера (УПП) с искусственной системой, включающей альбумин и билирубин. Заболевания печени стоят на 8-9-м местах среди причин смертности населения репродуктивного возраста в развитых странах, что делает разработку установок диализа актуальной проблемой. Создано несколько методов лечения, базирующихся на удалении белоксвязанных веществ из асцитической жидкости, плазмы крови и цельной крови с помощью этих адсорбционных материалов. Тем не менее, многие детали взаимодействия в системе альбумин-НБ-УПП до сих пор не ясны.

Для описания химических реакций процесса используется следующая система обыкновенных дифференциальных уравнений

здесь x1 - комплекс альбумин-билирубин (связь через первичный сайт), x2 - комплекс билирубин-УПП, x3 - комплекс альбумин-УПП, x4 - комплекс альбумин-билирубин (связь через первичный и вторичный сайты). Константы реакций и C0 ранее определялись с использованием программы dBsolve 5 по экспериментальным данным.

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

Таблица 1а. Экспериментальные данные №1 (n = 1).

В ходе расчетов производилась оптимизация до трех параметров задачи. Результаты расчетов приведены в таблицах 1а, б и на рис. 7а, б.

Таблица 1б. Экспериментальные данные №2 (n = 5).

Рис 7а. Экспериментальные данные №1. Штрихованная линия соответствует решению системы (15-18) с новыми значениями параметров k5, k7 и k9, сплошная линия -- k5, k8 и k12. Рис 7б. Экспериментальные данные №2. Штрихованная линия соответствует решению системы (15-18) с новыми значениями параметров k5, k8 и k12, сплошная линия -- k5, k8 и k10.

Как видно из рис. 8, значение функционала среднеквадратичного отклонения уменьшается с увеличением количества варьируемых параметров, а решение системы ОДУ (9-12) становится все ближе к экспериментальным данным. В ходе расчетов по оптимизации параметров было обнаружено, что один из параметров, а именно k8 становится меньше нуля. Это означает, что вторая реакция в системе (9-12) является автокаталитической. Данный факт привел к необходимости существенного изменения модели и ее уравнений. В настоящее время проводится работа по созданию улучшенной модели, которая бы учитывала автокаталитическое поведение билирубина.

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

Рис. 8. Оптимизация функционала среднеквадратичного отклонения.

В ЗАКЛЮЧЕНИИ формулируются основные результаты работы.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ

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

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

3. На основе оптимизации параметров в системе, описывающих процессы в гемодиализных колонках, сделан вывод о необходимости модернизации описания механизма химических реакций. Предположительно, при адсорбции альбумин-билирубинового комплекса на поверхности угольного пирополимера (УПП) необходимо учитывать активизацию соседних сайтов связывания на поверхности УПП и вводить в рассмотрение автокаталитичность реакций связывания.

ПУБЛИКАЦИИ ПО ТЕМЕ ДИССЕРТАЦИИ

1. Лобанов А. И., Украинец А. В. Исследование устойчивости стационарных структур тромбина в математической модели свертывания крови. // МАТЕМАТИКА, КОМПЬЮТЕР, ОБРАЗОВАНИЕ. Сборник научных трудов. -Москва-Ижевск.-2004. №11, ч.2, -C.829-836.

2. Украинец А. В. Хаотическая динамика в модели свертывания крови. // Современные проблемы фундаментальных и прикладных наук: Сборник трудов 47-й научной конференции МФТИ, Т. III / Моск. физ.-техн. ин-т.- М., 2004. -С.48-49.

3. Украинец А. В. Учет клеточной структуры в модели свертывания крови. // Современные проблемы фундаментальных и прикладных наук: Сборник трудов 48-й научной конференции МФТИ, Т.III / Моск. физ.-техн. ин-т.- М., 2005. -С.24-25.

4. Буравцев В. Н., Украинец А. В. Математическая модель роста тромбоцитарного тромба. // Аннотации докладов IX Всероссийского съезда по теоретической и прикладной механике. -Т. I.-Н. Новгород, 2006. -С.134.

5. Украинец А. В., Николаев А. В., Лобанов А. И. Об уточнении механизма реакции адсорбции на основе процедуры автоматической оценки параметров системы ОДУ. // Современные проблемы фундаментальных и прикладных наук: Сборник трудов 49-й научной конференции МФТИ, Т. II / Моск. физ.-техн. ин-т.- М., 2006. -С.66-67.

6. Buravtsev V. N., Lobanov A. .I., Ukrainets A. V. Mathematical model of platelet thrombus formation. // Proceedings of the International Conference «Numerical geometry, grid generation and high performance computing» / A. A. Dorodnicyn Computing Center Russian Academy of Sciences, Moscow, 2008. -P.33-36.

7. Куриленко И. А., Лобанов А. И., Украинец А. В. Численное исследование устойчивости некоторых автоволновых решений математической модели свертывания крови. // Биофизика -2009. -т.54, вып.1 -C.68-76.

8. Буравцев В. Н., Лобанов А. И., Украинец А. В. Математическая модель роста тромбоцитарного тромба. // Математическое моделирование -2009. -т.21, №3 -C.109-119.

В работах с соавторами лично соискателем выполнены: [1, 7] -- исследование устойчивости автоволновых режимов в модели ССК с использованием аппарата сопряженных уравнений; [4, 6, 8] -- исследование основных закономерностей формирования тромбоцитарного тромба в потоке вязкой жидкости; [5] -- построен алгоритм для автоматической оптимизации параметров системы ОДУ, сделан вывод о необходимости модернизации описания механизма химических реакций.

Размещено на Allbest.ru