CC BY
22
УДК 519.6
А.И. Лобанов, И.А. Куриленко, А.В. Украинец
Московский физико-технический институт (государственный университет)
Автоволновые решения и диссипативные структуры в двух математических моделях динамики свертывания крови
В статье приводится обзор и сравнение результатов численного моделирования динамики свертывания крови с использованием двух простых математических моделей.
Основу статьи составляют результаты, полученные на кафедре вычислительной математики МФТИ совместно с Гематологическим научным центром (ГНЦ РАМН). Приведены новые результаты по исследованию устойчивости автоволновых решений в математической модели с переключением активности тромбина. В качестве основного математического аппарата при проведении такого исследования использовалась сопряженная система, построенная с использованием формулы Тейлора.
Ключевые слова: автоволны, математическая модель, динамика свертывания крови, переключение активности тромбина, сопряженные уравнения.
I. Введение
Начиная с середины 90-х годов прошлого века большое количество теоретических и экспериментальных работ посвящено изучению одной из важнейших защитных систем организма — системы свертывания крови. Математическое моделирование процесса свертывания помогает осмыслить накопленный экспериментальный материал и выбрать правдоподобные гипотезы устройства этой сложной системы. Первая математическая модель свертывания [1] была предложена сразу же после высказывания гипотезы о каскадном устройстве свертывания [2]. Регулярное изучение кинетики системы свертывания крови методами математического моделирования началось с [3]. Среди многочисленных работ на эту тему отметим [4-6], где подробно проанализирована схема реакций каскада свертывания крови и приведена система ОДУ, описывающая кинетику свертывания. В настоящее время издаются специальные выпуски научных журналов, посвященные математическим моделям в динамике свертывания крови.
Каскад биохимических реакций свертывания охвачен петлями положительных и отрицательных обратных связей. Особого внимания заслуживает одна из реакций, в которой продукт нижней ступени
каскада может активировать фактор, расположенный на верхней ступени каскада. Ранее [7] было высказано предположение, что подобные реакции в системе свертывания могут приводить к автономному са-моподдерживающемуся образованию факторов свертывания при рассмотрении пространственной динамики системы. Эта гипотеза была проверена с помощью вычислительных экспериментов [8, 9]. Было показано, что введение в уравнения реакции активации XI фактора свертывания тромбином, не влияя на гомогенную кинетику свертывания из-за малой скорости, тем не менее драматически меняет характер пространственного роста тромба, делая его безостановочным. В пространственно-распределенной системе при этом наблюдается автоволновое распространение активных факторов.
Для получения в модели фибринового сгустка конечного размера разные авторы привлекают дополнительные гипотезы. Так, авторы [10] предполагают, что основным механизмом, позволяющим останавливать процесс свертывания, является взаимодействие с гидродинамическим потоком. Очевидным недостатком такого подхода является его противоречие экспериментальным данным. Показано [11-13], что в экспериментах in vitro в отсутствие конвективных потоков рост тромбов прекращается. При этом образующиеся тром-
бы могут иметь сложную структуру. Могут образоваться слоистые тромбы, в том числе, с конечным числом слоев. Таким образом, даже в однородной среде без гидродинамических потоков автоволна свертывания может останавливаться ранее, чем свернется вся имеющаяся в экспериментальной установке плазма крови. Включение в модель только известных реакций свертывания не приводит к описанию наблюдаемых явлений. Для объяснения результатов экспериментов in vitro необходимо вводить дополнительные гипотетические механизмы относительно устройства системы свертывания. Детальное обсуждение этой проблемы проведено в обзоре [14].
Первой гипотезой было допущение о наличии в системе реакций дополнительного белка, обеспечивающего самоускоря-ющийся характер производства ингибитора системы. Отметим, что в настоящее время в системе свертывания кандидат на роль такого белка неизвестен. Тем не менее гипотеза [7] позволила сформулировать феноменологическую модель, описывающую пространственные аспекты динамики свертывания крови [12, 15]. В настоящее время свойства этой математической модели хорошо изучены [16-20], она вошла составной частью в более сложные математические модели, описывающие динамику тромбообразования в потоке жидкости [21, 22]. Введение автокаталитиче-ского характера производства ингибитора приводит к тому, что автоволна свертывания из-за взаимодействия с ингибитором может останавливаться на конечном расстоянии от места активации [16-19].
Для практических целей и экспериментальных исследований одной феноменологической модели явно недостаточно. Необходимо иметь последовательность моделей, с разной степенью детализации описывающих динамику свертывания крови.
Модель [8, 9], базирующаяся на известных биохимических реакциях свертывания, показала, что вопрос об остановке пространственного роста фибринового сгустка не описывается данной системой уравнений. Позднее были получены экспериментальные данные [23, 24], свидетельствующие о том, что тромбин может существовать в нескольких формах: проко-агулянтной, ускоряя свертывание, и анти-коагулянтной, замедляя процесс собствен-
ного образования. Показано [25], что, вводя в модель [8, 9] механизм переключения между этими двумя состояниями тромбина, возможно описать экспериментально наблюдаемый характер пространственного роста тромба, включая его остановку. Подробное описание модели с последующей её редукцией для качественного анализа наблюдаемых динамических режимов тромбообразования изложены в [25]. Сравнение ряда свойств феноменологической двухволновой модели и математической модели с переключением активности тромбина проведено в [26], многие одномерные результаты, приведённые ниже, также процитированы из [26].
В [25, 26] исследовалась в основном одномерная постановка задачи. Для двухволновой модели многие двумерные результаты получены ранее численно. Для математической модели с переключением активности тромбина оставались невыяснённы-ми вопросы устойчивости некоторых типов решений в двумерном случае, например, бегущих импульсов тромбина, останавливающихся в однородной среде на конечном расстоянии от места активации. Тем более, что в пространственно-однородной (точечной) постановке эти режимы соответствовали хаотическим колебаниям (странному аттрактору, близкому по своим свойствам аттрактору Реслера). Поэтому представляет существенный интерес рассмотрение двумерной постановки задачи и исследование свойств устойчивости к малым возмущениям начальных данных решений уравнений математической модели с переключением активности тромбина.
Математические модели. Оговоримся сразу, что в отличие от ряда современных работ, таких, как [27], в рассматриваемые модели мы не включаем активацию тромбоцитов и ограничиваемся процессами, происходящими в плазме крови. Некоторые предварительные результаты расчётов по математической модели формирования тромбоцитарного тромба приведены в [28]. Оговоримся также, что результаты, полученные для систем типа «реакция-диффузия-конвекция», останутся за пределами данной статьи.
Уравнения первой из рассматриваемых моделей — феноменологической модели свертывания [9] — описывают изменение концентрации двух метаболитов — актива-
тора процесса свертывания и (тромбина) и ингибитора V (предположительно протеина С):
,2
ди ^. аи2
— = Б А и + —— дъ и + и0
— К1и —
(1)
Коэффициенты диффузии метаболитов предполагаются равными. Это связано с
тем, что и активатор, и ингибитор в системе — белки с примерно одинаковыми массами молекул. В пределах точности экспериментов измерения этих коэффициентов дают для них одинаковые значения. Характерные значения констант, входящих в систему (1-2), приведены в табл. 1. Отметим, что в более поздних работах значения коэффициентов системы подверглись уточнению [29].
Таблица 1
-1 а,мин /3,мин 1 _1 7,мин Ио,нМ г’о ,нМ С, нМ _1 ^1,мин _1 К‘2 ,МИН
2,0 0,0015 5,0 2,95 0,0525 5,0 0,05 0,35
Параметр а (точнее, а — К\) имеет смысл обратного АЧТВ (активированного частичного тромбопластинового времени) — одного из характерных времен, измеряемых во всех анализах крови. Остальные константы, в том числе пороговое значение концентрации ингибитора Vo и характерное значение концентрации активатора и0 оценены по экспериментальным данным и использовались для математического моделирования процессов свертывания в [16-22].
При полном перемешивании (Д = 0) модель переходит в систему обыкновенных дифференциальных уравнений (ОДУ). Такую систему будем называть точечной. Выражения для нуль-изоклин точечной модели могут быть выписаны:
«1 I и и
о, (3)
Заметим, что при выполнении условия а — к1 < С у системы имеется три положения равновесия. Но лишь одно из них — нулевое — оказывается устойчивым. При любых положительных значениях параметров системы это всегда устойчивый узел. Любое возмущение по активатору, превышающее пороговое значение, будет нарастать со временем. При наличии в системе ингибитора величина порога может меняться (не всегда в сторону уменьшения). Приближенно величину порога можно оценить как
к1и0
исг ----
а — К1
(2)
Выражение (2) даёт несколько завышенную величину пороговой концентрации активатора.
Поскольку в рассматриваемой системе распространяются две последовательных автоволны — активатора и ингибитора, то формирование тромба можно трактовать как результат взаимодействия двух этих автоволн. Для краткости будем называть в дальнейшем такую феноменологическую модель двухволновой (в [26] был использован более громоздкий термин — двухавто-волновая модель).
Предсказания модели [7, 15] подтолкнули к исследованию пространственной динамики свертывания. До её появления большинство моделей свертывания ограничивалось рассмотрением только гомогенной кинетики, то есть кинетики с перемешиванием. В работе [8] была построена количественная модель свертывания, опирающаяся на известную биохимическую схему реакций, в [9] исследовалась динамика пространственного роста тромба. Добавление к этой модели гипотезы о переключении активности тромбина и последующая редукция и упрощение (подробнее см. [25] привело к системе уравнений
ди ^ . . . 1 + К2 и . .
— = DZS.il + К\пш( 1 — и)------—----и, (3)
дЪ 1 1 ' 1 + Kзv
дv
— = 1)Аг + К5и2 - К6и, дъ
дт
-р— = иАи> + и — К^ик дъ
Уравнения (5-7), записанные в безразмерном виде, отражают динамику трёх переменных. Как и ранее, первое уравнение системы описывает динамику изменения ак-
аи
V
2
V
тиватора системы свертывания — тромбина (переменная и). Из записи уравнения (3) видно, что тромбин является автоката-литической переменной. Второе уравнение описывает изменение концентрации ингибитора свертывания — протеина С (пе-
ременная у). Третья переменная модели (т) соответствует XI фактору свертывания. При этом есть перекрёстная активация первой и второй переменных модели.
Характерные значения параметров модели (5-7) приведены в табл. 2.
Таблица 2
К\ 1<2 К з К4 К5 к 6
6,85 варьируется 2,36 0,087 17,0 варьируется
Особые точки системы (5-7) определяются нолями полинома
/(и)
и
А'4
К3К5
^■(1 -А'2)+
К4
и
2 К
---+ 1
) К
В системе всегда существует нулевая особая точка (и = 0). Она устойчива при любых значениях параметров. Кроме того, у рассматриваемого полинома всегда есть один отрицательный действительный корень, не имеющий физического смысла. Соответственно количество положительных действительных корней, определяющих особые точки системы, может меняться в зависимости от параметров модели от 0 до 2. Фазовое пространство для точечной системы, соответствующей системе (5-7), есть область (и ^ 0, у ^ 0, т ^ 0). Траектории системы никогда не покидают эту область.
Удобно рассматривать поведение системы (5-7) при изменении параметра К2, определяющего скорость производства активатора (тромбина), и параметра Кб, определяющего скорость инактивации ингибитора.
На рис. 1 показано изменение характера неотрицательных особых точек модели (5-7) при изменении констант К2 и Кб.В области «0» существует единственная особая точка — нулевая. При пересечении правой границы области «0» в системе появляются ещё два положения равновесия. Линии бифуркации особой точки с меньшим значением (по и), то есть расположенной ближе к началу координат, показаны пунктиром. Эта точка остаётся неустойчивой, меняя свой характер при увеличении параметра Кб следующим образом: седло-узел, неустойчивый фокус, седло-узел. Сплошными линиями показаны линии бифуркации дальней от нуля
особой точки. В областях «1» и «2» она является неустойчивой (1 — седло-узел, 2 — неустойчивый фокус), в областях 3 и 4 — устойчивой.
Рис. 1. В зависимости от параметров модели число положительных особых точек системы варьирует от 0 до 2. В заштрихованной области «0» в системе существует единственная особая точка — нулевая. Сплошные линии показывают бифуркационные значения дальней от нуля особой точки (области: «1» — седло-узел, «2» — неустойчивый фокус. «3» — устойчивый фокус, «4» — устойчивый узел). Линия, разделяющая области «2» и «3» — линия бифуркации Пуанкаре-Анд-ронова-Хопфа. Ближняя к нулю особая точка остаётся неустойчивой в данной области параметров. Пунктирные линии показывают её бифуркационные значения (слева направо просисходит переход между областями: седло-узел, неустойчивый фокус, устойчивый фокус, узел. (По работе [26])
Линия между областями «2» и «3» — линия бифуркации Пуанкаре-Андро-нова-Хопфа. Области 3 и 4 соответствуют бистабильности, когда в системе есть две устойчивых точки, разделенных неустойчивой. Уменьшение К2 приводит к расширению области «2». Численные эксперименты показали наличие здесь каскада бифуркаций удвоения периода и переход к хаотическим колебаниям по сценарию Фейгенбаума. При выбранных значениях
остальных параметров область существования предельных циклов и хаоса существуют только в очень узкой зоне значений Кб, примыкающей к линии бифуркации Хопфа слева.
Общим свойством моделей (1-2) и (5-7) является автокатализ при образовании активатора и производство ингибитора только в области, где уже имеется активатор. Это следует из вида уравнений обеих моделей. Модели также демонстрируют пороговое поведение при увеличении начальной концентрации активатора.
Дополним модели (1-2) и (5-7) уравнением, описывающим динамику образования мономеров фибрина, полагая, что последующая полимеризация является быстрым процессом:
d'£
dt
—— = ки
для (5-7) или, для модели (1-2),
^ = кв.
dt
В первых работах применялась грубая одностадийная математическая модель полимеризации фибрина. Это объяснялось тем, что рассматривались системы in vitro, для которых несущественно влияние процессов переноса, важно только то обстоятельство, что фибрин-полимер не диффундирует. По современным представлениям, процесс полимеризации в потоке будет отличатся от процессов полимеризации в неподвижной среде, для корректного описания потоков необходимо применять более сложные модели полимеризации.
Отметим одну особенность рассматриваемых задач. Поскольку все метаболиты системы свертывания — белковые макромолекулы с примерно одинаковой молекулярной массой, то их подвижность будет примерно одинаковой. Как следствие, все метаболиты будут иметь близкие коэффициенты диффузии. В моделях (1-2) и (5-7) коэффициенты диффузии полагаются равными для всех переменных каждой модели.
Результаты численных расчётов. Основные типы решений. Сравним основные типы решений систем уравнений в частных производных в одномерном плоском случае. И для (1-2, 8), и для (5-8) численно решалась смешанная задача. Для
всех метаболитов свертывания в качестве граничных условий принимались условия отсутствия потока вещества через границы области. В начальный момент времени концентрации всех метаболитов, кроме тромбина, полагались нулевыми. Начальное возмущение по тромбину во всех задачах представляло собой ступеньку, расположенную вблизи левой границы расчётной области. О применяемых для расчётов численных методиках см. в [16, 18, 25].
(1)-(3) в Ш плоском случае. Параметры приведены в табл. 1. Область «БВ» соответствует режиму бегущих импульсов с пульсирующей амплитудой, «Е» — режиму эховолны (см. ниже), «1», «2», «3», «4» —образованию слоистого тромба с соответствующим числом структурных элементов (полос). Единицы измерения [ио] = нМ, [^о] = 10_2 нМ (по данным [30]
Сначала зафиксируем два параметра для каждой задачи, на плоскости которых удобно классифицировать различные динамические режимы. Для двухавтовол-новой модели удобно выбрать плоскость (и0,у0) —характерные значение концентраций активатора и ингибитора — постоянные, входящие в систему (1-2). Для системы с учётом переключения активности тромбина удобно использовать плоскость (К2,Кб). На приведённых ниже рисунках по расчётным данным построены бифуркационные диаграммы на плоскости соответствующих параметров. На рис. 2 виден характер чередования режимов: образование уединенной фибриновой структуры (единственный момент обострения для активатора) сменяется образованием структур с конечным числом полос. Границы этих режимов достаточно изрезаны. Выше на графике приведена область параметров, в которой генерируются эховолны [16-19]. Под эховолной здесь и далее понимается автоволна (обычно активато-
ра), которая после деления первоначальной автоволны распространяется в направлении, противоположенном первоначальному направлению распространения «материнской» автоволны. При дальнейшем возрастании параметра у0 режимы с генерацией эхо-волн сменяется областью с появлением бегущего импульса с пульсирующей амплитудой. При дальнейшем увеличении пороговой концентрации ингибитора бегущие импульсы сменяются волнами переключения. Характер и чередование режимов были описаны ранее в работах [16-18, 30].
Похожий характер чередования режимов (различных типов решений) отмечен и в модели с переключением активности тромбина.
На рис. 3 в плоскости параметров (К2Кб) сплошными линиями показаны границы, разделяющие разные режимы распространения одного и того же начального возмущения. Для сравнения с динамикой точечной модели пунктиром нанесены линии бифуркации особых точек, соответствующие рис. 1.
Рис. 3. Бифуркационная диаграмма для точечной модели на плоскости параметров (К2, Кб) для пространственно-распределенной модели (5-7). Пояснения см. в тексте. Одномерная задача рассматривалась на отрезке длиной Ь = 10. Коэффициент диффузии полагался одинаковым для всех трёх переменных модели, О = 0,00026. Активация системы была смоделирована локальным повышением значения первой переменной на левой границе (и\ = 0,5 в 30 точках, прилегающих к левой границе, шаг по пространству Н = 0,001) [26]
Все обнаруженные режимы поведения распределенной системы условно можно разделить на четыре типа (А-0 на рис. 3). Рассмотрим их поочередно.
Область А — область затухания возмущения. Начальное возмущение на левой границе отрезка формирует импульс, который двигается от активирующей границы. При этом он распространяется все медленнее по мере того, как его амплитуда уменьшается, с последующим выходом системы на нулевое пространственнооднородное решение. По характеру образующихся фибриновых структур примерно соответствует области «1» на рис. 2 для двухавтоволновой модели.
В области параметров, соответствующих А на рис. 3, увеличение значения К2 приводит к тому, что размер тромба увеличивается, а увеличение значения Кб — к более длительному росту сгустка до момента остановки. Экспериментально было показано [12], что после активации свертывания тромб некоторое время растёт с примерно постоянной скоростью, а затем наблюдается резкая остановка. Этот режим количественно хорошо описывается полной нередуцированной моделью [24] и качественно сохраняется в упрощённой модели (5-8) рис. 4.
К.
Рис. 4. Распространение автоволны с последующей остановкой, имитирующее экспериментально наблюдаемый характер роста фибринового сгустка [12]. (К2 = 15, Кб = 0,044, область А). Пространственные профили распределения тромбина (переменная и) приведены с интервалом АЬ = 10
При К2 < 5 область А включает в себя часть области «3», в которой соответствующая точечная система имеет также ненулевую устойчивую особую точку и, следовательно, система становится триггерной. Выбор начальных условий определяет, в область притяжения какой именно особой точки попадет решение системы. Для того чтобы попасть в область притяжения ненулевой особой точки и получить волну переключения, необходимо увеличить раз-
мер области начального возмущения, одновременно снижая активирующую концентрацию. Аналогичная картина зависимости от амплитуды активирующей концентрации наблюдается и в точечной модели. При некоторых параметрах область притяжения верхнего ненулевого устойчивого состояния оказывается лежащей в середине отрезка возможных начальных концентраций по первой переменной. Значения меньшие и большие, то есть лежащие вне этого выделенного коридора, приводят к тому, что мы оказываемся в области притяжения нулевого состояния системы. В этих случаях система является дважды пороговой по концентрации активатора.
Область В — бегущие импульсы. В верхней части области В (при больших значениях К2) наблюдается классический автоволновой режим — распространение незатухающего бегущего импульса. В области В пространства параметров (К2,Кб) уменьшение значения К2 приводит к возникновению бегущих импульсов с пульсирующей амплитудой и появлению режимов с генерацией эхо-волн, которые аналогичны соответствующим режимам для двухавтоволновой модели. При К2 < 5 получена сложная зависимость решения от начальных условий. В этих условиях возможен переход к химической турбулентности. Под химической турбулентностью по-
нимается сложные апериодические колебания концентрации, неоднородные по пространству [31].
Область С. Граница области В справа совпадает с линией бифуркации Хопфа при значениях К2 > 5. Справа от нее, в области С, в ответ на возмущение формируется волна переключения. Она соответствует триггерной системе без диффузии, когда соответствующая точечная система имеет две устойчивых особых точки. Волны переключения относятся к хорошо изученным автоволновым режимам [32]. Аналогичные волны переключения существуют и в двухволновой модели, но не входили в круг интересов авторов исследований [15-19] из-за того, что существуют при значениях параметров, сильно отличающихся от физиологически обоснованных.
Область D. В некотором диапазоне К2 при переходе от области затухания начальных возмущений (А) к области бегущих импульсов (В) существует область О, в которой распространяющийся из зоны активации импульс после остановки не исчезает, а стабилизируется и существует неограниченно долго. рис. 5а иллюстрирует выход системы на такое стационарное пространственно-неоднородное решение. В двухволновой системе решения такого типа отсутствуют.
'! а '■ (6}
Рис. 5. Распространение автоволны тромбина с последующей остановкой и выходом на стационарное пространственно-неоднородное решение в области D (К2 = 11,0; Кб = 0,062, остальные параметры приведены в табл. 2). Профили даны с интервалом АЬ = 10. На рис. приведена пространственно-временная дмаграмма для тромбина (вверху) и набор сечений концентрации фибрина-полимера в последовательные моменты времени, по [26]
Для фибрина (а он в обеих моделях представляет собой след от тромбина) установление подобной структуры приводит к формированию структуры (тромба), для которой концентрация фибриновых звеньев увеличивается к краю тром-
ба (рис. 5б). Такой сценарий образования тромба может иметь физиологический смысл, так как приводит к более прочным по краям фибриновым тромбам, способным выдерживать достаточно большие ме-
ханические напряжения в потоках жидкости.
Подводя краткий итог сказанному в данном разделе, отметим, что за исключением режима формирования тромба с плотным краем (установление стационарной структуры по тромбину), все остальные режимы практически совпадают с известными в двухволновой модели. Более того, и порядок чередования динамических режимов в моделях при изменении параметра в целом совпадает для разных моделей.
Динамика формирования слоистых тромбов. Наиболее интригующим экспериментальным фактом при исследовании динамики свертывания крови оказался режим формирования колец фибрина [12]. При некоторых условиях активация плазмы крови в чашке Петри стеклянным бусинами сопровождалась следующей картиной. Вокруг активирующей бусины вырастал обычный фибриновый сгусток. Затем через какое-то время на некотором расстоянии от него формировалось новое кольцо фибрина. Экспериментально был зарегистрирован режим образования трёх таких колец [12]. Оказалось, что каждая из рассматриваемых моделей способна предложить несколько сценариев образования таких слоистых тромбов. В данной работе мы рассмотрим только те режимы, которые приводят к образованию сгустков с конечным числом полос, имитируя описанные выше экспериментальные данные. Отметим, что в обеих моделях также был получен ряд сценариев образования слоистых сгустков, формирование которых сопровождалось заполнением всего доступного пространства, то есть когда начальное возмущение распространялось до естественных границ среды.
Рассмотрим результаты расчётов для двумерной задачи в области параметров, которая в одномерном случае соответствует остановке бегущего импульса активатора на конечном расстоянии от места активации и образованию устойчивой структуры тромбина [8, 25, 26].
Одиночная стационарная структура тромбина численно исследовалась на устойчивость к малым возмущениям начальных данных в двумерном случае. В квадратной области задавалась ступенчатая функция, отвечающая начальному распределению тромбина. Проверялось, что в отсутствие возмущений
начальных данных процесс развивается так же, как и в одномерном случае. Затем вносилось малое возмущение (и(0,х,у) = и(0,х) ■ (1 + А сов(2п1у/Ь)), здесь первый сомножитель соответствует начальному условию для одномерной задачи, Ь — ширина области, А — амплитуда, в дальнейшем её удобно обозначать в процентах, I — целочисленный параметр) в начальные условия.
Возмущение с амплитудой меньше некоторого порогового значения не приводит к качественному изменению решения, формируется бегущий импульс тромбина, который затем останавливается и существует неограниченно долго. При превышении порогового значения система стремится к тривиальному решению. Расчёты проводились для значений параметров К2 = 12, Кб = 0, 058. Для высокочастотных возмущений (большое значение I) про-цее развивался как одномерный.
Формирование слоистого тромба с конечным числом полос в модели с переключением активности тромбина (5—8). Основу формирования слоистых сгустков в данном случае составляет режим движения автоволны с последующей остановкой. В процессе движения автоволны бегущий импульс активатора делится на два новых, как и в случае распространения эховолн [19]. Первый распространяется вперед, второй — назад. После этого импульсы останавливаются, как при формировании тромба с плотным краем. Образуется стационарная структура по тромбину с двумя максимумами. Рост фибринового тромба продолжается бесконечно. Динамика процесса показана на рис. 6.
При некоторых параметрах движущийся вперед импульс может дополнительно разделиться и тогда наблюдается установление трёх неподвижных пиков тромбина и соответственно образование трёх слоёв фибрина (рис. 7). При этом наблюдается сильная зависимость от начальных условий. Если начальное возмущение захватывает по пространственной координате в 1,5 раза большую область при той же амплитуде, то при данных значениях параметров мы получим два установившихся импульса, аналогичных изображённым на рис. 6. Теоретическое и численное исследование одномерного структурообразования в данной модели свертывания приведено в [33].
юоо
750
500
250
200 ■
100 ■
(а>
Рис. 6. Режим образования слоистого тромба с числом полос, равным двум (К = 8,2; Кб = 0,0742, на границе областей В и D). Пространственные профили приведены с интервалом Аі = 20 [26]
(б>
Рис. 7. Режим образования слоистого тромба с числом полос, равным трём (К = 8,0; Кб = 0,0745, на границе областей В и D). Пространственные профили приведены с интервалом по Аі = 40
Для системы с переключением активности тромбина также исследовалась численно устойчивость к малым возмущениям начальных данных. К начальной ступеньке добавлялись начальные возмущения, как описано выше.
Решение задачи практически несущественно к высокочастотным возмущениям (с малой длиной волны), поэтому все описанные ниже эффекты отмечены при числе длин волн возмущения (I) от 2 до 4. Задача решалась численно в квадратной области.
В двумерном случае при возмущении в 10% формируется одна стационарная структура тромбина (в одномерной задаче при тех же значениях параметров — две структуры, параметры расчёта соответствуют параметрам, при которых проведены расчёты, представленные на рис. 6). Но при данных параметрах характер решения системы существенно меняется при увеличении амплитуды возмущения.
Эти двумерные расчёты характеризовались отсутствием тромбина за передним фронтом бегущего импульса. Но при уве-
личении амплитуды начального возмущения появляется область с ненулевой концентрацией тромбина за передним фронтом. Концентрация тромбина имеет апериодические колебания и по пространственным переменным, и по времени. При времени Т = 460 теряется начальная симметрия решения (рис. 8), что затем приводит к разрыву переднего фронта (рис. 9).
В ходе дальнейшего расчёта распределение концентрации тромбина становится хаотическим, что, по всей видимости, соответствует режиму химической турбулентности. Любопытно, что те же значения параметров в одномерном случае, как отмечено выше, приводят к формированию тромба с числом структурных элементов фибрина, равным двум. Подобное хаотическое поведение наблюдалось и при амплитуде возмущения 50%, за исключением того, что передний фронт оставался непрерывным. При использовании других численных методов также была отмечена потеря симметрии решения. Более подробное описание численного исследования устойчивости таких решений приведено в [34].
FР
Рис. 8. Т = 460. Распределение тромбина (К2 = 8,2;
Кб = 0,0742, на границе областей В и D). Становится заметной асимметрия решения
Рис. 9. Т = 500. Распределение тромбина (К2 = 8,2;
Кб = 0,0742, на границе областей В и D). Разрыв переднего фронта
Формирование слоистого тромба с конечным числом полос в феноменологической модели (1-2, 8). В двухавтоволновой модели обнаружены два варианта образования структуры с двумя полосами фибрина. Первый (типовой) сценарий формирования такого тромба включает формирование бегущего импульса, два момента роста его амплитуды в режиме с обострением (в англоязычной литературе для обо-
значнения режимов с обострением употребляется термин blow up, поэтому ниже будет использоваться прямой перевод — вспышка) и достижение подпороговых значений концентрации активатора после второй вспышки в режиме с обострением (рис. 1О). По такому же сценарию формируются и тромбы, состоящие из большего числа полос.
Другой — резонансный — сценарий впервые был обнаружен при теоретическом исследовании двухавтоволновой модели и затем зафиксирован в расчётах [20]. Динамика такого формирования тромба приведена на рис. 11. После второй «вспышки» концентрации тромбина образуются две неоднородности по тромбину в виде локальных максимумов, распространяющихся в противоположных направле-
X
ниях. При этом тот локальный максимум концентрации тромбина, который движется к месту первой «вспышки», то есть туда, где ещё осталось небольшое количество ингибитора, вновь начинает расти. Возмущение, распространяющееся в область пространства, где ингибитора нет, затухает. В результате третья «вспышка» концентраций метаболитов происходит на месте первой или с небольшим сдвигом.
Рис. 10. Образование слоистого тромба с числом полос равным двум в феноменологической модели (1-2, 8). О = 610-4. Размер расчётной области 1 мм, шаг сетки Н = 0,01 мм, ширина начального возмущения по тромбину 0,18 мм, амплитуда 0,16 нМ. Остальные параметры приведены в табл. 1. Пространственные профили приведены с интервалом Аі = 40 (по [26])
Такой сценарий взаимодействия активатора и ингибитора может повторяться несколько раз. В таких режимах все пространство делится на две области, где возмущения дают «вспышки» (развиваются в режиме с обострением) и где они «дисси-пируют». Причём повторные «вспышки» происходят в той области, где от их предшественников осталось небольшое количество ингибитора. В расчётах зафиксировано до четырёх последовательных вспышек концентрации активатора в режиме с обострением. Этот режим приводит к формированию полос с гораздо большей плотностью фибриновых звеньев, чем в стандартном сценарии.
В двумерном случае структуры практически при всех значениях параметров имеют повторяющийся характер. Результирующая сформированная фибринововая структура представляет собой след от распространения по плоскости концентрационных подвижных «пятен» активатора. Подробнее двумерные решения в математической модели описаны в [17, 18, 30]. Отметим, что и данная задача в двумерном случае оказывается устойчивой по отношению к высокочастотным возмущениям.
Спиральные волны в модели с переключением активности тромбина.
В системах полулинейных уравнений реакция-диффузия в двумерном случае часто существуют решения в виде спиральных волн [32, 34]. В системе (5-7) такие решения найдены численно при значениях параметров К2 = 17, 0, Кб = 0, 04.
Для задачи с полным перемешиванием эти параметры соответствуют случаю существования в системе ОДУ неустойчивого фокуса и устойчивого предельного цикла, а в одномерной системе с диффузией — устойчивого бегущего импульса. Решения типа бегущего импульса относятся к грубым, они существуют в достаточно широком диапазоне параметров модели [26].
Для численных экспериментов со спиральными волнами в качестве начальных условий в одной половине расчётной области был взят профиль концентраций, соответствующий установившейся форме бегущего импульса. Во второй половине области концентрации всех веществ были приняты нулевыми.
В результате решения получалась устойчивая спиральная волна, центр которой двигался по окружности. Для реше-
ний типа спиральных волн использовались условия отсутствия потоков через границы области, соответствующие стенкам, нулевые концентрации всех веществ во входном сечении и свободные (неотражающие)
(3}
условия в выходном сечении. Взаимодействие этой волны с границами расчётной области не приводило к нарушению формы спиральной волны (рис. 12).
(6}
Рис. 11. «Резонансный» режим образования слоистого тромба с числом полос, равным двум в феноменологической модели (1-2, 8). ио = 3,03 нМ. Остальные параметры такие же, как для режима, представленного на рис. 10
Рис. 12. Решение в виде спиральной волны в ограниченной области. Концентрации активатора, значения параметров К2 = 17,0, Кб = 0,04. а) время і = 205. Первый виток спирали достигает границы расчётной области. б) время і = 400. Спиральную волну можно продолжить за границы области, её форма не меняется при взаимодействии со стенкой
Существенно иными свойствами обладает решение системы уравнений ре-акция-диффузия-конвекция. Результаты для системы с конвекцией описаны в [35].
Сопряженная задача. Система уравнений. Введём в рассмотрение наряду с задачей (5-7) сопряженную систему уравнений. Решение её позволяет выделить класс возмущений, к которым наиболее чувствительна основная задача. Поведение нормы решения сопряженной задачи позволяет судить об устойчивости решения основной задачи к данному классу возмущений (в норме Ь2) [36, 37].
С помощью сопряженной системы исследовалась устойчивость стационарных решений в одномерном случае [25], но в [25] сопряженная система строилась но основе формулы Лагранжа без использования формулы Тейлора для линеаризованной системы уравнений. Для нелинейной задачи сопряженная система может быть построена неединственным образом. Для исследования устойчивости получающихся стационарных решений задачи (5-7) сопряженная система нелинейных уравнений была построена с помощью формулы Тейлора [36].
Введём пространства функций X, У, такие что X = У = Ь2(О), где О = [0,ЬХ] х [0,Ьу] — расчётная область задачи. Пусть и Е О, V Е О, введём скалярное произведение
ЬУ Ьх
{и,у)ь2 =
щуі йхйу,
0 0
і=1
а соответствующая норма имеет вид
\\и\\ь2 = (и,и)1/2.
Обозначим: и1 = и, и2 = ,ш, и3 = V и
рассмотрим нелинейный оператор F, Г :
X -► У:
Г =
^ — БАщ — щ + Л4И2
“Iа — -ОАиз — + Л6«з
и1
\
В основу определения сопряженного оператора положено равенство [36]
потока через границы области) для оператора, сопряженного к Г'л(и), интегрированием по частям получим
С^М)* = аіа§ •
Для нелинейной части оператора Г справедливо
(Г'' (и)У,и*)
Г'^^ и*йО, =
Гц)* иїст=(у,{Ц(и))* и*)
откуда (Г(и))* = (Г(и))Т — транспонированная матрица Якоби системы (5-7). Таким образом, для сопряженного оператора получим
(А(и))*
(Г'(іи))* йі
Г (и)
Г'(іи)йіи = А(и)и.
Оператор (А(и))* = у § Г '(іи)йі^ будем
называть сопряженным оператором, соответствующим нелинейному оператору F.
Для оператора F справедливо представление
Г (и) = Гл(и) + Гі (и),
где р — линейный оператор, соответствующий дифференциальной части оператора Г (и), а Г ^ соответствует нелинейным функциям в правой части исходной системы уравнений.
Найдём сопряженный оператор (Г'(и))*. Для этого воспользуемся формулой Лагранжа (и* — произвольная функция, принадлежащая рассматриваемому пространству):
(Г'(и),ь,и*)ь2 = (ь,(Г'(и))*и*)ь2.
Для производной оператора Г получим равенство
Г'(и) = Г'„ (и) + Г' (и).
Поскольку X = X * = У = У * = Ь2(0), то с учётом граничных условий (отсутствие
(Г'л(іи))* йі+
Г'(іи)) * йі.
Сопряженная система будет иметь вид
с)'и
—7ГГ = £>Дг’1 - VI + 1’2 + К5щьз + У\д\,
дЬ
с)'и
---7^ = ОАу2 - К4У2 + ^102,
дЬ
~~^~ = ~ Аб1'з + ад,
где функции д1, д2 и д3 представляются интегралами
91 =
' К1и2а(1 — 2и1а + 2К2и1а — 3К2и^а) 1 + К3и3а
йа,
92 =
К1 ща(1 — и1а)(1 + К2 и1а) 1 + К3и3а
йа.
9з =
' К1К3и1и2а2(1 — и1а)(1 + К2и1а) (1 + К3и3а)2
йа.
Приведённые интегралы могут быть вычислены в явном виде. Поскольку получающиеся выражения достаточно громоздки, они здесь опущены.
3
V
і
1
1
0
1
1
0
0
0
+
1
1
и
1
0
и
Решение сопряженной задачи находится в обратном времени с использованием вычисленного решения основной задачи. Рассмотрение сопряженной системы позволяет выделить класс возмущений, которые приводят к наибольшим изменениям решения в финальный момент времени. Рассмотрим v(x,y) — дважды непрерывно дифференцируемую функцию, удовлетворяющую на границе расчётной области условию = 0. Очевидно, её можно представить в виде v(x,y) = v\\(x,y) + v±_(х,у), где
(у(х,у),и(х.у,Т)) , т)
(и(х.у,Т)Ах,У.Т))4 У' >' v±(x,y) = v(x,y) - Щ(х,у)-
В силу линейности сопряженной задачи решение сопряженной системы с начальными условиями v(x,y) представимо в ви-
де v(x,y,t) = v\\(x,y,t) + v±(x,y,t) для всех t < T.
Если начальное возмущение есть любая функция из v±_(x,y), то решение сопряженной задачи в любой момент времени будет ортогонально решению основной задачи. Умножим на v и и системы (1-3) и (6-8) соответственно. После интегрирования от 0 до T с учётом сопряжения и граничных условий получим
(vo,Uo) = (vt±,u(T )) = 0.
Аналогичное равенство устанавливается и при интегрировании от t до T с учётом сопряжения и граничных условий. Рассмотрим функционал
т. . \\u(x,y ,t)\\2 1 . . . . ..
J(t) =--------------= ~(u(x,y,t),u(x,y,t))
и его первую вариацию 8J = (u,8u(t)).
0— '—40
Рис. 13. Распределение первой переменной начального возмущения G,3
сопряженной задачи при t = 0 для амплитуды
Она всегда равна 0 на функциях из множества v±(x,y) = v(x,y) — v\\ ^,у). Максимальная вариация при Т получается для функции гф,у) =
но в силу линейности сопряженной задачи достаточно рассмотреть лишь
Щ ^,у) = u(x,У,T).
Тогда v\\(x,y,0) — возмущение, к которому наиболее чувствительно решение основной задачи в момент времени Т в норме Ь2.
При значениях параметров К2 = 8,0, Кб = 0,0745 поведение решения основной задачи существенно зависит не только от амплитуды начального распределения, но и от размера области активации [19].
Рассмотрим теперь двумерные решения. Как и прежде, на начальное условие в виде ступеньки накладыва-
о—ь
о
ется малое гармоническое возмущение u(0,x,y) = и(0^) ■ (1 + А со$,(2п1у/Ь)). Рассмотрим решение сопряженной задачи для тех же условий, для которых получен результат, представленный на рис. 8, 9. Норма возмущения в линеаризованной задаче убывает со временем, что свидетельствует об устойчивости решения основной задачи к данному классу возмущений. Этот результат согласуется с результатом, полученным в [25].
На рис. 13 показано конечное распределение v1 (для времени Ь = 0) в случае нестационарного (хаотического) решения. Динамика нормы решения сопряженной задачи также позволяет говорить о том, что решение основной задачи устойчиво по отношению к малым возмущениям, определяемыми решениями сопряженной задачи.
Ше
Рис. 14. Финальное распределение первой переменной сопряженной задачи для амплитуды возмущения начального распределения основной задачи 0,5
В случае 50% возмущения начального распределения тромбина при решении основной задачи (расчётная область и ширина начального распределения тромбина остались теми же) формируются две стационарные структуры. На рис. 14 изображено финальное распределение — пер-
вой компоненты решения сопряженной задачи. В этом случае возмущение начальных данных, приводящее к наибольшему изменению решения основной задачи, качественно отличается от распределения, показанного на рис. 13.
Иные результаты были получены при увеличении ширины области активации до
0,26, к =1, Р = 50%. Размеры расчётной области Ьх = 3, Ьу = 1. При данных начальных условиях в основной задаче снова наблюдается хаотическое поведение. Решение сопряженной задачи практически совпадает с решением, представленным на рис. 13.
В силу линейности сопряженной системы её решение после перемасштабиро-вания использовалось в качестве возмущения начальных условий основной задачи (именно к полученному классу начальных возмущений решение основной задачи должно быть наиболее чувствительно). При этом динамика решения основной задачи существенно меняется. В динамике тромбина вместо апериодических колебаний появляется упорядоченность. Одно из промежуточных распределений тромбина приведено на рис. 15. Любопытно формирование круговой структуры тромбина между двумя плоскими фронтами.
100 8060
01—-и—I—и—,---------------
О 50 100
Рис. 15. Распределение тромбина (в изолиниях) при расчёте основной задачи с начальными данными, полученными в результате наложения решения сопряженной задачи
Обсуждение и выводы. Рассмотреные математические модели пространственной
динамики свертывания крови отличаются степенью детализации процессов и сделанными при построении модели допущениями. В основу первой феноменологической двухволновой модели положена гипотеза о существовании автокаталитического ингибитора. Отметим, что эта модель была исторически первой и сыграла большую роль в исследовании характера роста фибринового сгустка [7, 15-22]. Другая модель привлекает гипотезу о переключении активности тромбина по отношению к его субстратам [25].
Теоретические предсказания феноменологической модели в своё время стимулировали экспериментальное исследование пространственной динамики свертывания крови [12]. Можно выделить несколько особенностей формирования фибринового сгустка при активации плазмы крови, помещённой в чашку Петри. Такой экспериментальной постановке соответствует система типа реакция-диффузия. Первым важным моментом является характер роста «нормального» фибринового сгустка. После небольшого запаздывания, необходимого для запуска каскада биохимических реакций свертывания, наблюдается рост фибринового сгустка с примерно постоянной скоростью в течение некоторого времени с последующей резкой остановкой этого роста [12]. Модель с переключением активности тромбина описывает такой режим (рис. 4).
Другим интересным фактом является возможность образования в однородной среде структур (полос или колец) фибрина. Интересно то, что этот режим был предсказан феноменологической моделью как результат распространения автоволны тромбина с пульсирующей амплитудой [7, 15, 18] и уже затем обнаружен в эксперименте [12].
Большинство динамических режимов свертывания описываются каждой из описанных моделей. В феноменологической (двухавтоволновой) модели не зафиксирован лишь сценарий формирования фибринового тромба с плотным краем. Данный динамический режим в системах типа реакция-диффузия был зарегистрирован впервые именно при исследовании математической модели с переключением активности тромбина и представляет особый интерес с точки зрения теории самоорганиза-
ции. Возможность создания неподвижного импульса на некотором удалении от активирующей границы может иметь биологическое значение и в применении к другим процессам, например, процессу дифферен-цировки в ходе морфогенеза.
В математической модели с переключением активности тромбина обнаружены как традиционные для систем реакция-диффузия решения типа бегущих импульсов, спиральных волн, так и решения с остановкой автоволны тромбина на конечном расстоянии от места активации. В одномерном случае устойчивость этих режимов исследована в [25]. Результат оказался несколько неожиданным — остановившийся импульс оказался устойчивым, хотя в точечной системе область параметров соответствовала хаотическим колебаниям (странному аттрактору). Проведённое исследование с использованием аппарата сопряженных уравнений показало, что в случае нескольких пространственных измерений в такой системе действительно наблюдаются хаотические режимы. При этом необходимо низкочастотное надпороговое возмущение системы в направлении, перпендикулярном направлению распространения волны. Вызвать пространственно-неоднородные хаотические колебания может и наличие в системе сдвигового течения. Интересно, что решения систем уравнений обеих математических моделей оказываются чувстви-телльными лишь к низкочастотным возмущениям.
Спиральные волны, соответствующие в точечной системе, являются грубыми. Наиболее чувствительны они к начальным возмущениям в окрестности центра спиральной волны. При этом в решении задачи вместо одной спиральной волны появляется другая спиральная волна (со смещённым центром спирали). Это соответствует численным результатам по поведению спиральных волн в других средах (например, при моделировании электрической активности сердца). Для задач такого моделирования существенным является то, что спиральную волну крайне трудно «разрушить» [38, 39].
Даже простейшие математические модели пространственной динамики свертывания крови демонстрируют весьма сложное поведение. Среди решений уравне-
ний модели существуют как традиционные для систем реакция-диффузия волны переключения, бегущие импульсы, спиральные волны и хаотические режимы (в многомерных случаях), так и решения, не зафиксированные ранее в моделях реакционно-диффузионного типа. Кроме того, характер развития неустойчивости для большинства таких моделей может определяться не только нормой возмущения, но и пространственным распределением такого возмущения. Иными словами, на одинаковые по величине превышения порога воздействия система может реагировать совершенно по-разному в зависимости от вида возмущения. Большое количество разнообразных неустойчивостей для такой системы является фактором, значительно усложняющим анализ математических моделей с учётом течений крови и фазовых переходов при образовании фибриновых тромбов.
Благодарности. Авторы выражают благодарность Г.Т. Гурии и Ф.И. Атаулло-ханову за то, что они вызвали интерес к проблематике свертывания, за постановки задач и плодотворные обсуждения. Многие рисунки в данной статье базируются на расчётах, проведённых в своё время Т. К. Старожиловой и В.И. Зарницыной, которым авторы также выражают свою благодарность. Мы признательны большому числу коллег из Гематологического научного центра РАН, Биологического факультета МГУ им. М. В. Ломоносова и МФТИ за плодотворные дискуссии и внимание к работе. Особую благодарность за внимание к работе авторы хотят выразить А.С. Холодову и А. А. Шананину.
Литература
1. Levin S.N. Enzyme amplifier kinetics // Science. — 1966. — V. 152 (3722). -P. 651-653.
2. Macfarlane R.G. An enzyme cascade in the blood clotting mechanism, and its function as a biochemical amplifier // Nature. — 1964. — V. 202(4931). -P 498-499.
3. Khanin M.A., Semenov V.V. A mathematical model of the kinetics of blood coagulation // J. Theor. Biol. — 1989. — V. 136. — P. 127-134.
4. Willems G.M., Lindhout T.,
Coenraad H., Hermens W.T., Hemker H.C. Simulation model for thrombin generation in plasma // Haemostasis. — 1991. — V. 21. — P. 197-207.
5. Kessels H., Willems G.M.,
Hemker H. C. Analysis of trombin generation in plasma // Comput. Biol. Med. — 1994. — V. 24. — P. 277-288.
6. Jones C.J., Mann K.G. A model for the tissue factor pathway to thrombin. II. A mathematical simulation // J. Biol. Chem. — 1994. — V. 269, № 37. — P. 23367-23373.
7. Атауллаханов Ф.И., Гурия Г.Т. Пространственные аспекты динамики свертывания крови. I. Гипотеза // Биофизика. — 1994. — Т. 39, № 1. — С. 87-96.
8. Zarnitsina V.I., Pokhilko A.V., Ataullakhanov F.I. A mathematical model for the spatio-temporal dynamics of intrinsic pathway of blood coagulation. I. The model description / / Thrombosis Research.
1996. — V. 84, № 4. — P 225-236.
9. Zarnitsina V.I., Pokhilko A.V., Ataullakhanov F.I. A mathematical model for the spatio-temporal dynamics of intrinsic pathway of blood coagulation. II. Results // Thrombosis Research. — 1996. — V. 84, № 5. — P. 333-344.
10. Барынин Ю.А., Старков И.А., Ха-нин М.А. Математические модели физиологии гемостаза // Изв. АН. Серия биологическая. — 1999. — № 1. — С. 59-66.
11. Атауллаханов Ф.И., Волкова Р.И., Гурия Г.Т., Сарбаш В.М. Пространственные аспекты динамики свертывания крови. III Рост тромба in vitro // Биофизика. — 1994. — Т. 40, № 6. — С. 1320-1328.
12. Ataullakhanov F.I., Guria G.T., Sarbash V.I., Volkova R.I. Spatio-temporal dynamics of clotting and pattern formation in human blood // Biochimica et Biophysica Acta. — 1998. — V. 1425. — P. 453-468.
13. Sinauridze E.I., Volkova R.I., Krasotkina Yu.V., Sarbash V.I., Ataullakhanov F.I. Dynamics of clot growth induced by thrombin diffusion into nonstirred citrate human plasma // Biochimica et Biophysica Acta. — 1998. — V. 1425. — P. 607-616.
14. Атауллаханов Ф.И., Зарницы-на В.И., Кондратович А.Ю., Лобанова Е.С., Сарбаш В.И. Особый класс автоволн (автоволны с остановкой) определяет пространственную динамику свертыва-
ния крови // Успехи физических наук. — 2002. — Т. 172, № 6. — С. 671-690.
15. Атауллаханов Ф.И., Гурия Г.Т., Сафрошкина А.Ю. Пространственные аспекты динамики свертывания крови. II. Феноменологическая модель // Биофизика. — 1994. — Т. 39, № 1. — С. 97-104.
16. Старожилова Т.К., Лобанов А.И., Гурия Г. Т. Численное исследование образования двумерных структур в модели возбудимой среды с активным восстановлением // Математическое моделирование. — 1997. — Т. 9, № 2. — С. 21-24.
17. Лобанов А.И., Старожилова Т.К., Гурия Г. Т. Численное исследование струк-турообразования при свертывании крови // Математическое моделирование. — 1997. — Т. 9, № 8. — С. 83-95.
18. Гурия Г.Т., Лобанов А.И., Ста-
рожилова Т.К. Формирование аксиальносимметричных структур в возбудимых средах с активным восстановлением // Биофизика. — 1998. — Т. 43, № 3.
С. 526-534.
19. Лобанов А.И., Старожилова Т.К. Качественное исследование начального этапа формирования неравновесных структур в модели типа «реакция-диффузия // Математическое моделирование. —
1997. — Т. 9, № 12. — С. 3-15.
20. Лобанов А.И., Старожилова Т.К., Черняев А.П. Резонансные явления в системах типа «реакция-диффузия // Математическое моделирование. — 1999. — Т. 11, № 7. — С. 75-82.
21. Чуличков А.Л., Николаев А.В., Лобанов А.И., Гурия Г.Т. Пороговая активация свертывания крови и рост тромба в условиях кровотока // Математическое моделирование. — 2000. — Т. 12, № 3. — С. 75-96.
22. Гузеватых А.П., Лобанов А.И., Гурия Г. Т. Математическое моделирование активации внутрисосудистого тромбо-образования вследствие развития стеноза // Математическое моделирование. — 2000. — Т. 12, № 4. — С. 39-60.
23. Berg D.T., Wileyand M.R., Grinnell B.W. Enhanced protein C activation and inhibition of fibrinogen cleavage by a thrombin modulator // Science. — 1996. — V. 273. — P. 1389-1391.
24. Dang Q.D., Vindigni A., Di Cera E. Identification of residues linked to the
slow-fast transition of thrombin // Proc. Natl. Acad. Sci. USA. — 1995. — V. 92. — P. 11185-11189.
25. Zarnitsina V.I., Ataullakhanov F.I., Lobanov A.I., Morozova O.L. Dynamics of stationary spatial-nonuniform patterning in the model of blood coagulation // Chaos: Interdisciplinary Journal of Nonlinear Sciences. — 2001. — V. 11, № 1. — P 57-70.
26. Лобанов А.И., Старожилова Т.К., Зарницына В.И., Атауллаханов Ф.И. Сравнение двух математических моделей для описания пространственной динамики процесса свертывания крови // Математическое моделирование. — 2002. — Т. 15, № 1. — С. 14-27.
27. Kuharsky A.L., Fogel A.L. Surfase-mediated control of blood coagulation: the role of binding site densities and platelet deposition // Biophysical J. — 2001. — V. 80, № 3. — P. 1050-1074.
28. Буравцев В.Н., Лобанов А.И., Украинец А.В. Математическая модель роста тромбоцитарного тромба // Математическое моделирование. — 2009. — Т. 21, № 3. — С. 109-119.
29. Anand M., Rajagopal K., Rajagopal K.R. A model for the formation and lysis of blood clots // Phatophisiol. Haemostasis and trombosis. — 2005. — V. 34. — P. 110-120.
30. Старожилова Т.К. Математическое моделирование автоволновых процессов и диссипативных структур в биологических системах: дис. к. ф.-м. н. — М., 2000. — 195 с.
31. Oono Y., Kohоmoto M. Discrete model of chemical turbulence // Physical
Review Letters. — 1985. — V. 55(27). — P. 2927-2931.
32. Васильев В.А., Романовский Ю.М., Яхно В.Г. Автоволновые процессы. — М.: Наука, 1987. — 240 с.
33. Lobanova E.S.Shnol E.E, Ataullakhanov F.I. Complex dynamics of the formation of spatially localized standing structures in the vicinity of saddle-node bifurcations of waves in the reaction-diffusion model of blood clotting // Physical review E 70, 032903 (2004).
34. Лоскутов А.Ю., Михайлов А.С. Введение в синергетику. — М.: Наука, 1990. — С. 44-62.
35. Куриленко И.А., Лобанов А.И., Украинец А.В. Численное исследование устойчивости некоторых автоволновых решений в математической модели свертывания крови // Биофизика. — 2009. — Т. 54, № 1. — С. 68-76.
36. Марчук Г.И., Агошков В.И., Шу-тяев В.П. Сопряженные уравнения и методы возмущений в нелинейных задачах математической физики. — М.: Физмат-лит, 1993.
37. Марчук Г.И. Методы вычислительной математики. — М.: Наука, 1989. — С. 443-450.
38. Лоскутов А.Ю., Михайлов А.С. Основы теории сложных систем. — Ижевск: РХД, 2007. — 612 с.
39. Алиев P. Р. Моделирование электрической активности сердца на компьютере // Медицина в зеркале информатики. — М.: Наука, 2008. — С. 81-100.
Поступила в редакцию 15.09.2009.
