CC BY
119
СПИСОК ЛИТЕРАТУРЫ
1. Экономико-математические методы и прикладные модели / Под ред. В.В. Федосеева. - М.: ЮНИТИ, 2001. - 391 с.
2. Занг В.-Б. Синергетическая экономика. Время и перемены в нелинейной экономической теории. - М.: Мир, 1999. - 334 с.
3. Хакен Г. Синергетика. Иерархии неустойчивостей в самоорганизующихся системах и устройствах. - М.: Мир, 1985. - 423 с.
4. Григорьев В., Козловских А., Ситникова О. Динамическая модель фьючерсного рынка // Рынок ценных бумаг. - 2004. -№ 24 (279). - C. 42-44.
5. Григорьев В.П., Козловских А.В., Марьясов Д.А. Качественное исследование системы дифференциальных уравнений модели
динамического хаоса и корреляция особых точек с трендами // Известия Томского политехнического университета. - 2006. -Т. 309. - № 2. - С. 12-17.
6. Шустер Г. Детерминированный хаос. - М.: Мир, 1998. - 240 с.
7. Мерфи Дж. Технический анализ фьючерсных рынков: теория и практика. - М.: Сокол, 1996. - 412 с.
8. Берже П., Помо И., Видаль К. Порядок в хаосе. О детерминистском подходе к турбулентности. - М.: Мир, 1991. - 368 с.
Поступила 10.06.2009 г.
УДК 517.9;577.3.01;577.38
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ОДНОМЕРНОЙ ПОПУЛЯЦИОННОЙ ДИНАМИКИ С НЕЛОКАЛЬНЫМ ВЗАИМОДЕЙСТВИЕМ
А.В. Борисов*, Р.О. Резаев, А.Ю. Трифонов, А.В. Шаповалов*
Томский политехнический университет
E-mail: [email protected] Томский государственный университет*
Численными методами исследовано влияние нелокальных эффектов на динамику популяции микроорганизмов в рамках диффузионной модели с квадратично-нелинейным нелокальным взаимодействием. Основное уравнение модели обобщает известное уравнение Фишера-Колмогорова-Петровского-Пискунова. Для описания нелокального взаимодействия использовались ядра в виде гауссова и равномерного распределений. Показано, что в обоих случаях образуются локальные максимумы численности популяции, что указывает на формирование популяционной структуры.
Ключевые слова:
Популяционная динамика, численные методы, модель Фишера, популяционные структуры. Key words:
Population dynamics, numerical techniques, Fisher model, population structure.
Введение
Популяционная динамика микроорганизмов дает пример нелинейной системы, в которой могут возникать эффекты самоорганизации, обусловленные коллективным поведением большого числа особей в условиях влияния управляющих внешних воздействий. Эффекты самоорганизации проявляются в формировании неоднородных пространственно-временных структур (популяционных агрегаций или паттернов) [1-3], популяционных волн и других подобных явлений [4]. Понимание закономерностей возникновения и динамики роста бактериальных структур важно для выявления основных механизмов контроля на начальном уровне зарождения и развития бактериальных инфекций в медицинской практике [5]. Систематические исследования популяций бактерий, проведенные микробиологами на протяжении прошлого века, привели к представлению о бактериальной культуре как о целостной единой системе [6, 7].
В работе численными методами исследовано влияние нелокальных эффектов на динамику по-
пуляции микроорганизмов в рамках диффузионной модели с квадратично-нелинейным нелокальным взаимодействием, обобщающей известную модель Фишера-Колмогорова-Петровского-Пис-кунова (ФКПП). Для описания нелокального взаимодействия использовались ядра двух типов: в виде гауссова и равномерного распределений. Показано, что в обоих случаях образуются локальные максимумы численности популяции, что указывает на формирование популяционной структуры.
1. Популяционная модель с нелокальной
квадратичной нелинейностью
В качестве базовой модели пространственно распределенной популяции, следуя традиционному подходу в динамической теории популяций [6-8], выберем модель, предложенную независимо Р. Фишером [9] и А.Н. Колмогоровым, Н.Г. Петровским, Н.С. Пискуновым [10]. Для простоты будем считать пространство одномерным. Данное предположение может быть реализовано в экспериментальных условиях, если область, в которой взаимодей-
ствуют особи микробной популяции, представляет собой узкую и длинную трубку. Вдоль трубки осуществляются процессы переноса, а в любом ее поперечном сечении происходит полное внутреннее перемешивание [7]. В модели Фишера-Колмогоро-ва-Петровского-Пискунова динамика зависящей от времени / и пространственной координаты х кинетической переменной и(х,0 (биомасса или число организмов данного вида в единице объема, в одномерном случае - на единицу длины) описывается эволюционным уравнением, в котором учитываются пространственная диффузия с постоянным коэффициентом диффузии Б, процесс производства бактерий с постоянным темпом роста (мальтузианским параметром) а и квадратичными по плотности потерями с коэффициентом Ь, обусловленными конкуренцией за ресурс. При сделанных предположениях уравнение модели ФКПП имеет вид
ди(х,0 = п^^ит + т(Хг0-Ъи\х,ъ (1)
dt
dx2
В модель ФКПП включены основные биологические закономерности, определяющие популя-ционную динамику, однако такие факторы, как неоднородность популяции по составу, в частности мутация, таксис (подвижность особей), влияние продуктов метаболизма на рост популяции, влияние неоднородности окружающей среды, в модели не учитываются. Модель ФКПП не достаточна для объяснения формирования популяционных агрегаций. Одно из направлений исследований, имеющих целью объяснить это явление, состоит в модификации модели ФКПП посредством включения в уравнение модели дополнительных слагаемых, ответственных за указанные выше дополнительные факторы популяционной динамики.
Следуя [1-6], рассмотрим обобщенную модель популяционной динамики ФКПП, включающую нелокальное квадратично-нелинейное взаимодействие. Основное уравнение модели записывается в виде
ди(х,/) 52и(х,^ ч ,
у ' '■ = п—+ а(х,0и(х,0 -
dt
dx2
-kU (x, t) J b(x, y, t)u (y, t)dy.
(2)
Здесь коэффициент прироста а(х,0 и коэффициент конкурентных потерь Ь(х,у,1) (функция влияния) являются функциями пространственной координаты х и времени ¡, что позволяет расширить область применения модели (1) на случай пространственно неоднородных условий воспроизводства популяции и потерь, создаваемых внешними факторами (например, различной освещенностью, градиентами температуры и т. д.), к - параметр нелинейности. Будем считать величины, входящие в ур. (2), безразмерными.
Уравнение (2) является примером уравнения типа «реакция - диффузия», или уравнением
РД-типа. Нелокальность в уравнении учитывает конечность пространственной области, окружающей особь, в которой она взаимодействует с соседними особями популяции, пространственное распределение продуктов метаболизма, таксис, неоднородность распределения субстрата и другие факторы окружающей среды в области, занимаемой популяцией.
2. Разностные схемы для нелокальной
популяционной модели
Исследуем уравнение (2) численными методами, используя соответствующие разностные схемы и алгоритмы численного решения. Для уравнения (2) с коэффициентами общего вида неизвестны точные решения и законы сохранения, что ограничивает возможности оценки адекватности построенных численных решений. Данная оценка проводилась посредством сопоставления результатов расчетов, полученных на основе явной и неявной разностных схем и метода конечных элементов.
Разностная схема строится стандартным образом (см., например, [11]). Введем равномерную сетку по времени / с шагом т, /=/т, /=0,...,/, на временном отрезке [0,7], Т=1т, а по координате х - с шагом к, х=к,}, /=-/, -/+1,...,0,.../, в симметричной области [-Х,Х], Ь=1к, 1,1 — целые числа. В узле (хД) сетки функцию и(х,0 обозначим как и(х;Д)=и/, аналогично а(х;Д)=а/, Ь(х,,у/,/;)=Ь;7'. Для производных в узле (х;Д) используем разностные аппроксимации
и, (х, ^) * и' = (и-+1 - и' )т~\
ихх(х!,) *Ки) = и+1 - Ц + и
Интеграл в уравнении (2) в узле (х,1) аппроксимируем выражением
г -1
I Ъ(х,у,^и(у,^йу*В) = ^Ъ)1и',Ь.
Простейший алгоритм численного решения ур. (2) реализуется на основе двухслойной разностной схемы
'= uii +т (DA xxuij + ajuj — ku i В i).
(3)
Условия применимости схемы (3) рассмотрены в [11] (ч. III, гл. 1, § 4). Схема (3) устойчива при выполнении условия rD/h2<1/2, что является ее недостатком, т. к. ограничение на шаг по времени приводит к необходимости увеличивать число шагов I.
Неявные разностные схемы использовались в следующих двух формах:
u:1 —r(DAu' + d+ V,+') = uii —TKui В' (4)
XX 1 / 1 ' 1 II
и'+ -т (БЛххи';1 + а'+ 1и'+ 1-ки'+ 'Б')=и). (5)
Схема (4) является двухслойной и решалась методом прогонки [11. С. 45]. При этом расчетный
и
шаг по времени т выбирался таким образом, чтобы вид решения качественно не изменялся при дальнейшем уменьшении т. Численное решение на основе схемы (5) строилось с помощью метода Ньютона [11].
Начальные условия и(х,0) для построения численных решений выбирались в форме функций, локализованных в малой окрестности с центром в начале координат х=0. Решения строились на временном интервале ?е(0,Г0). За расчетное время Т0 ненулевые значения искомой функции и(х,0 не должны достигать левой и правой границ области пространственной переменной х, -Ь<х<Ь. Это соответствует нулевым граничным условиям на бесконечности и(х,0채=0. Отметим, что наиболее простой способ построения численных решений ур. (2) с указанными условиями дает явная разностная схема (3).
При построении численных решений ур. (2) необходимо задать функцию конкурентных потерь Ь(х,у,1). Поскольку вид этой функции априори неизвестен, в данной работе рассмотрены два модельных примера: гауссово и равномерное распределения, что позволяет выявить характерные свойства динамики роста.
Сформулированная выше математическая постановка задачи согласуется с типичными условиями экспериментов по наблюдению динамики роста культуры бактерий. Согласно этим условиям, в ограниченную окрестность центра области (подложки) с нанесенным питательным веществом вносится небольшое количество культуры, дальнейший рост которой управляется законами попу-ляционной динамики. Наблюдение ведется до контакта области, занимаемой популяцией, с граница-
ми подложки. Численные решения ур. (2) проводились с помощью алгоритмов, построенных на основе разностных схем (3)-(5). Результаты расчетов дополнительно проверялись с помощью расчетов по методу конечных элементов с использованием программы ANSYS.
3. Решение для гауссова распределения
и асимптотических нулевых граничных условий
Рассмотрим пример численного решения уравнения (2) с гауссовой функцией влияния Ь(х,у,1). Начальное распределение численности популяции и функцию Ь(х,у,1) зададим в виде гауссовых распределений
и( х, 0) = ы0( х) = /0ехр( -х2/ ст02), (6)
Ь(х,у,() = ехр(-( х-у)2/а2). (7)
Начальное распределение (6) локализовано в окрестности начала координат х=0, область локализации определяется значением параметра а0 - дисперсией функции и0(х), а - дисперсия функции Ь(х,у,1).
Результаты численного решения уравнения (2) при условии (6), (7) приведены на рис. 1, из которого видно, что начальное возмущение не достигает пространственных границ расчетной области, что соответствует асимптотическим нулевым граничным условиям и(х,О채=0.
В динамике роста популяции можно выделить следующие стадии.
На первой стадии начальное распределение и0 вида (6) (рис. 1, а) преобразуется в симметричное распределение и1 с максимумом в начале координат х=0 (рис. 1, б). Дополнительно проведенные расче-
ты показали, что состояния иъ образующиеся на первой стадии при различных значениях амплитуды начального распределения /0 в (6), мало отличаются друг от друга.
На второй стадии центральный максимум распределения разделяется на два локальных максимума, которые движутся в разные стороны. Отметим, что амплитуды «вторичных» локальных максимумов распределения ип (рис. 1, в) меньше, чем амплитуда максимума состояния и1.
На третьей стадии в центре распределения ип (в начале координат) возникает локальный максимум, который, в свою очередь, разделяется на два локальных максимума меньшей величины, чем на второй стадии (рис. 1, г). Процесс образования центрального максимума и его разделения может продолжаться и далее. Анализ численных решений показывает, что амплитуды образовавшихся локальных максимумов со временем выходят на стационарный уровень, причем амплитуды максимумов, образовавшихся на более поздних стадиях, меньше, чем амплитуды максимумов, образовавшихся на более ранних стадиях.
4. Решение для равномерного распределения
и асимптотических нулевых граничных условий
Выберем функцию влияния Ь(х,у,1) в виде Ь(х, у, Г) = в (х - у + с/2) - в (х - у - с/2), (8)
где в(х) - ступенчатая функция Хевисайда; с -ширина прямоугольника. Численное решение ур. (2) с начальным условием (6) и функцией влияния вида (8) проиллюстрировано на рис. 2. Аналогично предыдущему примеру в динамике роста популяции выделены три стадии, первая из которых качественно совпадает с первой стадией динамики.
Вторая стадия существенно отличается от второй стадии предыдущего примера тем, что при разделении максимума распределения и1 образуется несколько максимумов, а возникновение последующих максимумов, в свою очередь, происходит при делении крайних максимумов. При этом смещение максимумов не наблюдается, а происходит лишь их сужение, причем минимальные значения стремятся к нулю. Из вида распределения ип на второй стадии (рис. 2, в) очевидно, что амплитуды максимумов модулируются функцией (огибающей локальных максимумов функции ип), медленно меняющейся по сравнению с изменением функции ип. Огибающая имеет максимум в начале координат.
Третья стадия (рис. 2, г) характерна тем, что минимумы, образовавшиеся между максимумами, при этом достигают нулевого значения. Этот процесс продолжается и образуется пространственно периодическая структура, период которой определяется значением параметра с. В отличие от третьей стадии предыдущего примера (рис. 1, г), амплитуды образовавшихся локальных максимумов со временем достигают стационарного значения.
Описанные выше процессы обусловлены балансом между протекающими процессами диффузии, конкурентных потерь с учетом нелокального взаимодействия и прироста популяции с темпом роста а. Процесс продолжается до достижения колонией бактерий границы занимаемой ею области. Отметим, что для конечной области количество образующихся максимумов ограничено.
По достижении границ пространственной области динамика роста колонии бактерий претерпевает следующие изменения. В случае симметричного относительно начала координат интервала и
периодических граничных условий функция u(x,t) на границах растет. С увеличением времени функция u(x,t) переходит в стационарное состояние: u(x,t)^J^ust(x). В случае несимметричного интервала и периодических граничных условий поведение колонии бактерий отличается лишь промежуточными стадиями, предельное стационарное распределение щ^) сохраняется при неизменных условиях роста популяции (пространственной области xe[-Z,Z] и параметров D,a,к,b(x,y,t)).
Заключение
Построены численные решения уравнения (2), соответствующие нулевым граничным условиям на бесконечности. Анализ численных решений, проведенных для различных начальных распределений u0(x) и двух видов функций конкурентных потерь Ь^у^), показал, что процесс эволюции проходит несколько характерных стадий (см. рис. 1, 2). В первом случае характерным является последовательное возникновение локальных максимумов в центре распределения и последующее их разделение на два локальных максимума, движущихся в противоположные стороны. Во втором случае возникает пространственно периодическая структура с одинаковыми по форме и параметрам пиками. Такая динамика указывает на процесс формирования популяционной агрегации. Аналогичные результаты получены в [4] для периодических гра-
СПИСОК ЛИТЕРАТУРЫ
1. Budrene E.O., Berg H.C. Complex patterns formed by motile cells of E. coli // Nature. - 1991. - V. 349. - № 2. - P. 630-633.
2. Ben Jacob E., Garik P. The formation of patterns in non-equilibrium growth // Nature. - 1990. - V. 343. - № 2. - P. 523-530.
3. Fujikawa H., Matsushita M. Bacterial fractal growth in the concentration field of nutrient // J. Phys. Soc. of Japan. - 1991. - V. 60. -№ 1. - P. 88-94.
4. Цыганов М.А., Бикташев В.Н., Бриндли Дж., Холден А.В., Иваницкий Г.Р. Волны в кросс-диффузионных системах - особый класс нелинейных волн // Успехи физических наук. -2007. - Т. 177. - № 3. - С. 275-300.
5. Fuentes M.A., Kuperman M.N., Kenkre V.M. Nonlocal interaction effects on pattern formation in population dynamics // Phys. Rev. Lett. - 2003. - V. 91. - № 158104. - P. 158104-1-158104-4.
6. Murray J.D. Mathematical Biology. I. An Introduction (Third Edition). - N.Y., Berlin, Heidelberg: Springer-Verlag, 2001. - 551 p.
ничных условий и функции влияния, вид которой отличается от (7) и (8). Эти закономерности и особенности связаны с пространственной неоднородностью лимитирующих рост факторов (конкуренцией за субстраты, продуктами метаболизма), с взаимодействием различных видов или групп клеток в процессе роста [1-3, 5-6].
Выводы
1. Численными методами исследовано влияние нелокальных эффектов на динамику популяции микроорганизмов в рамках диффузионной модели с квадратично-нелинейным нелокальным взаимодействием.
2. Показано, что основное уравнение модели обобщает известное уравнение Фишера-Кол-могорова-Петровского-Пискунова. Для описания нелокального взаимодействия использовались ядра в виде гауссова и равномерного распределений.
3. Установлено, что при решении уравнения типа «реакция - диффузия» образуются локальные максимумы численности популяции, что указывает на формирование популяционной структуры.
Работа выполнена при частичной финансовой поддержке гранта Президента РФ НШ-871.2008.2; АВЦП Министерства образования и науки РФ № 2.1.1/3436.
7. Романовский Ю.М., Степанова Н.В., Чернавский Д.С. Математическая биофизика. - М.: Наука, 1984. - 304 с.
8. Ризниченко Г.Ю., Рубин А.Б. Математические модели биологических продукционных процессов. - М.: Изд-во МГУ, 1993. - 302 с.
9. Fisher R.A. The wave of advance of advantageous genes // Annual Eugenics. - 1937. - V. 7. - № 3. - P. 255-369.
10. Колмогоров А.Н., Петровский Н.Г., Пискунов Н.С. Исследование уравнения диффузии, соединенной с возрастанием вещества, и его применение к одной биологической проблеме // Бюллетень МГУ. Сер. А. Математика и Механика. - 1937. -Т. 1. - № 6. - С. 1-16.
11. Самарский А.А., Гулин А.В. Численные методы. - М.: Наука, 1989. - 432 с.
Поступила 02.06.2009 г.
