ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ РЕАКЦИОННО-ДИФФУЗИОННЫХ ПРОЦЕССОВ В ЗАДАЧЕ ИССЛЕДОВАНИЯ КОЛЛЕКТИВНОГО ПОВЕДЕНИЯ МИКРООРГАНИЗМОВ

Математика, Прикладная математика

УДК 51-76:519.63

П.А. Юшкевич, Д.А. Коршик, А.Г. Масловская

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ РЕАКЦИОННО-ДИФФУЗИОННЫХ ПРОЦЕССОВ В ЗАДАЧЕ ИССЛЕДОВАНИЯ КОЛЛЕКТИВНОГО ПОВЕДЕНИЯ МИКРООРГАНИЗМОВ

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

Работа выполнена в рамках проекта «Математическое моделирование процессов типа «реакция-диффузия» в приложениях к биологическим и физическим системам» по государственному заданию Министерства науки и высшего образования РФ высшим учебным заведениям, код проекта: 1.13421.2019/13.2.

Введение

На современном этапе разработка и развитие математических моделей, описывающих реакци-онно-адвективно-диффузионные системы, представляют как фундаментальный научный, так и практический интерес. В математической форме подобные системы могут быть описаны с помощью аппарата обыкновенных дифференциальных уравнений и уравнений с частными производными [1-9]. Примерами приложений, использующих нелинейные эволюционные уравнения в частных производных, являются уравнение Шредингера [4, 9], описывающее вероятность состояния и динамику частиц в квантовой механике; уравнение Нагумо, моделирующее распространение нервных импульсов в биологических системах [4, 7]; уравнение Блэка - Шоулза для определения стоимости азиатских и европейских опционов [5]; формализм теории Ландау для представления фазовых переходов, сверхпроводимости и процессов волоконной оптики [3]; уравнения «реакция-диффузия-конвекция (адвекция, дрейф)» в физике и химии (для описания химических реакций и диффузионных явлений) [4, 8, 9], в теории тепломассопереноса [1], в геологии, в задачах экологического прогнозирования и мониторинга окружающей среды [6].

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

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

В данном междисциплинарном подходе бактерия рассматривается не как простейший и примитивный организм, а вводится описание коллективного поведения этого сообщества, обладающего «quorum-sensing» - чувством кворума (способностью коллективно действовать на внешние возбудители) [10-12]. Особую актуальность в биологии это направление приобретает в связи со способностью бактерий адаптироваться к воздействию антибиотиков и необходимостью человека прогнозировать и управлять реакцией сообщества бактерий на внешние воздействия [12]. В аспекте математического и компьютерного моделирования данный подход требует всестороннего анализа и развития (разработка модифицированных моделей, учитывающих наличие обратной связи и эффект запаздывания, разработка вычислительных схем реализации моделей, проектирование системы компьютерного моделирования с учетом набора управляющих параметров в данной предметной области, интерпретация и анализ результатов).

Механизм, описывающий «чувство кворума» при коммуникации колоний бактерий, математически может быть описан в виде системы обыкновенных дифференциальных уравнений [13-17]. Базовые математические модели допускают формализацию с учетом законов популяционного роста, наличия отрицательной обратной связи, присутствия эффекта запаздывания и др. Последнее может быть ассоциировано с активацией в системе фермента «лактоназа» [14]. Однако во многих естественных ситуациях необходимо принять во внимание пространственное распределение ячеек и получающееся неоднородное распределение их отклика на внешнее воздействие. В зависимости от того, как моделируется бактериальная динамика, определяется пространственное распределение питания, которое не только зависит от пространственного распределения ячеек, но и само влияет на бактериальный рост. Бактерии и химически активные ферменты задаются в модели с помощью удельных весов и концентраций. Модифицированные математические модели описываются в терминах уравнений вида «реакция-диффузия», формализованных с помощью уравнений с частными производными параболического типа [17].

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

Постановка задачи моделирования

Математическую модель сформулируем в виде одномерной начально-граничной задачи для уравнения типа «реакция-диффузия» [17]:

^^ = DAU(x,t)-yU(x,t) + F(x,U), 0<x<L, 0<t<~t, (1)

U(0,t) = 0, U(L,t) = 0, (2)

U(x, 0) = 0, (3)

где U(xj) - концентрация аминокислоты в точке х в момент времени /, моль/л; D - коэффициент диффузии аминокислоты, мкм2/час; у - абиотическое (обусловленное внешними факторами) снижение скорости образования аминокислоты, 1/час; а - базовая производительность аминокислоты без масштабирования, моль/(л-час); функция, описывающая источник генерации аминокислоты и реакционное слагаемое

F(x,U) = ¿/(С/)хмJ (*) > f(U) = « + РU2-5I(U™ +U25), моль/(ед,час);

X\\x-xk\\ ~~ дельта-функция Дирака для каждой ячейки фиксированного размера.

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

Вычислительная схема

Сконструируем вычислительную схему для реализации математической модели в постановке (1)-(3) на основе идеологии метода конечных разностей. Введем в рассмотрение конечно-разностную сетку:

«А = {*/ = (& - 1) К / = = 0 - 1) X, ] = йм).

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

явная схема порядка точности 0(//2,т) - проста в реализации, но обладает только условной устойчивостью;

неявная схема порядка точности 0(//2,т) - требует решения системы линейных алгебраических уравнений на каждом шаге по времени, абсолютно устойчива;

схема Кранка - Николсон порядка точности 0(//2,т2) - явно-неявная схема, обладающая безусловной устойчивостью.

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

(М/И)(у+1) И

-(¿+1,7+1)

Рис. 1. Конечно-разностный шаблон пятиточечной конечно-разностной схемы.

ь2 _ 1—1

. 4£>Т „

и Г1 +

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

Известно, что неявные схемы лишены указанного недостатка (потери монотонности), в связи с чем в данной работе для численной реализации модели в постановке (1)-(3) предлагается использовать двухслойную неявную схему второго порядка точности 0(//2,т2), конечно-разностный шаблон которой представлен на рис. 1.

Поставим в соответствие непрерывной задаче (1)-(3) ее дискретный аналог. В результате преобразований для всех внутренних узлов области будем иметь:

и£=411{ -иГ1+2Е/т,

где / = 2,Ж-1, 7 = 1,М-1.

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

Задача замыкается заданием начального условия

и] = 0 для всех 1 = (6)

и краевых условий

и( = 0, и^ =0 для всех у = \\М . (7)

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

Решение модельной тест-задачи

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

Переменные и параметры математической модели

Переменная/ параметр Описание Значение / единица измерения

Щх, 0 концентрация аминокислоты моль/литр

Т время наблюдения процесса 10 час

в коэффициент диффузии аминокислоты 320 мкм2/час

функция, описывающая источник генерации аминокислоты и реакционное слагаемое ед. актив, ферм, /литр (или катал/литр)

У абиотическое (обусловленное внешними факторами) снижение скорости образования аминокислоты у = 0.005545, 1/час

а базовая производительность аминокислоты без масштабирования а = 4.37 • 10"4, моль/(литр-час)

Р производительность аминокислоты (учет обратной связи) без масштабирования Р = 4.37 • 10~3, моль/(литр-час)

дельта-функция Дирака для каждой ячейки фиксированного размера безразмерная величина

ил пороговое значение концентрации аминокислоты (в функции источника) 70 нмоль/л

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

рия Фурье (4): /Лс = —^«10 час.

Таким образом, при временах X » 10 час. можно наблюдать поведение системы в стационарном режиме, при времени порядка 10 час. система будет находиться в «динамическом» режиме. Установим параметры управления вычислительным процессом: т=1 час, //=0.1 мкм. Результат эксперимента - модельное распределение концентрации аминокислоты как функция координаты и времени -представлено на рис. 2.

При переходе в стационарный режим (время наблюдения 100-200 час.) максимальное значение концентрации аминокислоты будет наблюдаться в центральной ячейке и соответствовать 7-10"3 моль/л.

Рис. 2. Распределение концентрации аминокислоты в течение всего времени наблюдения процесса (три симметричные ячейки с линейным размером 1 мкм).

Заключение

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

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

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