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

ВВЕДЕНИЕ

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

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

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

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

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

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

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

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

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

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

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

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

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), = ?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 мин.

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

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

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

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

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

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

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

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

) 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.). Ниже показано окно, в котором пользователь может задать все необходимые величины для численного моделирования: длину временного интервала, максимальное значение шага разбиения по времени, точность, тип решателя и прочее.

.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. МОДЕЛЬ СОЛОДЯННИКОВА

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Блок артерий

Блок вен

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

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

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

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, A) X(x, A);(x, A) = A1 (A22 - x2) (A16 + A3x1 + A4x3) (?1(x, A) - ?2(x, A)) - A1A15;(x, A) X(x, A);(x, A) = - x(x, A);(x, A) X(x, A);(x, A) = A (?1(x, A) - x6);(x, A) (A18 + A7A15 + A10A28 + A6?2(x, A)) (?2(x, A) - ?3(x, A));(x, A) = A17(x6 - ?1(x, A)) - (A16 + A3x1 + A4x3) (?1(x, A) - ?2(x, A));(x, A) = - (A16 + A3x1 + A4x3) (?1(x, A) - ?2(x, A));(x, A) X(x, A) 0,

где= ?1(x, A) = A20x5 - (A8A20 + A9A19) x1 + A9x1x5 - A8A9x - A19A20;= ?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(x2) = +S3 (1-k1)

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

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

) Скачок скорости изменения объема желудочка (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. КОМПЛЕКСНАЯ МОДЕЛЬ ГЕМОДИНАМИКИ

.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 с.

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

.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.

.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.

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

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

.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 densa3.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.75achemoАктивность хеморецепторов0.25CadhКонцентрация антидиуретического гормона4*10-3 частиц/лCalКонцентрация альдостерона85 нг/лCanpНормализованная концентрация натриуретического пептида1CatКонцентрация ангиотензина20 нг/лCgcfКоэффициент фильтрации клубочковых капилляров0.00781CkКонцентрация калия в крови5 мЭк/лCrНормализованная концентрация ренина1CsodКонцентрация натрия в крови144 мЭк/лKbarКоэффициент, связывающий базовое сопротивление артерий с васкуляризацией (обеспеченностью сосудами)16.6 мм рт ст* мин/лKvdКоэффициент снижения числа сосудов0.00001MsodОбщее количество натрия2160 мЭкNadhНормализованная концентрация антидиуретического гормона1NalsНормализованная скорость секреции альдостерона1n?-dtНормальное значение реабсорбции натрия дистальным канальцем0.5n?-cdНормальное значение реабсорбции натрия собирательной трубкой (доля от проходящего количества)0.93n?-ptНормальное значение реабсорбции натрия проксимальным канальцем (доля от проходящего количества)0.8NalНормализованная концентрация альдостерона1NrsНормализованная скорость секреции ренина1NrsnaНормализованная почечная симпатическая активность1PBГидростатическое давление в капсуле Боумена-Шумлянского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Почечная симпатическая нервная активность1RtpОбщее периферическое сопротивление20 мм рт ст* мин/лRrПочечное сосудистое сопротивление83.33 мм рт ст* мин/лRvrСопротивление венозному возврату1.4 мм рт ст* мин/лNadhsНормализованная скорость секреции антидиуретического гормона1TadhКонстанта времени для секреции антидиуретического гормона6 минTalКонстанта времени для секреции альдостерона30 минTrКонстанта времени для секреции ренина15 минvasВаскуляризация1vasdСкорость снижения васкуляризации0.00001vasfСкорость увеличения васкуляризации0.00001VbОбъем крови5 лVecfОбъем внеклеточной жидкости15 лПРИЛОЖЕНИЕ 2

function [M]=Runge2(T, tau)Fi1=Fi1(x1,x5,A20,A8,A9,A19)=A20*x5-(A8*A20+A9*A19)*x1+A9*x1*x5-A8*A9*x1*x1-A19*A20;;Fi2=Fi2(x1,x4,x5,A11,A21,A23,A30)=(A11*A23-A11*A30)*x1-A11*x1*x4-A11*x1*x5-A21*x4-A21*x5+A21*(A23-A30);;Fi3=Fi3(x4,A24,A25,A26,A29)=A29*(x4-A26)*(A24*(x4-A26)+A25);;X_1=X_1(x1,x5,A12,A13,A16,A4,A3,A14,A28,A20,A8,A9,A19,x4,A11,A21,A23,A30)_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;;X_2=X_2(x1,x2,x3,A1,A22,A16,A15,A4,A3,x5,A20,A8,A9,A19,x4,A11,A21,A23,A30)_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);;X_3=X_3(x1,x2,x3,A1,A22,A16,A15,A4,A3,A2,x5,A20,A8,A9,A19)_3=(-A2/A1)*X_2(x1,x2,x3,A1,A22,A16,A15,A4,A3,x5,A20,A8,A9,A19,x4,A11,A21,A23,A30);;X_4=X_4(x6,A17,x1,x5,A20,A8,A9,A19)_4=A17*(Fi1(x1,x5,A20,A8,A9,A19)-x6);;X_4_2=X_4_2(A18,A17,A15,A7,A10,A28,A6,x1,x4,x5,A11,A21,A23,A30,A24,A25,A26,A29)_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));;X_5=X_5(x1,x3,x6,A17,A16,A4,A3,x5,A20,A8,A9,A19,x4,A11,A21,A23,A30)_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));;X_5_2=X_5_2(x1,x3,A16,A4,A3,x5,A20,A8,A9,A19,x4,A11,A21,A23,A30)_5_2=-(A16+A3*x1+A4*x3)*(Fi1(x1,x5,A20,A8,A9,A19)-Fi2(x1,x4,x5,A11,A21,A23,A30));;=round(T/tau);=0.0008;A6=3.3;A11=0.00000185;=0.5;A7=4.0;A12=2.0;=0.626;A8=120;A13=6.0;=0.00665;A9=0.3;A14=0.0015;=22;A10=0.8;A15=4;=0.0733;A25=0.023;=11.7;A26=120;=120;A27=125;P0=125;=80;A28=70;=0.35;A29=0.007;=0.0006694; A30=0.2;=0.17;A23=3500;A24=0.01;=0.6; k2=20;=22;S2=0.25;S3=0.2;=1.1;=0.14;=0;=156;=360;=125;=x1;=x4;=0;=[x1 x2 x3 x4 x5 x6 x7 x8 x9];=0;t1=0;=0;(N>=2000) nn=round(N/2000);Mas=zeros(round(N/nn),9);time=zeros(round(N/nn));nn=1; Mas=zeros(N,9);end;k = 1:N

t=t+tau;

t1=t1+tau;=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);(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);(1,3)=tau*X_3(X0(1),X0(2),X0(3),A1,A22,A16,A15,A4,A3,A2,X0(5),A20,A8,A9,A19);c==0 (1,4)=tau*X_4(X0(6),A17,X0(1),X0(5),A20,A8,A9,A19);(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,,A30); end;(c==1) (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);(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;j = 2:4(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)+(j-1,5))*(j==4),A12,A13,A16,A4,A3,A14,A28,A20,A8,A9,A19,(X0(4)+(j-1,4)/2)*((j==2)+(j==3))+(X0(4)+m(j-1,4))*(j==4),A11,A21,A23,A30);(j,2)=tau*X_2((X0(1)+m(j-1,1)/2)*((j==2)+(j==3))+(X0(1)+(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)+(j-1,3)/2)*((j==2)+(j==3))+(X0(3)+m(j-1,3))*(j==4),A1,A22,A16,A15,A4,A3,(X0(5)+(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+m(j-1,5))*(j==4),A20,A8,A9,A19,(X0(4)+(j-1,4)/2)*((j==2)+(j==3))+(X0(4)+m(j-1,4))*(j==4),A11,A21,A23,A30);(j,3)=tau*X_3((X0(1)+m(j-1,1)/2)*((j==2)+(j==3))+(X0(1)+(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)+(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)+(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+m(j-1,5))*(j==4),A20,A8,A9,A19);c==0 (j,4)=tau*X_4(X0(6),A17,(X0(1)+m(j-1,1)/2)*((j==2)+(j==3))+(X0(1)+(j-1,1))*(j==4),(X0(5)+m(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+(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)+(j-1,1))*(j==4),(X0(3)+m(j-1,3)/2)*((j==2)+(j==3))+(X0(3)+(j-1,3))*(j==4),X0(6),A17,A16,A4,A3,(X0(5)+m(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+(j-1,5))*(j==4),A20,A8,A9,A19,(X0(4)+m(j-1,4)/2)*((j==2)+(j==3))+(X0(4)+(j-1,4))*(j==4),A11,A21,A23,A30);;c==1 (j,4)=tau*X_4_2(A18,A17,A15,A7,A10,A28,A6,(X0(1)+(j-1,1)/2)*((j==2)+(j==3))+(X0(1)+m(j-1,1))*(j==4),(X0(4)+(j-1,4)/2)*((j==2)+(j==3))+(X0(4)+m(j-1,4))*(j==4),(X0(5)+(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+m(j-1,5))*(j==4),A11,A21,A23,A30,A24,A25,A26,A29);(j,5)=tau*X_5_2((X0(1)+m(j-1,1)/2)*((j==2)+(j==3))+(X0(1)+(j-1,1))*(j==4),(X0(3)+m(j-1,3)/2)*((j==2)+(j==3))+(X0(3)+(j-1,3))*(j==4),A16,A4,A3,(X0(5)+m(j-1,5)/2)*((j==2)+(j==3))+(X0(5)+(j-1,5))*(j==4),A20,A8,A9,A19,(X0(4)+m(j-1,4)/2)*((j==2)+(j==3))+(X0(4)+(j-1,4))*(j==4),A11,A21,A23,A30);j=1:5(j)=X0(j)+(m(1,j)+2*m(2,j)+2*m(3,j)+m(4,j))/6;(j)=R(j);;(mod(k,nn)==0) Mas(round(k/nn),:)=X0; time(round(k/nn))=t;(round(k/nn),5)=Fi1(X0(1),X0(5),A20,A8,A9,A19);

end;(c==0) && (k>1)(k)=X0(4)-(1-k1)*X0(8)-k2;(S12(k)*S12(k-1)<0) && (c==0)(9)=t1;=0;=1; end; end;c==1=t1-1/(X0(7))+X0(9);(S21>0) &&(c==1)s=S2/(X0(7))+S3*(1-k1);(6)=X0(6)+S1*(X0(9)-T0s);(7)=X0(1);(8)=X0(4);(9)=0;=0;=0;end;end;; F=figure; figure(F); (time,Mas(:,5),'blue'); hold; pause;('C:\t.txt',Mas1(:,5),'precision', '%.6f','newline','pc');

end