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

Размещено на http://www.allbest.ru/

ВВЕДЕНИЕ

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

В работе рассматриваются 3 модели, построенные на основе различных физических, биологических и химических законов.

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

Другой подход реализован в моделях Шумакова, Иткина и построенной на их основе модели Солодянникова [2]. Как пишут об этой модели авторы, её главная особенность в том, что она позволяет изучать нелинейные колебательные (в частности, периодические) процессы в кровеносной системе. Модель является самонастраивающейся. С механической точки зрения, система кровообращения в модели Солодянникова представляет собой сложную гидродинамическую систему, включающую сердце, разветвленную сеть труб и резервуаров - артериальных, венозных сосудов, капиллярных сосудов, в которых происходит передача транспортируемых кровью веществ органам и тканям. Математической идеализацией такого объекта является динамическая система дифференциальных уравнений.

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

Эта модель включает в себя 55 основных артерий тела человека, характеризующихся собственными параметрами, такими как длина, поперечное сечение, удаленность от сердца и эластичность стенок.

Основными задачами данной работы являлись:

1) получение систем уравнений моделей Карааслана и Солодянникова,

2) исследование существования и единственности решений этих систем, их устойчивости. 3) реализация моделей Карааслана и Солодянникова для проведения численных расчетов.

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

5) поиск возможностей объединения моделей Карааслана и Солодянникова с гидродинамической моделью с целью получения комплексной модели сердечно-сосудистой системы.

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

1. МОДЕЛЬ КАРААСЛАНА

1.1 Модели предшественники

Модель Карааслана [1] и её предшественники - модели Гайтона (1972 г.) Утамшинга (1985 г.) [5] и Колемана-Холла (1992 г.) описывают долговременные физиологические процессы, регулирующие артериальное давление. В них детально описана роль почки и почечных гормонов.

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

Математическая модель Утамшинга дает достаточно детальное представление о долговременной динамике почечной системы, а также включает базовую сердечно-сосудистую динамику. Новизна этой модели заключалась на тот момент в введении зависимости секреции ренина от тока ионов натрия в Macula densa.

В модели Колемана-Холла более детально учтен вклад почки (секреция ренина, тубуло-гломерулярные процессы и механизмы, а также влияние на них ангиотензина) с фокусомна токе натрия через Macula densa.



Рис. 1.1 - Сравнение почечных блоков моделей Гайтона (слева) и Колемана-Холла (справа)

На схеме ниже представлено сравнение почечных блоков в моделях Гайтона (Рис. 1.1. слева) и Колемана-Холла (Рис. 1.1. справа).

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

Таблица 1.1 - Сравнение модели Карааслана с моделями-предшественниками

Факторы, учитываемые моделью	Модель Гайтона (1972г.)	Модель Утамшинга (1985г.)	Модель Колемана-Холла (1992г.)	Модель Карааслана (2005г.)
Факторы, регулирующие секрецию ренина	--	Ток ионов натрия через Macula densa	Ток ионов натрия через Macula densa	Нерв в районе Macula densa. Почечная симпатическая активность
Факторы, регулирующие выработку ангиотензина	Концентрация ионов натрия. Почечный кровоток	Концентрация ренина	Концентрация ренина	Концентрация ренина
Факторы, регулирующие секрецию альдостерона	Концентрация ионов калия. Концентрация ионов натрия. Концентрация ангиотензина. Артериальное давление	Концентрация ангиотензина	Концентрация ангиотензина	Концентрация ионов калия. Концентрация ионов натрия. Концентрация ангиотензина. Артериальное давление
Факторы, регулирующие секрецию антидиуретического гормона	Давление крови в правом предсердии. Концентрация ионов натрия. Активность автономной нервной системы.	Концентрация ионов натрия. Избыток внеклеточной жидкости	--	Давление крови в правом предсердии. Концентрация ионов натрия. Активность автономной нервной системы
Факторы, регулирующие секрецию предсердного натриуретического пептида	Натриуретический фактор	--	--	Давление крови в правом предсердии
Параметры, зависящие от почечной симпатической нервной активности	Сопротивление афферентной артериолы	--	--	Сопротивление афферентной артериолы. Скорость секреции ренина. Скорость реабсорбции натрия в канальцах почки

1.2 Схема блоков модели Карааслана

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



Рис. 1.2 - Схема блоков модели Карааслана. Номера блоков указаны в круглых скобках

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

Для описания механизмов регуляции концентраций гормонов используется подход, предложенный Гайтоном. Этот подход основан на упрощенной модели, связывающей концентрацию гормона с уровнем его секреции. Уровень секреции задается стимулами к секреции, специфичными для каждого гормона. Скорость разрушения пропорциональна концентрации гормона в данный момент времени. Скорость изменения концентрации гормона есть разница между скоростью секреции и скоростью разрушения. Предположим, что в стационарном состоянии концентрация гормона равна С0, тогда концентрацию в любой момент времени можно представить как С0*x(t), где x(t) - относительная концентрация гормона. Учитывая, что в стационарном состоянии скорость секреции и разрушения равны, относительный уровень секреции xs можно представить в виде vs/(a*C0), где vs - скорость секреции. Тогда скорость изменения концентрации принимает вид:

Здесь - постоянная, соответствующая времени, за которое концентрация достигнет примерно 63% от стационарного уровня.

Всего в модели насчитывается 61 алгебраическое и 8 дифференциальных уравнений. Также к работе прилагается список 85 параметров и констант с их начальными значениями (Табл. П.1.) Приведём ниже блоки, в которых присутствуют дифференциальные уравнения в интегральной форме, задающие динамику процесса.

Блок 14: Объем внеклеточной жидкости

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

Блок 20: Обеспеченность кровеносными сосудами (васкулярность)

Когда в сердечно-сосудистой системе долгое время поддерживается высокий уровень кровотока, количество и диаметр кровеносных сосудов медленно увеличивается. И наоборот, если кровоток понижен, васкулярность ткани уменьшается. Математически васкулярность представляется как временной интеграл результирующей скорости прироста васкулярности, которая является разностью между скоростью формирования и разрушения сосудов. Скорость разрушения пропорциональна васкулярности. Скорость формирования зависит только от величины минутного объёма.

Блок 25: Эффект множителя автономной нервной системы (хеморецепторы + барорецепторы)

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

Блок 26: Концентрация антидиуретического гормона (вазопрессина)

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



Блок 30: Суммарное количество натрия

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

Блок 32: Концентрация ренина

Секреция ренина осуществляется гранулярными клетками юкстагломерулярного аппарата. Стимулами к секреции являются: уменьшение потока натрия через macula densa, увеличение активности почечного симпатического нерва и уменьшение давления крови вблизи юкстагломерулярного аппарата.

Блок 34: Концентрация альдостерона

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

1.3 Система дифференциальных уравнений модели Карааслана

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

x1 = Vecf(t), x2 = vas(t), x3 = Nadh(t), x4 = abaro(t),

x5 = дra(t), x6 = Msod(t), x7 = Cr(t), x8 = Nal(t).

В результате, получим систему из восьми дифференциальных уравнений:

Здесь f1 = Цu, f2 = vasf, f 3 = Nadhs, f 4 = aauto, f 5 = Pra, f 6 = Цu-sod, f 7 = Nrs, f 8 = Nals.

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

После преобразований производных и из формул (1.1), (1.2) получим следующую систему.



Начальные условия для этой системы: x1=14.2; x2=1.1; x3=1.225; x4=0.992; x5=0; x6=2027; x7=0.921; x8=0.914.

Итак, имеем задачу Коши для системы из 8 дифференциальных уравнений. Пусть , где F - пространство переменных системы. После тщательного анализа исходной системы алгебро-дифференциальных уравнений модели Карааслана было выявлено, что непрерывны и существуют и также непрерывны. Следовательно, правые части системы вместе со своими частными производными непрерывны. Значит, используя теорему о существовании единственности [6], получим, что в решение системы существует и единственно.

Полученная система является автономной, нелинейной. В силу сложности системы, аналитическое исследование задачи на устойчивость провести не удалось. Но при изучении графиков численного решения установлено, что при варьировании входных параметров в небольших пределах, решение также меняется в небольших пределах. На Рис.1.3. представлены графики артериального давления при различных начальных значения параметра x1 - объема внеклеточной жидкости (от 14.2 до 13 л) - и шаге по времени t. Сначала артериальное давление было посчитано при x1= 14.2 л и шаге разбиения по времени t = 1 мин (Рис 1.3.А)). Затем начальное значение x1 было уменьшено до 13 л (Рис 1.3.Б)) И наконец, шаг разбиения по времени был увеличен до 10 мин. (Рис1.3.В)). Как видно из графиков, в любом случае система приходит к равновесному давлению 102.5 мм рт ст. При рассмотрении других величин и варьировании других параметров наблюдается аналогичная картина, что и показывает устойчивость нашей системы.

А)

Б)

В)

Рис. 1.3 - Постоянство среднего артериального давления в системе при различных значениях параметра x1: А) при x1= 14.2 л и шаге t = 1 мин; Б) при x1= 13 л и шаге t = 1 мин; В) x1= 13 л и шаге t = 10 мин.

1.4 Методы численного решения системы уравнений модели

Численное решение данной задачи было реализовано на пакете BioUML

BioUML - Biological Universal Modeling Language - это интегрированный расширяемый пакет программ, созданный с использованием языка программирования Java, который адаптирует подход визуального моделирования для формального описания и симуляции комплексных биологических систем

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

1) структура графа (graph structure) - структура системы описывается компартментализованным графом;

2) уровень базы данных (database level) - каждый элемент графа может содержать ссылку на какой-либо объект базы;

3) математическая модель (mathematical model) - любой элемент графа может быть элементом математической модели. BioUML поддерживает следующие математические элементы: переменные, формулы, уравнения, события, состояния и переходы.

Для работы с моделями BioUML использует две альтернативные моделирующие машины:

1) Java-машину - программа автоматически создает и компилирует Java код на основе визуальной модели (диаграммы) биологической системы.

2) MATLAB-машину - программа автоматически создает код для MATLAB и запускает предустановленный на рабочем месте MATLAB для симуляции модели с использованием библиотеки JMatlink (Patterson and Spiteri, 2003).

Обе моделирующие машины проходят 100% семантических SBML тестов (Finney, 2004). Пусть X=() - вектор неизвестных. В системе BioUML все уравнения разбиты на три блока.

блок скалярных уравнений вида: , i=1…k. Скалярные уравнения просто выражают переменные через переменные X.

блок алгебраических уравнений вида , j=1,…,m. Здесь - набор переменных, которые не входят ни в скалярные ни в дифференциальные уравнения.

блок дифференциальных уравнений вида , l =(k+m)…n.

Схема решения такова: временной интервал, на котором мы хотим получить решение, бьётся на N промежутков. На каждом шаге по времени сначала решаются алгебраические уравнения методом Ньютона. Затем разрешаются скалярные соотношения с использованием значений переменных алгебраических уравнений и значений X полученных на предыдущем шаге по времени. Используя решения алгебраических и скалярных соотношений, можно найти правые части дифференциальных уравнений . Далее приращение l=(k+m),…,n находится методом Рунге-Кутта. В итоге, . В результате получается вектор значений переменных X. Программа переходит к следующему шагу по времени.

Все уравнения модели Карааслана были введены в BioUML (Рис.1.4.). Ниже показано окно, в котором пользователь может задать все необходимые величины для численного моделирования: длину временного интервала, максимальное значение шага разбиения по времени, точность, тип решателя и прочее.

1.5 Численные эксперименты, графики экспериментов

Далее был проведен ряд экспериментов для проверки работоспособности модели. Была смоделирована солевая диета. В эксперименте входной поток соли варьировался следующим образом: сначала задавалось нормальное потребление соли,

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

затем через два дня потребление соли скачком увеличивается примерно в 2 раза, ещё через три с половиной дня потребление соли скачком падает почти до нуля (Рис.1.5.).

Были получены графики, показывающие динамику изменения следующих параметров: баланса соли, то есть разность входного объёма соли и объёма соли, выводящегося с мочой (Рис.1.6.); объёма крови (Рис.1.7.); концентрации антидиуретического гормона (Рис.1.8.); концентрации ренина (Рис.1.9.) и почечной симпатической нервной активности (Рис.1.10.).

Рис. 1.5 - Динамика поступления соли натрия

Рис. 1.6 - Динамика баланса соли натрия



Рис. 1.7 - Динамика изменения объёма крови

Рис. 1.8 - Динамика изменения концентрации антидиуретического гомона (АДГ)



Рис. 1.9 - Динамика изменения концентрации ренина

Рис. 1.10 - Динамика изменения почечной симпатической нервной активности

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

Относительная активность почечного симпатического нерва -

Коэффициент пропорциональности в уравнении , где - сопротивление афферентной артериолы, - влияние почечного симпатического нерва, поток натрия через Мacula densa;

Гидростатическое давление в капсуле Боумена-Шумлянского - ;

Онкотическое давление гломерулярного капилляра ;

Коэффициент фильтрации .

Значение каждого из приведённых параметров увеличивалось и уменьшалось на 50%. Были получены соответствующие графики для каждого из этих параметров. Ниже (Рис.1.11.) приведены графики динамики артериального давления, полученные при варьировании относительной активности почечного симпатического нерва.

Рис. 1.11 - Влияние относительной активности почечного симпатического нерва на артериальное давление

2. МОДЕЛЬ СОЛОДЯННИКОВА

2.1 Основные блоки и биохимические законы модели

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

Система кровообращения может находиться в качественно различных состояниях (фазах). Можно выделить две основные фазы: сокращение желудочков и изгнание крови (систола j=1), пассивное и активное расслабление желудочков и заполнение их кровью (диастола j=2). Вместе с тем наличие системы клапанов означает наличие в системе пороговых нелинейностей. Математической идеализацией такого объекта является динамическая система класса дифференциальных уравнений, когда в различных фазах объект описывается различными системами уравнений при j=1 и j=2.

Рис. 2.1 - Механическая схема модели Солодянникова

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

1) контур управления величиной ударного выброса (закон Франка - Стерлинга)

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

Здесь - конечно-систолический объём желудочка. - конечно-диастолический объём желудочка. -инотропные коэффициенты сердца.

2) контур управления тканевым метаболизмом необходим для поддержания требуемого значения кровотока через капилляры. Управляющий сигнал в нём - кислородный долг .

Здесь - мгновенная доставка кислорода тканям; - содержание кислорода в венозной крови; - кислородный долг; - мгновенная потребность организма в кислороде; - содержание кислорода в артериальной крови; Q - мгновенный поток через капилляр.

3) контур нейрогуморального управления требуется для поддержаниея на необходимом уровне частоты сердечных сокращений. Управляющий сигнал в нём - нейрогуморальный фактор.

Здесь Q - мгновенный поток через капилляры; - объём артерий; - артериальное давление

4) контур управления длительностью систолы отвечает за обеспечение постоянного соотношения между временем заполнения и временем сокращения желудочка.

Также в модели участвуют:

Блок артерий

Блок вен

Блок описания периферической системы (капилляры)

Блок описания процессов, происходящих в желудочке.

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

2.2 Системы дифференциальных уравнений модели Солодянникова

Модель содержит систему из 24-ёх алгебро-дифференциальных уравнений, которая сводится к двум системам из 6 дифференциальных уравнений для каждой фазы. Переменные этих уравнений: - нейрогуморальный фактор; - содержание кислорода в венозной крови; - кислородный долг; - объём желудочка; объём артерий; - систолическое внутрижелудочковое давление. Также вводятся дополнительные переменные: - "запоминает" значение переменной в момент 2-1-перехода (1 - систола, 2 - диастола); - "запоминает" значение в момент 2-1-перехода; - запоминает" значение переменной времени в момент 1-2-перехода. Так же модель содержит список из 30 параметров A1 - A30 с их значениями. Каждый параметр характеризует некоторое индивидуальное для каждого организма свойство и может быть измерен. Даны ещё коэффициенты S1, S2, S3 и инотропные коэффициенты сердца - k1, k2 с их значениями. Далее приведены правые части систем , где i=1…6 - номер уравнения в каждой фазе; j=1,2 - номер фазы (систола, диастола).

X(x, A) = (A16 + A4x3 + A3x1) (ц1(x, A) - ц2(x, A)) - A12A14 (ц1(x, A) - A28) - - A12x1;

X(x, A) X(x, A);

X(x, A) = A1 (A22 - x2) (A16 + A3x1 + A4x3) (ц1(x, A) - ц2(x, A)) - A1A15;

X(x, A) X(x, A);

X(x, A) = - x(x, A);

X(x, A) X(x, A);

X(x, A) = A (ц1(x, A) - x6);

X(x, A) (A18 + A7A15 + A10A28 + A6ц2(x, A)) (ц2(x, A) - ц3(x, A));

X(x, A) = A17(x6 - ц1(x, A)) - (A16 + A3x1 + A4x3) (ц1(x, A) - ц2(x, A));

X(x, A) = - (A16 + A3x1 + A4x3) (ц1(x, A) - ц2(x, A));

X(x, A) X(x, A) 0,

где

PA = ц1(x, A) = A20x5 - (A8A20 + A9A19) x1 + A9x1x5 - A8A9x - A19A20;

PB = ц2(x, A) = (A11A23 - A11A30) x1 - A11x1, x4 - A11x1x5 - A21x4 - A4x5 +A21 (A23 - A30);

РД = ц3(x, A) = A29(A24(x4 - A26)2 + A25(x4 - A26)).

Здесь PA, PB, РД - артериальное, венозное и диастолическое внутрижелудочковое давления соответственно.

Таким образом, мы получили две системы для систолы и диастолы:

Начальные значения переменных:

;;;;;.

Имеем 2 системы дифференциальных уравнений для систолы и диастолы. В начале каждой фазы по значениям с предыдущей фазы для системы ставится задача Коши.

Пусть , где F -пространство переменных . Из приведённых уравнений не сложно показать, что правые части систем (2.1) непрерывны вместе со своими частными производными. Значит, по теореме о существовании единственности (ССЫЛКА НА КНИГУ) решение данных систем существует и единственно в области .

2.3 Переход между фазами модели

Условия перехода от описания системой с j=1 (систола) к описанию системой с j=2 (диастола) имеют вид:

s12(x1, k) = x - (1-k1 )x - k2 = 0 (2.2.)

Это следствие закона Франка - Старлинга для объема желудочка:

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

Условия перехода от описания системой с j=2 (диастола) к описанию системой с j=1 (диастола) имеют вид:

модель карааслан солодянников гемодинамика

- длительность сердечного цикла. Тогда длительность диастолы или t2 определяется из t2-T+T_s=0, что и отражено в уравнении (2.3.)

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



T = T(x2) = +S3 (1-k1)

При переходах происходит следующее:

1) Скачкообразное изменение систолического давления в желудочке:

P_c=P_c0+S1*(T_s-T_s0). Здесь P_c0 - начальное систолическое давление, T_s0 - номинальное систолическое давление - некоторая доля от T: T_s0=S2*T+S3*(1-k1).

2) Скачок скорости изменения объема желудочка (x4) и объема артерий (x5).

Изменение систолического давления в разных циклах (2.4.) осуществляется с тем, чтобы реализовать механизм саморегуляции системы: если время длительности систолы становится большим, становится большой разность T_s-T_s0. Значит, в следующем сердечном цикле увеличится P_c и уменьшится длительность систолы.

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

2.4 Реализация численного решения

Модель Солодянникова была реализована с помощью программы, написанной на языке программирования MATLAB7 (Приложение 2). Данная программа предназначена для проведения тестовых расчетов по модели Солодянникова, моделирования различных патологий и ситуаций в кровеносной системе человека. Моделирование патологий осуществляется с помощью варьирования значений параметров A1 - A30, которые задаются в программе. Данная программа использует метод Рунге - Кутта [7] для решения систем дифференциальных уравнений. Временной интервал, на котором мы хотим получить решение, бьется на N промежутков (цикл по времени - цикл по переменной k). На каждом промежутке используя правые части систем (2.1.)методом Рунге _ Кутта находятся приращения переменных системы на текущем временном промежутке. Это делается следующим образом. Пусть - значения переменных на предыдущем шаге по времени, при фиксированном j. Тогда приращения этих переменных на текущем временном шаге находятся по формулам

;

;

,

На каждом шаге по времени происходит проверка условий (2.2.) и (2.3.). Если одно из этих условий выполняется, то происходит соответствующий переход к расчетам системы противоположной фаза (от расчетов систолы к расчетам диастолы или наоборот). При этом вектор значений переменных модели X0 (см. приложение 2) передаётся в качестве начальных значений для решения системы противоположной фазы.

Входными данными для программы являются длинна временного интервала T, на котором мы хотим найти решение и шаг разбиения по времени tau. Так же внутри программы задаются значения всех параметров A1 - A30 и коэффициентов S1, S2, S3, k1, k2, а так же начальные значения переменных .

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

Правильность работы программы можно проверить командой `Runge2(10,0.05)' при этом должен появиться следующий график артериального давления (Рис. 2.2).

Рис. 2.2 - Нормальное значение артериального давления

2.5 Моделирование патологий, графики экспериментов

С помощью модели Солодянникова были смоделированы кардиогенный шок и физическая нагрузка.

Кардиогенный шок - это острая левожелудочковая недостаточность крайней степени тяжести, развивающаяся при инфаркте миокарда. Уменьшение ударного и минутного объема крови при шоке настолько выражено, что не компенсируется повышением сосудистого сопротивления, вследствие чего резко снижаются артериальное давление и системный кровоток, нарушается кровоснабжение всех жизненно важных органов. Моделирование этой ситуации осуществлялось уменьшением инотропного коэффициента k1 от значения 0,6 до 0,3 на шестидесятой секунде эксперимента.

Ниже на графиках представлена динамика изменения артериального давления (Рис.2.3.), объёма левого желудочка (Рис.2.4.) и кислородного долга (Рис.2.5.)

Как видно на графиках на шестидесятой секунде эксперимента наблюдается скачкообразное падение давления и рост объёма желудочка. Далее среднее значение давления падает почти линейно и пульс (частота колебаний давления) уменьшается

Рис. 2.3 - Динамика артериального давления при кардиогенном шоке



Рис. 2.4 - Динамика объёма левого желудочка при кардиогенном шоке

Рис. 2.5 - Динамика кислородного долга при кардиогенном шоке

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

Под увеличением физической нагрузки понимается механическая работа мускулатуры, ведущая к усилению процессов обмена, что приводит к увеличению потребления организмом кислорода. Для увеличения объёмов поставки кислорода к органам и тканям происходит увеличение частоты сердечных сокращений, повышение артериального давления, увеличение ударного и минутного объема крови. Физическая нагрузка моделировалась увеличением потребности тканей в кислороде: На двадцатой секунде эксперимента значение параметра A15 было увеличено в 7 раз.

На графиках изображены артериальное давление (Рис. 2.6), тканевый кровоток (Рис. 2.7), частота сердечных сокращений (Рис. 2.8)

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

Рис. 2.7 - Динамика тканевого кровотока при физической нагрузке



Рис. 2.8 - Динамика частоты сердечных сокращений при физической нагрузке

Начиная с двадцатой секунды амплитуда давления возрастает со значения 120 на 80 мм.рт.ст. до значения 140 на 90 мм.рт.ст. на пятидесятой секунде. Частота сердечных сокращений то же возрастает. Растёт и среднее значение тканевого кровотока, чтобы удовлетворить возросшие потребности организма в кислороде.

3. КОМПЛЕКСНАЯ МОДЕЛЬ ГЕМОДИНАМИКИ

3.1 Построение комплексной модели гемодинамики

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

Эта модель включает в себя 55 основных артерий тела человека (Рис. 3.1.), характеризующихся собственными параметрами, такими как длина, поперечное сечение, удаленность от сердца и эластичность стенок.

Рис. 3.1 - Схема 55 артерий человека - "артериальное древо" - из гидродинамической модели.

Данная модель позволяет рассчитать поток крови и площадь сечения в каждой точке каждого из 55 сосудов. А так же давление крови в любой точке этих сосудов. Модель была построена следующим образом. Трёхмерные уравнения Навье-Стокса были приведены к одномерным. Была получена система из 2-ух уравнений в частных производных для потока крови и площади сечения сосуда. В уравнения вошли производные по времени и по длине сосуда. Затем производная по времени была заменена своим разностным аналогом. Далее был разработан метод ортогональной прогонки для решения этой системы. Написана программа на системе MATLAB7.

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

При создании гидродинамической модели остро стоял вопрос о задании зависимости давления от времени в первом сосуде на выходе из сердца. Данная проблема была успешно решена с помощью моделей Карааслана и Солодянникова. В результате была реализована комплексная модель гемодинамики. К управляющим контурам модели Солодянникова была добавлена почечная регуляция модели Карааслана. Модель артериального резервуара была дополнена гидродинамической моделью из 55 артерий (Рис. 3.2).

Рис. 3.2 - Схема комплексной модели гемодинамики, созданной на основе моделей Карааслана, Солодянникова и гидродинамической модели

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

3.2 Эксперименты

Чтобы более наглядно проиллюстрировать возможности данной модели, рассмотрим пример. Представим себе человека, идущего по пустыне под палящим солнцем. Испарение забирает из тела в два раза больше влаги, чем выводится обычно. У путника ограничен запас воды, поэтому он пьёт 30% от обычного потребления воды. После трёх дней пути странник приходит к оазису, где напивается воды (пьёт в четыре раза больше, чем обычно) и отдыхает в тени два дня (потоотделение становится обычным), после чего снова отправляется в путь. На графике (Рис.3.3.) изображена динамика среднего артериального давления у нашего путника за весь период путешествия. Хорошо видно, что после трёх дней пути (~4300 мин) среднее давление сильно падает до 70 мм рт ст. Затем за 2 дня вырастает выше нормы - до 120 мм рт ст, а затем снова падает.

Рис. 3.3 - Динамика среднего артериального давления путешественника. Цифрами обозначено: 1) отправная точка, 2) приход в оазис, 3) уход из оазиса, 4) дополнительная точка для исследования при тахикардии

Допустим, путник страдал тахикардией (патологический симптом, характеризующийся увеличением частоты сердечных сокращений в покое от 90 ударов в минуту). Для моделирования этой патологии длина каждого сердечного цикла в модели Солодянникова задавалась равномерно распределённым на отрезке (0.8 сек; 1.6 сек) случайным числом. График, отражающий динамику артериального давления в этом случае, приведен на Рис. 3.4.

Проследим динамику изменения артериального давления в середине первого (восходящая ветвь аорты) и тринадцатого (правая внешняя сонная артерия) сосудов в течение путешествия, т.е. в точках 1-4 (см. Рис.3.3.). Для этого, в соответствии со значениями среднего давления в этих точках, понизим давление, полученное при моделировании тахикардии, на выходе из сердца в точке 2 (см. Рис.3.3.) на 30% и в точке 4 на 12%; а в точке 3 повысим давление на выходе из сердца на 20%. В точке 1 давление остаётся неизменным. Полученные зависимости давления от времени подадим на вход в гидродинамическую модель. Графики давлений в интересующих нас сосудах во время эксперимента представлены на Рис. 3.4.

Рис. 3.4 - График артериального давления страдающего тахикардией человека

А)

Б)

В)

Г)

Рис. 3.5 - Динамика давления в середине правой внешней сонной артерии (более тёмная линия) в сравнении с давлением в середине восходящей ветви аорты (более светлая линия) у страдающего тахикардией: А) в точке путешествия 1 (см. Рис.3.3.); Б) в точке путешествия 2; В) в точке путешествия 3; Г) в точке путешествия 4

ЗАКЛЮЧЕНИЕ

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

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

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

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

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

ЛИТЕРАТУРА

1. Karaaslan F, Denizhan Y, Kayserilioglu A, Gulcur HO. (2005) Long-term mathematical model involving renal sympathetic nerve activity, arterial pressure, and sodium excretion. Ann Biomed Eng. V.33. P.1607-30.

2. Солодянников Ю.В. Элементы математического моделирования и идентификация системы кровообращения. Самара: Самар. ун-т, 1994. 315 с.

3. Guyton, A.C., T.G. Coleman, and H. J. Granger. (1972) Circulation: Overall regulation. Ann. Rev. Physiol. V.34. P.13-46.

4. Coleman T.G., and J.E. Hall. A mathematical model of renal hemodynamics and excretory function. In: Structuring Biological Systems: A Computer Modelling Approach, edited by S. S. Iyengar. Boca Raton, Florida: CRC Press, 1992, pp. 89-124.

5. Uttamsingh R.J., M.S. Leaning, J.A. Bushman, E.R. Carson, and L. Finkelstein. (1985) Mathematical model of the human renal system. Med. Biol. Eng. Comput. V.23. P.525-535.

6. Понтрягин Л.С. Обыкновенные дифференциальные уравнения. Ижевск: НИЦ “Регулярная и хаотическая динамика”, 2001. 400 с.

7. Березин И.С., Жидков Н.П. Методы вычислений. 2 изд., перераб. М.: Гос. Изд-во физ.-мат. лит., 1962. Т.2. С.292-326.

8. Lamponi D.N One dimentional and multiscale models for blood flow circulation. Lausanne, EPFL, 2004.

ПРИЛОЖЕНИЕ 1

Таблица П.1 - Условные обозначения, принятые в модели Карааслана (значения параметров без размерности даны в абсолютных величинах)

Символ	Расшифровка	Нормальное равновесное значение
бmap	Эффект среднего артериального давления на почечную симпатическую нервную активность	1
бrap	Эффект давления в правом предсердии на почечную симпатическую нервную активность	1
вsna	Эффект почечной симпатической нервной активности на сопротивление афферентной артериолы	1
дra	Эффект давления в правом предсердии на нормализованную скорость секреции антидиуретического гормона	0
гat	Эффект концентрации ангиотензина на реабсорбцию натрия проксимальным канальцем	1
гfilsod	Эффект нагрузки профильтровавшегося натрия на его реабсорбцию проксимальным канальцем	1
гrsna	Эффект почечной симпатической нервной активности на реабсорбцию натрия проксимальным канальцем	1
еaum	Эффект автономного множителя	1
зcd-sodreab	Реабсорбция натрия собирательной трубкой (доля от проходящего количества)	0.93
зdt-sodreab	Реабсорбция натрия дистальным канальцем (доля от проходящего количества)	0.5
зpt-sodreab	Реабсорбция натрия проксимальным канальцем (доля от проходящего количества)	0.8
лanp	Эффект натриуретического пептида на скорость реабсорбции натрия собирательными трубками	1
лdt	Эффект тока натрия, прошедшего дистальный каналец, на скорость его реабсорбции в собирательных трубках	1
мadh	Эффект концентрации антидиуретического гормона на скорость реабсорбции воды канальцами	1
мal	Эффект концентрации альдостерона на скорость реабсорбции воды канальцами	1
нmd-sod	Эффект тока натрия через Macula densa на скорость секреции ренина	1
нrsna	Эффект почечной симпатической нервной активности на скорость секреции ренина	1
оat	Эффект ангиотензина на скорость секреции альдостерона	1
оk/sod	Эффект отношения концентраций калия/натрия на скорость секреции альдостерона	1
оmap	Эффект среднего артериального давления на скорость секреции альдостерона	1
Уtgf	Сигнал тубулогломерулярной обратной связи	1
Цcd-sodreab	Абсолютная скорость реабсорбции натрия собирательной трубкой	1.674 мЭк/мин
Цco	Сердечный выброс (минутный объем)	5 л/мин
Цdt-sod	Поток натрия, выходящий из дистального канальца	1.8 мЭк/мин
Цdt-sodreab	Абсолютная скорость реабсорбции натрия дистальным канальцем	1.8 мЭк/мин
Цfilsod	Нагрузка профильтровавшегося натрия	18 мЭк/мин
Цgfilt	Скорость клубочковой фильтрации	0.125 л/мин
Цmd-sod	Скорость тока натрия через Macula densa	3.6 мЭк/мин
Цpt-sodreab	Абсолютная скорость реабсорбции натрия проксимальным канальцем	14.4 мЭк/мин
Фrb	Скорость тока крови в почке	1.2 л/мин
Фsodin	Потребление соли	0.126 мЭк/мин
Фt-wreab	Скорость реабсорбции воды в канальцах	0.124 л/мин
Фu	Скорость тока мочи	0.001 л/мин
Фu-sod	Поток натрия в составе мочи	0.126 мЭк/мин
Фvr	Венозный возврат	5 л/мин
Фwin	Потребление воды	0.001 л/мин
шal	Эффект альдостерона на реабсорбцию натрия дистальным канальцем	1
аauto	Активность автономной системы	1
аbaro	Активность барорецепторов	0.75
achemo	Активность хеморецепторов	0.25
Cadh	Концентрация антидиуретического гормона	4*10-3 частиц/л
Cal	Концентрация альдостерона	85 нг/л
Canp	Нормализованная концентрация натриуретического пептида	1
Cat	Концентрация ангиотензина	20 нг/л
Cgcf	Коэффициент фильтрации клубочковых капилляров	0.00781
Ck	Концентрация калия в крови	5 мЭк/л
Cr	Нормализованная концентрация ренина	1
Csod	Концентрация натрия в крови	144 мЭк/л
Kbar	Коэффициент, связывающий базовое сопротивление артерий с васкуляризацией (обеспеченностью сосудами)	16.6 мм рт ст* мин/л
Kvd	Коэффициент снижения числа сосудов	0.00001
Msod	Общее количество натрия	2160 мЭк
Nadh	Нормализованная концентрация антидиуретического гормона	1
Nals	Нормализованная скорость секреции альдостерона	1
nе-dt	Нормальное значение реабсорбции натрия дистальным канальцем	0.5
nз-cd	Нормальное значение реабсорбции натрия собирательной трубкой (доля от проходящего количества)	0.93
nз-pt	Нормальное значение реабсорбции натрия проксимальным канальцем (доля от проходящего количества)	0.8
Nal	Нормализованная концентрация альдостерона	1
Nrs	Нормализованная скорость секреции ренина	1
Nrsna	Нормализованная почечная симпатическая активность	1
PB	Гидростатическое давление в капсуле Боумена-Шумлянского	18 мм рт ст
Pt	Фильтрационное давление в сети капилляров	16 мм рт ст
Pgh	Клубочковое гидростатическое давление	62 мм рт ст
Pgo	Клубочковое осмотическое давление	28 мм рт ст
Pma	Среднее артериальное давление	100 мм рт ст
Pmf	Среднее наполняющее давление	7 мм рт ст
Pra	Давление в правом предсердии	0 мм рт ст
Ra	Артериальное сопротивление	16.6 мм рт ст* мин/л
Raa	Сопротивление афферентной артериолы	31.67 мм рт ст* мин/л
Rba	Базовое артериальное сопротивление	16.6 мм рт ст* мин/л
Rbv	Базовое венозное сопротивление	3.4 мм рт ст* мин/л
Rea	Сопротивление эфферентной артериолы	51.66 мм рт ст* мин/л
rsna	Почечная симпатическая нервная активность	1
Rtp	Общее периферическое сопротивление	20 мм рт ст* мин/л
Rr	Почечное сосудистое сопротивление	83.33 мм рт ст* мин/л
Rvr	Сопротивление венозному возврату	1.4 мм рт ст* мин/л
Nadhs	Нормализованная скорость секреции антидиуретического гормона	1
Tadh	Константа времени для секреции антидиуретического гормона	6 мин
Tal	Константа времени для секреции альдостерона	30 мин
Tr	Константа времени для секреции ренина	15 мин
vas	Васкуляризация	1
vasd	Скорость снижения васкуляризации	0.00001
vasf	Скорость увеличения васкуляризации	0.00001
Vb	Объем крови	5 л
Vecf	Объем внеклеточной жидкости	15 л

ПРИЛОЖЕНИЕ 2

function [M]=Runge2(T, tau)

function Fi1=Fi1(x1,x5,A20,A8,A9,A19)

Fi1=A20*x5-(A8*A20+A9*A19)*x1+A9*x1*x5-A8*A9*x1*x1-A19*A20;

end;

function Fi2=Fi2(x1,x4,x5,A11,A21,A23,A30)

Fi2=(A11*A23-A11*A30)*x1-A11*x1*x4-A11*x1*x5-A21*x4-A21*x5+A21*(A23-A30);

end;

function Fi3=Fi3(x4,A24,A25,A26,A29)

Fi3=A29*(x4-A26)*(A24*(x4-A26)+A25);

end;

function X_1=X_1(x1,x5,A12,A13,A16,A4,A3,A14,A28,A20,A8,A9,A19,x4,A11,A21,A23,A30)

X_1=((A12*A13)*(A16+A4*x3+A3*x1)*(Fi1(x1,x5,A20,A8,A9,A19)-Fi2(x1,x4,x5,A11,A21,A23,A30)))/x5-A12*A14*(Fi1(x1,x5,A20,A8,A9,A19)-A28)-A12*x1;

end;

function X_2=X_2(x1,x2,x3,A1,A22,A16,A15,A4,A3,x5,A20,A8,A9,A19,x4,A11,A21,A23,A30)

X_2=(A1*(A22-x2)*(A16+A3*x1+A4*x3)*(Fi1(x1,x5,A20,A8,A9,A19)-Fi2(x1,x4,x5,A11,A21,A23,A30))-A1*A15);

end;

function X_3=X_3(x1,x2,x3,A1,A22,A16,A15,A4,A3,A2,x5,A20,A8,A9,A19)

X_3=(-A2/A1)*X_2(x1,x2,x3,A1,A22,A16,A15,A4,A3,x5,A20,A8,A9,A19,x4,A11,A21,A23,A30);

end;

function X_4=X_4(x6,A17,x1,x5,A20,A8,A9,A19)

X_4=A17*(Fi1(x1,x5,A20,A8,A9,A19)-x6);

end;

function X_4_2=X_4_2(A18,A17,A15,A7,A10,A28,A6,x1,x4,x5,A11,A21,A23,A30,A24,A25,A26,A29)

X_4_2=(A18+A7*A15-A10*A28+A6*Fi2(x1,x4,x5,A11,A21,A23,A30))*(Fi2(x1,x4,x5,A11,A21,A23,A30)-Fi3(x4,A24,A25,A26,A29));

end;

function X_5=X_5(x1,x3,x6,A17,A16,A4,A3,x5,A20,A8,A9,A19,x4,A11,A21,A23,A30)

X_5=A17*(x6-Fi1(x1,x5,A20,A8,A9,A19))-(A16+A3*x1+A4*x3)*(Fi1(x1,x5,A20,A8,A9,A19)-Fi2(x1,x4,x5,A11,A21,A23,A30));

end;

function X_5_2=X_5_2(x1,x3,A16,A4,A3,x5,A20,A8,A9,A19,x4,A11,A21,A23,A30)

X_5_2=-(A16+A3*x1+A4*x3)*(Fi1(x1,x5,A20,A8,A9,A19)-Fi2(x1,x4,x5,A11,A21,A23,A30));

end;

N=round(T/tau);

A1=0.0008;A6=3.3;A11=0.00000185;

A2=0.5;A7=4.0;A12=2.0;

A3=0.626;A8=120;A13=6.0;

A4=0.00665;A9=0.3;A14=0.0015;

A5=22;A10=0.8;A15=4;

A16=0.0733;A25=0.023;

A17=11.7;A26=120;

A18=120;A27=125;P0=125;

A19=80;A28=70;

A20=0.35;A29=0.007;

A21=0.0006694; A30=0.2;

A22=0.17;A23=3500;A24=0.01;

k1=0.6; k2=20;

S1=22;S2=0.25;S3=0.2;

x1=1.1;

x2=0.14;

x3=0;

x4=156;

x5=360;

x6=125;

x7=x1;

x8=x4;

x9=0;

X0=[x1 x2 x3 x4 x5 x6 x7 x8 x9];

t=0;t1=0;

c=0;

if (N>=2000) nn=round(N/2000);Mas=zeros(round(N/nn),9);time=zeros(round(N/nn));

else nn=1; Mas=zeros(N,9);end;

for k = 1:N

t=t+tau;

t1=t1+tau;

m=zeros(4,5); m(1,1)=tau*X_1(X0(1),X0(5),A12,A13,A16,A4,A3,A14,A28,A20,A8,A9,A19,X0(4),A11,A21,A23,A30);

m(1,2)=tau*X_2(X0(1),X0(2),X0(3),A1,A22,A16,A15,A4,A3,X0(5),A20,A8,A9,A19,X0(4),A11,A21,A23,A30);

m(1,3)=tau*X_3(X0(1),X0(2),X0(3),A1,A22,A16,A15,A4,A3,A2,X0(5),A20,A8,A9,A19);

if c==0

m(1,4)=tau*X_4(X0(6),A17,X0(1),X0(5),A20,A8,A9,A19);

m(1,5)=tau*X_5(X0(1),X0(3),X0(6),A17,A16,A4,A3,X0(5),A20,A8,A9,A19,X0(4),A11,A21,

A23,A30); end;

if (c==1)

m(1,4)=tau*X_4_2(A18,A17,A15,A7,A10,A28,A6,X0(1),X0(4),X0(5),A11,A21,A23,A30,A24,A25,A26,A29);

m(1,5)=tau*X_5_2(X0(1),X0(3),A16,A4,A3,X0(5),A20,A8,A9,A19,X0(4),A11,A21,A23,A30); end;

for j = 2:4

m(j,1)=tau*X_1((X0(1)+m(j-1,1)/2)*((j==2)+(j==3))+(X0(1)+m(j1,1))*(j==4),

(X0(5)+m(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+

m(j-1,5))*(j==4),A12,A13,A16,A4,A3,A14,A28,A20,A8,A9,A19,(X0(4)+

m(j-1,4)/2)*((j==2)+(j==3))+(X0(4)+m(j-1,4))*(j==4),A11,A21,A23,A30);

m(j,2)=tau*X_2((X0(1)+m(j-1,1)/2)*((j==2)+(j==3))+(X0(1)+

m(j-1,1))*(j==4),(X0(2)+m(j-1,2)/2)*((j==2)+(j==3))+(X0(2)+m(j-1,2))*(j==4),(X0(3)+

m(j-1,3)/2)*((j==2)+(j==3))+(X0(3)+m(j-1,3))*(j==4),A1,A22,A16,A15,A4,A3,(X0(5)+

m(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+m(j-1,5))*(j==4),A20,A8,A9,A19,(X0(4)+

m(j-1,4)/2)*((j==2)+(j==3))+(X0(4)+m(j-1,4))*(j==4),A11,A21,A23,A30);

m(j,3)=tau*X_3((X0(1)+m(j-1,1)/2)*((j==2)+(j==3))+(X0(1)+

m(j-1,1))*(j==4),(X0(2)+m(j-1,2)/2)*((j==2)+(j==3))+(X0(2)+m(j-1,2))*(j==4),(X0(3)+

m(j-1,3)/2)*((j==2)+(j==3))+(X0(3)+m(j-1,3))*(j==4),A1,A22,A16,A15,A4,A3,A2,(X0(5)+

m(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+m(j-1,5))*(j==4),A20,A8,A9,A19);

if c==0

m(j,4)=tau*X_4(X0(6),A17,(X0(1)+m(j-1,1)/2)*((j==2)+(j==3))+(X0(1)+

m(j-1,1))*(j==4),(X0(5)+m(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+

m(j-1,5))*(j==4),A20,A8,A9,A19);

m(j,5)=tau*X_5((X0(1)+m(j-1,1)/2)*((j==2)+(j==3))+(X0(1)+

m(j-1,1))*(j==4),(X0(3)+m(j-1,3)/2)*((j==2)+(j==3))+(X0(3)+

m(j-1,3))*(j==4),X0(6),A17,A16,A4,A3,(X0(5)+m(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+

m(j-1,5))*(j==4),A20,A8,A9,A19,(X0(4)+m(j-1,4)/2)*((j==2)+(j==3))+(X0(4)+

m(j-1,4))*(j==4),A11,A21,A23,A30);

end;

if c==1

m(j,4)=tau*X_4_2(A18,A17,A15,A7,A10,A28,A6,(X0(1)+

m(j-1,1)/2)*((j==2)+(j==3))+(X0(1)+m(j-1,1))*(j==4),(X0(4)+

m(j-1,4)/2)*((j==2)+(j==3))+(X0(4)+m(j-1,4))*(j==4),(X0(5)+

m(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+m(j-1,5))*(j==4),A11,A21,A23,A30,A24,A25,A26,A29);

m(j,5)=tau*X_5_2((X0(1)+m(j-1,1)/2)*((j==2)+(j==3))+(X0(1)+

m(j-1,1))*(j==4),(X0(3)+m(j-1,3)/2)*((j==2)+(j==3))+(X0(3)+

m(j-1,3))*(j==4),A16,A4,A3,(X0(5)+m(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+

m(j-1,5))*(j==4),A20,A8,A9,A19,(X0(4)+m(j-1,4)/2)*((j==2)+(j==3))+(X0(4)+

m(j-1,4))*(j==4),A11,A21,A23,A30);

for j=1:5

R(j)=X0(j)+(m(1,j)+2*m(2,j)+2*m(3,j)+m(4,j))/6;

X0(j)=R(j);

end;

if (mod(k,nn)==0) Mas(round(k/nn),:)=X0; time(round(k/nn))=t;

Mas(round(k/nn),5)=Fi1(X0(1),X0(5),A20,A8,A9,A19);

end;

if (c==0) && (k>1)

S12(k)=X0(4)-(1-k1)*X0(8)-k2;