Моделирование сосуществования конкурирующих систем хишник–жертва в пространственно неоднородной области Текст научной статьи по специальности «Математика»

УДК 519.63

МОДЕЛИРОВАНИЕ СОСУЩЕСТВОВАНИЯ КОНКУРИРУЮЩИХ СИСТЕМ ХИЩНИК-ЖЕРТВА В ПРОСТРАНСТВЕННО НЕОДНОРОДНОЙ ОБЛАСТИ

© 2013 г. А.В. Будянский

Будянский Александр Владимирович - аспирант, Южный федеральный университет, ул. Мильчакова, 8 а, г. Ростов-на-Дону, 344090, e-mail: [email protected].

Budyansky Alexander Vladimirovich - Post-Graduate Student, Southern Federal University, Milchakov St., 8 a, Rostov-on-Don, 344090, e-mail: [email protected].

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

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

We consider a model of two populations ofpreys that are closely related as well two populations ofpredators. We employ the method of lines and study the impact of migration on scenarios of local competition and coexistence of species. Conditions on system parameters are determined for appearance of nontrivial family of steady states.

Keywords: population dynamics, nonlinear parabolic equations, cosymmetry, method of lines.

Изменение климата, урбанизация и развитие сельского хозяйства приводят к миграции животных, что в свою очередь порождает конкурентную борьбу близкородственных видов. Это делает актуальным моделирование сосуществования видов в одной экологической нише [1, 2]. Принцип Гаузе [3] утверждает, что в случае конкуренции за единый ресурс происходит вытеснение менее приспособленной популяции. Имеются примеры [4] того, что при определенных условиях близкие виды способны сосуществовать. В [5] дано обобщение принципа Гаузе и показано, что при нелинейности миграционных потоков конкуренция биологических видов может не приводить к вытеснению менее приспособленной популяции. Конкуренция за обобщенный ресурс близкородственных популяций была рассмотрена в [6]. При моделировании эволюции сложных экосистем необходимо учитывать трофистические взаимодействия видов, в частности, отношения между популяциями хищников и жертв [1, 2, 7-9].

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

Модель динамики конкурирующих популяций

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

ды

дЫг = -V-q, + Л, i = 1 *4, V = (dx,dyf

dt

(1)

с миграционными потоками Ч и функциями роста (убыли) £ для популяций жертв и хищников

Чг = к + Щ^Р + Эг-,3-гмг^м3-г , 4 12

£г = Щгиг£0 - иг 2 ¡уи у , /о = 1--2 и,, г = 1,2,

У =3 р] =1

2

=-к¿Уиг + иг 2 Ру Vиу, (2)

у=1

2

/г = иг 2 Щуиу - ¡¡щ, г = 3,4. (3)

j=1

Для описания динамики близкородственных популяций аналогично [6] используется единая функция обобщенного ресурса (емкости среды) р(х, у). Изменение плотности популяций жертв определяется логистическим законом и убылью из-за присутствия хищника (слагаемые с коэффициентами ¡у) [1]; различия в естественном приросте - параметрами роста щ; прирост плотности популяции хищника за счет потребления им популяций жертв - слагаемыми с коэффициентами Щу , а естественная убыль - слагаемыми ¡^ [2].

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

порядка). Неоднородность обобщенного ресурса по ареалу проявляется градиентами Vp в выражениях

миграционных потоков жертв, а также поточечно через функции /. Модель также учитывает кросс-диффузию, влияние неравномерности распределения одного вида на миграцию другого [11]. Так, перенос плотности популяции и^ (х, у, t), вызванный неравномерностью распределения плотности популяции иу(х,у,t), описывается слагаемым ^уи^иу .

Рассматривался ареал в виде прямоугольника ^ = [0, а\х[0, Ъ\, на границах которого ставятся условия отсутствия потоков: 4л(0, УД ) = Чл{а, У, t)= 412 (x,0,t )= Ча (ХЪ,1: )= 0 > i = 1 — 4 . (4)

Система дополняется начальными распределениями для плотностей популяций

ui (х, y,ö) = uö (x, y,ö), i = 1 ^ 4 .

(5)

Рассматриваемая задача относится к классу ко-симметричных систем [8]. В этом случае возможно возникновение непрерывных семейств стационарных распределений популяций [6, 12].

Анализ показывает, что система (1)-(5) обладает косимметрией вида

L = ("Ль Пэ. ^4,), П2 = q2u1>

Пэ = П4 = ?4u3 (6)

при дополнительных условиях на вещественные параметры q j, которые получаются из условия ортогональности вектора L правой части системы (1)-(4).

2 4

J Z[(-V-qt + fН]+S[(-V-дг + fК-ddy = 0. (7) Qi=1 i=3

Равенство (7) выполняется для любых функций ui (x, y, t) при дополнительных связях между вещественными параметрами qj и коэффициентами системы

k1?2 = -k2?1 > М-1?2 = > l1 jq2 = -l2jSb j = 3,4 > (8)

k3q4 = -k4q3, Цэj?4 = ?3 > j =1,2 >

l3q4 = -l4q3 .

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

q1 = -k2,11, q3 = -k4,11 .

Численный метод

Для численного решения задачи (1)-(5) применяется метод прямых с дискретизацией на основе смещенных сеток. По переменным x и y вводятся равномерные сетки: xr = rhx, ys = shy, r = -1,0,К , nx +1, s = -1,0,К , ny +1, hx = a / nx, hy = b / ny; uirs - значение плотности распределения популяции ui в узле (xr, ys). Для вычисления потоков вводятся вспомогательные сетки: xr+^2 = -hx /2 + rhx, r = 1,К ,nx;

Xs+1/2 = ~hy /2 + shy, s = 1,K , Пу. Компоненты

ков qii и qi2 определяются в узлах | ys I и

Xr, У,+ 1/

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

i 1 \ +1,s ror,s i , \ (d1®)r+12, s =-1-, (d2ro);

2Ю)г, z+12

®r,s+1 ®r, s

h

У

и операторы вычисления среднего

®r+1, s +®r,s \ ®r,s+1 + ®r, s

M,

V+12,« =-~-, (82ю) г, г+12 ^

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

= [^Чл + ¿2^2 + / \га > г = 0,К , пх > « = 0,К , пу >

^ )г+12,« =

= [{к }циг " К К^ЬР5^« - {Рг,3-г }11^и3-г5иг \

^ )г,«+12

r+12, s

[{ki J22d2ui "К)22d2P^2ui "{ßi,3-i }29du3-i8U

I

i = 1,2 (q«1)

+ 12, s

{ki}nd1ul" ^^ {ßij j11 dui8u

(q,2),

2 )r, s+1/2

J=1 2

-+12,

{ki }22d1ui " T {ßj }22dui^uj

J =1

(9)

r, s+1/2

i = 3,4.

Дискретные аналоги краевых условий записываются с помощью законтурных узлов

(Яи )-1 / 2,« =(1п )1/2,« , (Яг1 Хх -1/2,« =(1п к +1/2,« ,

(10)

s = 0,K , n

У '

пото-

(Ч12 ) г,-1/2 (4 12 ) г,1/2 ' (4г2)г,Иу-1/2 =(4г2)г,Иу +1/2 ' Г = 0,К ,Пх , ' = 1 -4 .

Из (5) получаются начальные условия для (9), (10): Щ,Г8 = и° (ХГ, У«) . ' = 1 - 4 . Г = 0,К , Пх >

« = 0,К , Пу.

Система (8), (9) записывается в виде

#= Е (г), У (0) = У0. (11)

Здесь вектор У = (и1,и2,и3,и4содержит значения переменных в узлах сетки

и1 = Цдь и/,12,К , uг,1пy ,uг,21, и;,22,К , и1,пхпу

/ = 1 — 4 . У0 - вектор начальных данных. Для интегрирования по времени задачи Коши (11) применяется метод Рунге-Кутта четвертого порядка.

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

Косимметрией системы (11) является дискретный аналог векторного поля Ь, получаемый из (6) в результате дискретизации. Поскольку нулевые решения задачи (11) (вектор У = 0) аннулируют косимметрию, любое

ненулевое стационарное решение У (т.е. Р(У*) = 0) не обнуляет косимметрию и, таким образом, принадлежит однопараметрическому семейству равновесий. В спектре устойчивости равновесия У* будет нулевое значение, которое соответствует нейтральному направлению вдоль семейства. Если остальные спектральные величины лежат в левой полуплоскости, то равновесие У* устойчиво. Это отвечает устойчивости в трансверсальном к семейству многообразии [10].

Численные результаты

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

. Функция обобщенного

хищников

*3, v2

k з =

ресурса р(х, у) соответствует наличию на ареале двух благоприятных зон (рис. 1).

Начальные распределения и0,размещались в локализованных зонах (рис. 1б, в). Расчеты проводились до выхода на устойчивые стационарные распределения. В вычислениях параметры диффузии и роста-убыли фиксированы:

к! ^ 0 ], к 2 =Г0Д)4 0 ] ,

1 1 0 0,03/ 2 1 0 0,04) 0,06 0 ^ (0,08 0 0 0,06/ 4 = ^ 0 0,08)

Щ3 = 3 Щ4 = 4 , Щ31 = 0,9, Щ41 = 1,2, Щ32 = 1,2, Щ42 = 1,6 ; ¡13 = ¡14 = 0,9, ¡23 = ¡24 = 1,2 , ¡3 = 0,9, ¡4 = 1,2.

Данный набор параметров удовлетворяет условию косимметрии (8). В ходе эксперимента менялись миграционные коэффициенты.

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

стационарное сосуществование популяций хищников и жертв. При этом имеется преобладание популяции и2 над и1 на ареале (рис. 2а, б). Это обусловлено более высоким коэффициентом роста и близким расположением к благоприятным для роста видов зонам в начальный момент времени. Аналогичная ситуация складывается с популяциями хищников (рис. 2в, г). При этом финальные профили распределений популяций жертв по ареалу по форме напоминают функцию обобщенного ресурса р(х, у) . Распределения популяций хищников имеют экстремум в наиболее благоприятной зоне.

Для полученных финальных распределений рассчитан спектр устойчивости. Наличие практически

нулевого значения в спектре -1,17 -10-6 говорит о том, что данное решение принадлежит непрерывному семейству стационарных режимов и, стартуя из иных начальных условий, можно получить другие финальные распределения, входящие в непрерывное семейство равновесий [10].

Сосуществование конкурирующих видов имеет место и при отсутствии косимметрии, например, при учете неравномерности распределения ресурса и соседних популяций. Наблюдается фиксация ненулевых уровней близкородственных популяций. На рис. 3 показаны финальные распределения для следующего набора миграционных параметров: (0,05 0 ^ (0,02 0

а1 =1 „ а2 =1

0 0,05

0 0,02

Pl2 =

- 0,03 0

0

- 0,03

P21 =

- 0,03 0

0

- 0,03

p31 =р32 =р41 :

Рис. 1. Функция обобщенного ресурса (а), начальные распределения хищников (в) и жертв (б)

(0,02 0 ^ _(0,05 0 ^ ^ 0 0,02) Р42 = [ 0 0,05) . Сравнение рис. 2 и 3 показывает, что для выбранных параметров неравномерность р(х, у) сказывается

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

в

в г

Рис. 2. Финальные распределения при отсутствии направленной миграции

Рис. 3. Финальные распределения при наличии направленной миграции

В результате этого популяция и2 (х, у, /) получает возможность роста вне максимумов и^х, у, /). Такое размежевание происходит из-за отрицательности диагональных элементов матриц Р12 и р21. Профили плотностей популяций хищников также характеризуются большей концентрацией в наиболее благоприятной для роста зоне.

Выводы

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

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

Автор выражает благодарность научному руководителю В.Г. Цибулину за внимание к работе и Ю.В. Тютюнову за полезные замечания.

Исследование проводилось при финансовой поддержке гранта РФФИ № 11-01-02.

Литература

1. Мюррей Дж. Математическая биология. Пространственные модели и их приложения в биомедицине. Т. 2. М., 2011. 1104 с.

2. Базыкин А.Д. Математическая биофизика взаимодействия популяций. М., 1985. 182 с.

б

a

в

г

3. Гаузе Г.Ф. Борьба за существование. М., 2002. 175 с.

4. Бигон М., Харпер Дж., Таунсенд К. Экология. Особи, популяции и сообщества. М., 1989. 477 с.

5. Белотелое Н.В., Лобанов А.И. Популяционные модели с нелинейной диффузией // Мат. моделирование. 1997. Т. 9, № 12. C. 43-56.

6. Будянский А.В., Цибулин В.Г. Моделирование пространственно-временной миграции близкородственных популяций // Комплексное исследование и моделирование. 2011. Т. 3, № 4. С. 477-488.

7. Говорухин В.Н., Моргулис А.Б., Тютюнов Ю.В. Медленный таксис в модели хищник-жертва // Докл. РАН. 2000. T. 372, № 6. C. 730-732.

8. Tyutyunov Yu., Titova L., Arditi R. A minimal model of pursuit-evasion in a predator-prey system // Mathematical Modelling of Natural Phenomena. 2007. Vol. 2, № 4. P. 122-134.

9. Тютюнов Ю.В., Жадановская Е.А., Ардити Р., Мед-винский А.Б. Пространственная модель развития устойчивости насекомых-вредителей к трансгенной инсектицидной сельскохозяйственной культуре на примере кукурузного стеблевого мотылька // Биофизика. 2007. Т. 52, вып. 1. С. 95-113.

10. Юдович В.И. Косимметрия, вырождение решений операторных уравнений, возникновение фильтрационной конвекции // Мат. заметки. 1991. Т. 49, № 5. С. 142-148.

11. Цыганов М.А., Бикташев В.Н., Бриндли Дж., Хол-ден А.В., Иваницкий Г.П. Волны в кросс-диффузионных системах - особый класс нелинейных волн // Успехи физ. наук. 2007. Т. 177, № 3. С. 275-300.

12. Ковалева Е.С., Цибулин В.Г., Фришмут К. Семейство стационарных режимов в модели динамики популяций // Сиб. журн. индустриальной математики. 2009. Т. 12, № 1. C. 98-107.

Поступила в редакцию_10 июля 2012 г.