Математическое моделирование в медицине и биологии на основе моделей механики сплошных сред

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

Разработка этого аппарата предполагает:

• —    построение замкнутой механикоматематической модели процесса, описывающей поведение биологической среды на основе системы уравнений в частных производных механики сплошных сред (МСС);

• —    разработку замыкающих систему МСС реологических соотношений, описывающих поведение той или иной среды (для гидродинамики это уравнения состояния, для механики деформируемого твёрдого тела (МДТТ) — соотношения между компонентами тензоров напряжений и деформаций);

• —    корректную математическую постановку задачи, то есть представление замкнутой системы МСС, постановку необходимых для её решения начальных и граничных условий, условий на контактных границах (если они есть);

• —    разработку или реализацию вычислительных методов, адаптированных к специфике решения конкретной задачи;

• —    разработку алгоритма численного решения задачи и его программную реализацию;

• —    численное решение задачи и визуализацию полученных результатов.

Разумеется, при исследовании биомедицинских проблем встречаются процессы, для математического описания которых используются аппараты обыкновенных дифференцированных уравнений (ОДУ), систем алгебраических нелинейных уравнений, разностные отображения, теории бифуркций, хаоса и порядка. Примеры успешного использования таких математических аппаратов представлены в [1] для прогнозирования развития болезни, в [2–5] — для решения задач нелинейной динамики в биологии, химической кинетике и др.

Кроме того, при изучении некоторых медицинских процессов необходимо численно решать жёсткие системы ОДУ, например, при моделировании протекания химических реакций, что представляет собой самостоятельную проблему, которой посвящена обширная литература [6–8].

Развитие численных методов решения задач МСС началось с проблем газодинамики (обтекание тел, спускаемых в плотных слоях атмосферы, точечного взрыва). Затем при помощи этих методов решались задачи физики плазмы, МДТТ и мн. др. Известно, что некоторые математические методы развивались под влиянием биомедицинских проблем, например, методы ма- тематической статистики, уравнение Воль-терра, нелинейные разностные отображения, теория хаоса и порядка, конечные автоматы, нейросети, методы решения жёстких ОДУ.

Постановки биологических и медицинских задач, которые приводят к необходимости численного решения систем уравнений в частных производных (МСС, уравнения параболического и эллиптического типа), появились относительно недавно. Они представлены в работах [9–11]. Для численного решения таких задач использовались методы, ранее применявшиеся для решения задач газогидродинамики [12–16]. Реологические соотношения для биологических сплошных сред разрабатывались в работах [17–19]. Круг рассматриваемых задач из этой области достаточно широк. Так, в работах [10, 11] рассматривались биохимическая и электрокардиографическая модели инфаркта миокарда, анализ которых выявил механизм его формирования и некоторые закономерности его течения. Сопоставление результатов расчётов с клиникой острого инфаркта миокарда выявило, что они позволяют отличить инфаркт миокарда лёгкого клинического течения от инфаркта тяжёлого клинического течения. Тяжесть острого инфаркта, как показали вычислительные эксперименты, определяется скоростью некротизации.

Механическая (в рамках МСС) модель сердца рассматривалась в работах [20–22]. Подобное рассмотрение представляет большой интерес для клинической практики, так как изучение работы сердечной мышцы в условиях частичного изменения её механических характеристик позволит понять, как изменится функционирование всего сердца и кровеносной системы в условиях частичного некроза. Однако разработка полной физико-математической трёхмерной модели, учитывающей электрохимические эффекты, их связь с напряжениями в сердечной мышце и численное решение соответствующей динамической пространственной задачи в рамках МСС, по-видимому, является перспективной проблемой.

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

ТРУДЫ МФТИ, 2009, Том 1, № 1 водится в работах [23, 24]. Основываясь на численном эксперименте с использованием математической модели Мак’Аллистера–Нобла–Тасена, описывающей динамику распространения импульсов в волокнах Пуркинье, авторы этих работ обнаружили различные режимы распространения таких импульсов. В частности, оказалось, что при различных режимах возбуждения возможно появление таких явлений, как отражения, автоколебания, аннигиляция, солитоноподобные режимы. Интересно, что подобные явления имеют место в нервных волокнах, возбуждение импульсов в которых описывается в рамках известной математической модели Ходжкина–Хаксли [25, 26].

Описание простейших математических моделей работы систем кровообращения и сердца можно найти, например, в [27, 28]. Функции кровеносной системы человека, которая состоит из малого и большого кругов кровообращения, очень важны и разнообразны, поэтому их моделирование, как в нормальных, так и в патологических условиях, представляет одну из важнейших задач медицины.

По-видимому, на сегодняшней день наиболее адекватными представляются динамические модели пульсирующих течений несжимаемой жидкости в системе растяжимых трубок. В работах [29–31] использовалась квазиодномерная гидравлическая модель несжимаемой жидкости в деформируемом кровеносном сосуде переменного сечения, обобщённая на случай иерархической ветвящейся системы кровеносных сосудов факториальной структуры. Подобная же иерархическая (или сетевая) квазиодномерная модель использовалась и для математического описания работы дыхательной системы на участке трахея — бронхи [31, 33].

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

В [30] рассматривался нестационарный квазипериодический режим кровообраще- ния головного мозга. С учётом взаимодействия с большим кругом кровообращения в этой работе рассчитана динамика изменения параметров крови в различных отделах головного мозга, подтверждена гипотеза о связи изменений объёма крови в процессе сердечного цикла с движением спинномозговой жидкости, возрастания амплитуды пульсовой волны, распространяющейся в сосудах Виллизиева круга и масштабных артерий в грудной и брюшной полости.

По инициативе нейрохирургов Главного военного клинического госпиталя им. Н.Н. Бурденко и Института скорой помощи им. Н.В. Склифосовского была поставлена задача о расчёте последствий черепно-мозговых травм. Надо сказать, что в рамках модели сплошной среды реальной представляется постановка задач о расчёте травм костей, грудной клетки, суставов, появлении гематом в мягких тканях тела, что также рассматривалось в работе. Экспериментальным данным, описывающим последствия черепно-мозговых травм, посвящены хорошо известные в нейрохирургии работы [34, 35]. Изучение областей повреждения мозга при различных динамических воздействиях на него путём математического моделирования в рамках МСС может дать важную информацию о последствиях травм. В газете «Известия» от 13.08.2003 опубликована статья «Бомба для мозгов», в которой утверждается, что в настоящее время составлена трёхмерная компьютерная карта мозга человека, над которой в течение нескольких лет в шести странах мира работали 20 тыс. нейрохирургов и невропатологов. Данные для этой карты получены в результате сканирования мозга 7 тыс. человек. Идея такой деятельности крайне важна, поскольку теперь с высокой степенью точности известно, какие области мозга управляют теми или иными его функциями. Понятно, что, определив области повреждения мозга при различных динамических нагрузках на него (или при черепно-мозговых травмах), мы сможем определить те области функции организма, которые могут оказаться поврежденными или измененными. В нейрохирургической практике хорошо известен следующий факт: области поражения мозга при черепно-мозговых травмах не всегда совпадают с обла- стями, прилежащими к месту удара. Примером этому является феномен «противо-удара»: при ударе затылком область повреждения мозга локализуется в лобной части головного мозга. Объяснение этому явлению можно дать только путём проведения численного исследования сложнейших волновых процессов, образующихся в неоднородной механической конструкции, которая представляет собой систему череп–мозг (разумеется, речь идёт только о механической системе, о принципах функционирования мозга, как об управляющей системе, речь не идет). Кроме того, механико-математические модели реакции рассматриваемой механической системы на динамические нагрузки можно использовать как инструмент для изучения различных параметров травмы (место приложения удара, сила удара, геометрия головы, атрофия мозга и др.) на вероятную степень риска. Эти задачи ставились в работах [36, 37]. В [38] численно изучался процесс распространения ударного импульса через 5-слойную биоконструкцию черепа человека, которая, как показали расчёты, оказалась идеальной противоударной защитой для мозгового вещества при естественных динамических воздействиях (разумеется, при действии боевых поражающих средств необходимо использовать дополнительные средства защиты головного мозга).

В работе [39] приводятся примеры расчётов многослойных защитных конструкций, которые могут использоваться для защиты как техники, в которой находится человек, так и тела (бронежилеты, каски, др.). Разработанные авторами вычислительные методы и алгоритмы позволяют изучать также динамические процессы, происходящие в теле человека, при воздействии динамических нагрузок. Такая задача оказывается весьма актуальной, поскольку появляется возможность априорного предсказания областей и степеней поражения тела при каких-либо опасных нагрузках. Области приложения подобных задач очевидны (кроме медицинской): спортивные травмы, травмы при производстве опасных видов работ в экстремальных условиях (например, горные или взрывные работы), боевых и полицейских операциях, автомобильных, авиационных, иных техногенных катастрофах.

К проблеме математического моделирования травматологических процессов относится и задача о залечивании ран. Численному изучению этого процесса посвящены работы [40, 41], в которых получено количественное описание динамики залечивания резаной раны кожного покрова человека.

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

Перспективным направлением вычислительной медицины представляется компьютерная реализация виртуальных хирургических операций и предсказания их последствий. Разумеется, это очень сложное направление, которое только начинает появляться, а постановка некоторых задач до конца не ясна (например, моделирование работы хирурга с помощью скальпеля). Однако реализация некоторых виртуальных операций является реальной задачей. Так, в работе [43] представлено численное моделирование операций литотрипсии (дробление почечных камней акустическим волнами, инициируемыми искровым разрядом или лазерным импульсом). Цель таких исследований — найти режимы работы литотриптора (длительность и интенсивность импульса, количество импульсов), при которых фрагменты разрушенного камня были бы достаточно малыми для их выведения из организма естественным путём. Для этого численно исследовалась картина распространения акустического импульса в теле и в камне, а также решалась задача его разрушения.

Другой пример — моделирование офтальмологической операции экстракции (удаления) катаракты. Суть операции в том, чтобы с помощью лазера или ультразвукового факоэмульсификатора разрушить помутневший хрусталик (точнее его плотное ядро) так, чтобы не повредить сетчатку и роговицу глаза, и вывести мутные хрусталиковые массы. Эта задача условно подразделяется на три части: первая — расчёт импульсного воздействия на хрусталик, вторая — распространение акустического импульса в стекловидном теле до сетчатки и расчёт динамического воздействия на неё (поскольку в результате взаимодействия импульса с сетчаткой последняя может расслаиваться), третья — вымывание мутных хрусталиковых масс из передней камеры глаза. Численное решение последней задачи, которое сводится к решению уравнения Пуассона, в двухмерном случае рассматривалось в работе [44], причём с использованием как прямоугольных, так и треугольных расчётных сеток. Расчёты проводились с целью определения застойных зон в передней камере глаза при вымывании мутных масс и оптимизации рабочих режимов хирургических инструментов. Проблема распространения импульса через хрусталик и стекловидное тело к сетчатке представляет сложную динамическую задачу, поскольку глаз является неоднородной механической системой с рядом поверхностей раздела сред. Численное её решение представлено в работе [45]. Моделирование этого процесса сводится к решению динамической системы уравнений гиперболического типа (уравнений МДТТ). Исследования выполнялись в тесном контакте с сотрудниками Центра микрохирургии глаза им. академика С.Н. Федорова, который является соавтором этих работ.

Важным приложением вычислительной медицины являются проблемы предсказания динамики развития онкологических заболеваний, то есть развития опухолей, в том числе с учётом кровообращения. Для их численного решения используются уравнения гидродинамики, жёсткие системы ОДУ, уравнения типа реакция–диффузия. Численное моделирование этих процессов проводится в работах [46, 71].

Отметим, что с помощью нелинейных уравнений параболического типа (реакция–диффузия) проводится также и численное моделирование процессов структурообразования в активных биосре- дах, колониях бактерий, микроорганизмов (например, Esherichia coli, Distyostelium discoicleum). Этим задачам посвящены работы [47–50], в которых численно решаются двух- и трёхмерные динамические задачи об образовании таких структур. Задачи структурообразования при свертывании крови и тромбообразовании рассматривались в работах [50, 51].

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

Для решения задач гидродинамики используются уравнения Эйлера или На-вье–Стокса, для задач МДТТ — система уравнений Ламэ (в напряжениях или деформациях).

Эти системы можно найти, например, в монографиях [52–54].

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

Численные методы решения эллиптических уравнений (уравнения Лапласа или Пуассона) подробно рассматриваются в известных монографиях [56–59]. Среди методов, представленных в литературе для решения этого класса задач, отметим наиболее широко используемые: в вычислительной практике — метод верхней релаксации, трёхслойный метод Чебышева, метод Дугласа–Ренфорда, попеременно-треугольный метод, многосеточный метод Р.П. Федоренко. В работе [60] предлагаются монотонные разностные схемы для численного решения этого типа уравнений на нерегулярных сетках. С использованием такого подхода рассматривалась задача о гидродинамическом течении в передней стенке глаза при экстракции катаракты.

Широкий класс процессов в биологии и медицине моделируется при помощи нелинейных уравнений параболического типа (реакция–диффузия). Исследованию свойств решений уравнений этого типа посвящены работы [4, 61, 62 и др.]. Для численного решения таких уравнений обычно используются двух- и трёхслойные разностные схемы типа Кранка–Николь-сона, Нумерова, Толстых с последующим проведением итераций по нелинейности. Эти вычислительные методы описаны в [12, 58, 59].

По-видимому, наибольшие трудности при численном решении рассматриваемых задач вызывают системы уравнений в частных производных гиперболического типа (МСС). Разработке вычислительных методов решения систем этого типа посвящены работы [63–67, 69, 70]. Наиболее эффективными для численного решения такого класса задач являются сеточно-характеристические методы [65–67, 69, 70]. Эти методы оказываются не только наиболее точными, но и позволяют наиболее корректно строить вычислительные алгоритмы на границах области интегрирования и границах раздела сред. Одним из наиболее важных показателей качества численного решения является его близость к точному вблизи областей с большими градиентами решений. Это вызвано немонотонным, в одномерном случае, поведением численных решений, получаемых с помощью так называемых немонотонных схем (в линейном случае это, в соответствии с теоремой С.К. Годунова, схемы порядка точности выше первого) и «размыванием» разрывов схемами первого порядка точности. По этой причине многие работы, посвященные этой теме, направлены на регуляризацию решений уравнений гиперболического типа [59, 60, 66, 70]. В работе [70] приводятся определения и критерии монотонности разностных схем для этого типа уравнений в частных производных. В приложении к задачам медицины эти методы использовались при изучении процессов литотрипсии, офтальмологических операций, черепно-мозговых травм, деформации сердечной мышцы [22, 37, 43, 44, 45].

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

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

На рис. 1 представлена схема черепномозговой системы, которая выбиралась для моделирования последствий черепномозговых травм. Расчётные сетки (в горизонтальной проекции) как четырёхугольные, так и треугольные показаны на рис. 2 (соответственно a, b, c).

Рис. 1. Схема черепно-мозговой системы, ко- торая выбиралась для моделирования послед- ствий черепно-мозговых травм

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

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

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

Рис. 3. Области максимальных сжимающих и растягивающих напряжений

Рис. 4. Сравнение расчётной и полученной при томографических исследованиях областей поражения головного мозга

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

Сравнение расчётной и полученной при томографических исследованиях областей поражения головного мозга видно на рис. 4 (томография проведена в Главном военном клиническом госпитале имени Бурденко).

На рис. 5 показана схема глаза человека, которая использовалась при моделировании динамических процессов, происходящих при проведении офтальмологических операций по экстракции катаракты при помощи лазера или факоэмульсифика-тора (ультразвуковой инструмент).

Рис. 6. Расчётные сетки, поля скоростей и области возможных поражений глаза при лазер- ном разрушении хрусталика

Рис. 7. Расчёт давления в глазу человека при проведении лазерной операции и расчётная сетка

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

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

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

Схематичное изображение сердца человека (в разрезе) представлено на рис. 9.

Рис. 9

Расчётная сетка, используемая при мо- делировании процесса начала сжатия сердца (систолы) и начального выталкивания крови, а также поля скоростей в сердце, дана на рис. 10.

Рис. 10. Расчётная сетка, используемая при моделировании процесса начала сжатия сердца (систолы) и начального выталкивания крови, а также поля скоростей в сердце

Столкновение двух импульсов в волокнах Пуркинье в автоколебательном режи- ме в четыре момента времени показаны на рис. 11 (в расчётах были получены также режимы аннигиляции и солитоноподобные режимы).

Рис. 11. Столкновение двух импульсов в волокнах Пуркинье в автоколебательном режиме в четыре момента времени

Рис. 13. Структура, сформированная пассивными клетками E-coli в бактериальной колонии, отражающая механизмы межклеточной

Динамика залечивания кожной раны (пореза) представлена на рис. 12. Здесь показаны трёхмерные картины распределения плотности коллагена, являющиеся «подложкой», на которой растут клетки кожи, в два момента времени: t = 0 ( р ^ 1 в области поражения и р = 1 на неповрежденной коже; ρ — плотность коллагена); t = 40 дней (плотность коллагена заметно повышается внутри раны вследствие процесса хемотоксиса). Величина плотности коллагена представлена светлой частью поверхностей, изображённых на рисунке.

Рис. 12. Динамика залечивания кожной раны (трёхмерные картины распределения плотности коллагена в два момента времени)

регуляции

Рис. 14. Распределения концентрации ЦАМФ и её дальнейшее «закручивание» в спираль в двухмерном случае

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

Рис. 15. Формирование трёхмерной волны в активной среде с распределенным вдоль оси коэффициентом возбудимости

Спиралевидное двух- и трёхмерное распространение волн возбуждений в нелинейной активной среде, в частности, в колонии одноклеточных организмов Dictyostelium discoideum, представлено на рис. 14 и 15. Особенность этих микроорганизмов состоит в том, что благодаря своей подвижности и хемотоксису в определённых внешних условиях (при недостатке субстрата) амебы объединяются в псевдоорганизм, который благодаря согласованному движению составляющих его клеток ведёт себя, как одно целое. В активной среде, сформированной этими клетками, образуется спиральная волна ЦАМФ (циклический аденозинмонофосфат) с центром в точке агрегации.

Для численного решения также использовалась система уравнений параболического типа как в двух, так и в трёхмерном случаях. Начальные распределения концентрации ЦАМФ и её дальнейшее «закручивание» в спираль, в двухмерном случае показаны на рис. 14. Формирование трёхмерной волны в активной среде с распределенным вдоль оси коэффициентом возбудимости показано на рис. 15. Понятно, что вихреобразная волна в передней части закручивается быстрей, чем в задней из-за изменения этого коэффициента по координате, что приводит к её неравномерной закрутке.

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

Автор выражает благодарность за внимание к циклу работ, о котором шла речь в статье, действительному члену РАН О.М. Белоцерковскому, члену-корреспонденту РАН А.С. Холодову, профессору Г.Г. Малинецкому, профессору А.А. Полежаеву, профессору А.И. Лобанову, профессору Г.Ю. Резниченко.