CC BY
26
Известия высших учебных заведений. Прикладная нелинейная динамика. 2021. Т. 29, № 5 Izvestiya Vysshikh Uchebnykh Zavedeniy. Applied Nonlinear Dynamics. 2021;29(5)
Научная статья УДК 530.182
DOI: 10.18500/0869-6632-2021-29-5-751-764
Нелинейная динамика системы хищник-жертва на неоднородном ареале и сценарии локального взаимодействия видов
В.Г.Цибулин1, Т.Д. Ха1'2, П. А. Зеленчук1М
1Южный федеральный университет, Ростов-на-Дону, Россия 2 Вьетнамско-Венгерский индустриальный университет, Ханой, Вьетнам E-mail: [email protected], [email protected], [email protected] Поступила в редакцию 12.01.2021, принята к публикации 15.04.2021, опубликована 30.09.2021
Аннотация. Цель настоящей работы - изучить влияние различных локальных моделей в уравнениях диффузии-адвекции-реакции на пространственные процессы сосуществования хищников и жертв в условиях неоднородного распределения ресурса жертвы. Рассматривается система нелинейных уравнений параболического типа, учитывающая диффузию, таксис и локальное взаимодействие хищника и жертвы на одномерном ареале. Методы. Исследование системы проводится с помощью анализа динамических систем на фазовой плоскости, а также вычислительного эксперимента на основе метода прямых и схемы смещённых сеток. Результаты. Изучено поведение системы хищник-жертва при различных вариантах описания локального взаимодействия, учитывающих гиперболический закон роста жертвы и эффект Холлинга II рода при неравномерности распределения пищевого ресурса для жертвы. Установлены парадоксальные сценарии взаимодействия жертвы и хищника для ряда вариантов трофической функции. Проанализировано формирование стационарных и нестационарных решений при учёте диффузии и направленной миграции видов. Заключение. Предложена учитывающая неоднородность ресурса трофическая функция, которая не приводит к парадоксальной динамике.
Ключевые слова: модель хищник-жертва, нелинейная динамика, неоднородный ареал, диффузия, таксис, трофическая функция.
Благодарности. Работа выполнена при поддержке гранта Правительства РФ № 075-15-2019-1928.
Для цитирования: Цибулин В. Г., Ха Т.Д., Зеленчук П. А. Нелинейная динамика системы хищник-жертва на неоднородном ареале и сценарии локального взаимодействия видов//Известия вузов. ПНД. 2021. T. 29, № 5. С. 751-764. DOI: 10.18500/0869-6632-2021-29-5-751-764
Статья опубликована на условиях Creative Commons Attribution License (CC-BY 4.0).
©Цибулин В. Г., Ха Т.Д., Зеленчук П. А., 2021
751
Article
DOI: 10.18500/0869-6632-2021-29-5-751-764
Nonlinear dynamics of the predator-prey system in a heterogeneous habitat and scenarios of local interaction of species
V.G. Tsybulin1, T.D. Ha1'2, P. A. Zelenchuk1M
1 Southern Federal University, Rostov-on-Don, Russia 2Vietnam-Hungary Industrial University, Hanoi, Vietnam E-mail: [email protected], [email protected], [email protected] Received 12.01.2021, accepted 15.04.2021, published 30.09.2021
Abstract. The purpose of this work is to study the influence of various local models in the equations of diffusion-advection-reaction on the spatial processes of coexistence of predators and prey under conditions of a nonuniform distribution of the carrying capacity. We consider a system of nonlinear parabolic equations to describe diffusion, taxis, and local interaction of a predator and prey in a one-dimensional habitat. Methods. We carried out the study of the system using the dynamical systems approach and a computational experiment based on the method of lines and a scheme of staggered grids. Results. The behavior of the predator-prey system has been studied for various scenarios of local interaction, taking into account the hyperbolic law of prey growth and the Holling effect with nonuniform carrying capacity. We have established paradoxical scenarios of interaction between prey and predator for several modifications of the trophic function. Stationary and nonstationary solutions are analyzed considering diffusion and directed migration of species. Conclusion. The trophic function that considers the heterogeneity of the resource is proposed, which does not lead to paradoxical dynamics.
Keywords: predator-prey system, nonlinear dynamics, heterogeneous habitat, diffusion, taxis, trophic function.
Acknowledgements. This work was supported by a grant from the Goverment of the Russian Federation No. 075-15-2019-1928.
For citation: Tsybulin VG, Ha TD, Zelenchuk PA. Nonlinear dynamics of the predator-prey system in a heterogeneous habitat and scenarios of local interaction of species. Izvestiya VUZ. Applied Nonlinear Dynamics. 2021;29(5):751-764. DOI: 10.18500/0869-6632-2021-29-5-751-764
This is an open access article distributed under the terms of Creative Commons Attribution License (CC-BY 4.0).
Введение
Моделирование систем типа хищник-жертва или хозяин - паразит является важным для математической биологии, современной экологии, и актуально при решении проблем медицины, пищевой промышленности, борьбы с вредителями и пр. [1,2]. Первые математические модели базировались на обыкновенных дифференциальных уравнениях, описывающих локальное взаимодействие - усреднённую реакцию конкурирующих или антагонистических видов друг на друга. Для учёта пространственного распределения видов далее стали применяться уравнения диффузии, благодаря которым появилась возможность моделировать миграционные процессы [3,4]. К настоящему времени предложено большое число систем, относящихся к классу уравнений «диффузия-адвекция-реакция», в которых рассматриваются различные эффекты популяционно-го взаимодействия [5,6].
Актуальным является построение моделей, учитывающих неоднородность среды обитания. В условиях сильной диффузии и/или адвекции влияние пространственной неоднородности ареала мало заметно и в большинстве задач не учитывается [7]. В случае миграционных потоков малой интенсивности важен корректный учёт членов реакции [8]. Неоднородность ресурса, как правило, учитывается только при описании локального роста жертвы и не входит в трофическую функцию (функциональный отклик хищника), см. [9-13].
Это, как будет показано далее, приводит к парадоксальным ситуациям. Так, например, на некоторых участках ареала жертва оказывается «индифферентной» к распределению своего ре-
сурса, а хищник может отсутствовать в местах, где жертвы «достаточно» много. Между тем, данные полевых наблюдений показывают [1], что между видами всегда имеется корреляция в распределении вдоль ареала даже при малых потоках. Действительно, хищнику не обязательно совершать много хаотических движений в поисках жертвы (сильная диффузия) или преследовать её по запаху через весь ареал (быстрый таксис), достаточно затаиться и ждать в местах с большим количеством ресурса жертвы, куда она неизбежно придёт. Конечно, эта стратегия не является единственной, но учёт тенденции к такого рода поведению хищника, на наш взгляд, является вполне разумным.
В настоящей работе рассматриваются различные варианты описания локального взаимодействия с использованием трофической функции Холлинга II рода [14] и учётом неоднородности ресурса жертвы. На основе бифуркационного анализа локальной модели (миграция и диффузия отсутствуют) и численного анализа полной модели установлен вид трофической функции, не приводящий к парадоксальным ситуациям при малых миграционных потоках.
1. Уравнения динамики хищника и жертвы на неоднородном ареале
Математическая модель пространственно-временного взаимодействия жертвы с плотностью и(х, ¿) и хищника с плотностью ь(х, может быть записана в виде системы уравнений [1,4]
дИ . гр ■ дЧ2 . „ пл
и =-----+ «^1, У =-----+ УГ2, (1)
ОХ ОХ
где точка обозначает дифференцирование по времени ¿, а функции ^ и ^ описывают соответственно миграционные потоки и локальное взаимодействие (реакцию) видов. Потоки ^ задаются следующим образом [5,6]:
ди <9ф1 дь <9ф2
Я1 = ~к1— + и^—, Ц2 = + V——. (2)
ох ох ох ох
Здесь первое слагаемое характеризует диффузию, а второе слагаемое отвечает за направленную миграцию - таксис, определяемый с помощью функций ф^, которые могут быть выражены в виде [15,16]
ф1 = ар - впи - 012V, ф2 = 02Ш - $22У. (3)
Функция ф1 состоит из трёх частей, которые определяют различные виды направленной миграции: таксис жертвы на ресурс р = р(х) и от мест с избыточным скоплением особей своего вида (-011-и), а также от хищника (-012^). Функция ф2 состоит из двух частей: 021^ - это таксис хищника, направленный на жертву, (-022^) - таксис, направленный от мест, где охотятся другие хищники.
В уравнениях (2) и (3) коэффициенты кг, а, 0^ (г,] = 1, 2) являются постоянными положительными величинами, значения которых определяются из данных наблюдения. Функция р(х) > 0 описывает неравномерное распределение ресурса жертвы вдоль ареала. Локальное взаимодействие видов выражается следующим образом [4,11]
« = <-> - ТТЬ « = -2 + (4)
Первое слагаемое в задаёт рост популяции жертвы, причём функция / (и) имеет вид [3]
(1 - ?) ■
/ (и) = (Ьа + 61^(1 - ) , 0 < 8о, 81 < 1 (5)
и характеризует закон роста (80, 81 - константы). Так, при 80 = 1, 81 =0 будем иметь логистический, а при 80 = 0, 81 = 1 - гиперболический закон роста, возможны также и их комбинированные варианты [10]. Первое слагаемое в функции отвечает за естественную убыль хищника. Оба вторых слагаемых в (4) описывают эффект хищничества с помощью трофической функции Холлинга II рода [11,14]. Положительные коэффициенты 61 и 62 характеризуют соответственно убыль жертвы и прирост хищника в результате их контакта. Неотрицательный параметр С позволяет учесть инертность хищника, проявляемую им при поиске, поглощении и переработке жертвы [11]. В случае С = 0 трофическая функция приобретает классический вид Лотки-Вольтерры.
В уравнениях (1)-(4) все коэффициенты могут быть функциями от ж и ¿, но в данной работе предполагается только пространственная зависимость коэффициентов трофической функции хищника 62 и С, причём эта зависимость соотносится с функцией ресурса жертвы.
Число параметров в системе (1)-(4) можно уменьшить, проведя замену переменных £ , V и вводя новые коэффициенты (предполагается независимость параметров а1, Ъ1 от х)
М , «2
а^ ^ I,--> V, К = —,
а1 а1
в = —.
а1
(6)
Коэффициент К далее используется в качестве бифуркационного параметра, который характеризует отношение коэффициентов убыли хищника и роста жертвы. В результате получается
система
дя 1 , и = —-—+ и дх
/ (и) —
1 + Си
V = —
д 92
д х
+ У
—К +
в и
1 + Си
Рассматривается кольцевой ареал, система (1)-(3) дополняется условиями периодичности:
(7)
и (о, г) = и (1, г), ql (о, г) = ql (1, г), V (0,1)=у (1,1), <?2 (о, *) = ^ (1, ¿),
(8)
и начальными распределениями плотностей популяций:
и ( х, 0) = и0 ( х) , ( х, 0) = 0 ( х)
(9)
2. Численный метод решения неоднородной начально-краевой задачи
Для численного решения задачи (1)-(5) применяется метод прямых с дискретизацией на основе смещенных сеток [12]. По переменной х вводится равномерная сетка:
хг = гН, г = 1,
м,
Н = N
(10)
Плотность распределения популяции щ в узле хг далее обозначается через . При вычислении
потоков используется вспомогательная сетка: хг, 1
' + 2
системы (1)-(5) по прост тегрируется по отрезку х 1 ,хг+1
= г Н + Щ ,г = 1,
, N. Для аппроксимации
эанственной координате применяется метод баланса: уравнение (1) ин-
а для потоков ^ интегрирование проводится по отрезку
[хг,хг+1]. При этом используются операторы
( йт)г =
1 — ыг-1
(8ы)г =
1
(11)
Н
2
и условия периодичности шм+1 = ■1,т 1 = +1. В результате получается следующая система ОДУ: 2 2
йг = [-йд 1 + /1 ]г , ¡1,г = иг[ Ц1/0,г -
/о,г = (бо + б^ Н1 - —
1 + Сг и.
г = 1, ...,Ы,
(1 - ю • * = * Г?
+2 йх
р(х)
1
Уг = \-dq2 + ¡2]г , ¡2,г = Уг[ -X +
-Ву '^'ТГ
1 + Сг и.
:) ■
= 1,
д1 г+1 = [-к1йи + айрби + 011йиби + 012йу би]г+1
Ч2,Г+1 = [-^2<1у + 021 йибу + 022йь бу ]г+ 1 . Построенная конечномерная модель с дискретными аналогами условий периодичности
им+1 = щ, д1,1 = д1М +1, ум+1 = У1, д2> 1 = д2^м +1,
может быть записана в виде
!¥ = Ф (Ш), Ш (0) = Ш0,
(12)
(13)
(14)
(15)
(16)
(17)
здесь Ш = (и1,и2,..., им, и1, и2,..., им) - вектор значений переменных в узлах сетки. Начальные данные Ш0 для системы (12)-(16) получаются из (4).
Для интегрирования системы (17) по времени использовался метод Рунге-Кутты четвертого порядка.
3. Локальное взаимодействие хищника и жертвы
Рассмотрим систему (7) при отсутствии диффузии и таксиса
и = и
!(и) -
1 + Си
V = V
-X +
Ви
1 + Си
(18)
Для любых значений параметров система (18) имеет неустойчивое нулевое равновесие и = V = 0 и равновесие без хищника
и = р, V = 0. (19)
Для равновесия (19) получается следующее характеристическое уравнение
2 (Ср2б1 + Србо + СрХ -Вр + б1р + бо + X) (СрХ -Вр + Х)(б1р + бо)
о2 +
Ср + 1
-о +
Ср + 1
0. (20)
С учётом положительности функции ресурса р(х) и коэффициентов бj, С, X получаем, что равновесие (19) устойчиво при X > Хсг, где
Хст -
В
Ср + 1
(21)
г
2
При X < Xcr возникает равновесие, отвечающее сосуществованию хищника и жертвы,
I С2 X2p + CX2 + В2р - 2 С В Xp-В \ В Xp-CX2p - X2 U* = В-CX, V* =-(В - CX)3p-В8° + (В - СX)3p BSl- (22)
С уменьшением X это равновесие теряет устойчивость в результате бифуркации Пуанкаре-Андронова-Хопфа, соответствующее критическое значение будем обозначать X*.
Из (21) видно, что при постоянных В и С критическое значение Xcr явно зависит от p(x). Логично предположить обратную пропорциональность коэффициентов В и С от p(x). Далее анализируются различные варианты зависимости параметров В и С от функции ресурса ( x), см. Таблицу (п, с - постоянные). В Таблице также приведены аналитические выражения для равновесия (u*, v*) и границ его устойчивости Xcr, X*.
В первых трёх случаях критические значения Xcr и X* явно зависят от функции ресурса. Только вариант IV приводит к независящим от x значениям Xcr, X*.
Для иллюстрации распределения критических значений Xcr, X* и плотностей видов u*, v* на отрезке [0,1] были проведены вычисления для функции ресурса с единичными максимумом и минимумом p(x) = 1 + I sin(2nx). На рис. 1 приведены графики Xcr, X*, прямые линии отвечают варианту IV, кривые соответствуют вариантам I, II и III из Таблицы.
Использование вариантов с Xcr, зависящим от p(x), может привести к тому, что для соответствующей модели на ареале могут получаться как равновесия с сосуществованием хищника и жертвы, так и решения без хищника. Это имеет место для случаев I, II и III. Например, для варианта III графики Xcr и X* находятся в противофазе распределению p(x). Это, в частности, влечёт сосуществование хищников и жертв при малой величине ресурса p(x) и отсутствие хищника при больших p(x). В этом случае возникает контраст в распределении жертвы до точки бифуркации X = Xcr и после, см. рис. 2, где приведены распределения жертвы и* и хищника v* для всех случаев при X = 1.1. Возникает парадоксальная ситуация - на части ареала распределение хищника «повторяет» функцию ресурса p(x), а жертва к нему «индифферентна» (рис. 2, a). Также имеются участки, где при сравнимых плотностях жертвы хищник может отсутствовать (рис. 2, b). Наконец, возможен вариант, когда хищник сосуществуют с жертвой при минимальных уровнях функции ресурса, а при максимальных - отсутствует (рис. 2, c).
В случае модели IV критическое значение Xcr не зависит от p(x), что приводит к ненулевым распределениям u*(x), v*(x), причём с уменьшением параметра X происходит снижение уровня плотности жертвы и рост плотности хищника.
Таблица. Стационарное решение сосуществующих видов и границы его устойчивости для различных вариантов определения параметров трофической функции
Table. The stationary solution of co-existing species and the limits of its stability for various
options for determining the parameters of trophic function
В С и* V* X* Xcr
I II III IV П с n -п — с p П с X П - сХ pX Xn(np - cpX - X) p(n - cX)3 Xp2n(pn - cX - X) 1 + у/c2p2 + cp + 1j n c(cp + 1) pn (-1 + Vc2 + c+1) np p + 1 np
pn - сХ pX (pn - cX)3 p^n(n - cpX - X) , c(c+1) N 1 + \Jc2p2 + cp + 1 j n c+1 n
п- p pX (n - pcX)3 Xpn(n - cX - X) cp(cp + 1) n (-1 + Vc2 + c+1) p + 1 n
п- (n - cX)3 c(c+ 1) c+1
1.6 1.6
1.3 1 1.3
1 1 х 7 \ \ х X 1
\2/
0.7 0.7
а
Рис. 1. Графики зависимости от х критических значений Xcr (штриховые кривые) и X* (сплошные кривые): a - сравнение вариантов I (кривые 1, 2) и IV (прямые 7, 8), b - II (кривые 3, 4) и IV, c - III (кривые 5, 6) и IV;
П = 4, с = 2.5, р(х) = 1 + 2 sin(2nx)
Fig. 1. Displaying the variation with x of critical values Xcr (dashed curves) and X* (solid curves): a - the comparison of variants I (curves 1, 2) and IV (direct lines 7,8), b - II (curves 3, 4) и IV, c - III (curves 5, 6) and IV; n = 4, с = 2.5,
p(x) = 1 + 2 sin(2rnx)
a
Рис. 2. Стационарные распределения хищника г>(ж) (сплошная) и жертвы и(х) (штриховая): a - вариант I, b - II, c - III, d - IV; при X = 1.1, n = 4, с = 2.5, функция ресурса р(х) = 1 + | sin(2rnx) (штрихпунктир)
Fig. 2. Stationary distributions of the predator г>(ж) (solid) and the prey u(x) (dashed): a - variant I, b - II, c - III, d - IV; for X = 1.1, n = 4, с = 2.5, the function of resource is p(x) = 1 + 5 sin(2rnx) (dashed dotted lines)
0
1
0
1
0
1
0
1
4. Пространственная динамика хищника и жертвы
В бездиффузионном приближении для 0 < X < X* равновесие (и*, V*) становится неустойчивым, и возникает предельный цикл. Для демонстрации пространственно-временной миграции при разных вариантах определения членов, описывающих эффект Холлинга, на рис. 3-8 представлены результаты вычисления динамики при начальных распределениях, отвечающих возмущению стационарного состояния популяции хищника: и,0(х) = и*(х), у0(х) = у*(х)+ е. На рис. 3 приведены распределения и(х, ¿), у(х, ¿) для варианта II при п = 4, с = 2.5, X = 1.32, е = 0.1. В бездиффузионном приближении (ку = а = = 0) на правой части ареала (х > 0.5) возмущения для хищника в основном затухают, стремясь к нулю. На большой части ареала устанавливается стационарное распределение плотностей жертвы и хищника, а посередине левой части (х « 0.25) формируется колебательный режим. При учете диффузии (к1 = 0.001, к2 = ^, а = Р12 = Р21 = 0) колебательный режим формируется на всем ареале, при этом плотность хищника в правой части ареала существенно меньше, чем в левой части, напоминая распределение в бездиффузионном приближении. Аналогичные результаты получаются для вариантов I и III.
С увеличением диффузионных коэффициентов происходит выравнивание уровней плотностей распределения хищников (рис. 4). При малых ^ сохраняется разница в уровнях плотности хищников для участков ареала, отвечающих максимуму и минимуму функции ресурса. С увеличением диффузионных параметров эта диспропорция сокращается.
При достаточно больших X происходит выход на стационарное распределение. Для модели IV при X = 1.05 получается стационарное решение, а колебательный режим формируется при X < 1. Для первых трех указанных в Таблице вариантов имеются диапазоны значений диффузионных параметров, при которых реализуются колебательные режимы. С увеличением параметра X колебания сохраняются при малых значениях диффузии, а для больших значений формируются стационарные распределения.
На рис. 5 приведены результаты вычисления таких состояний при X = 1.1 и к1 = 0.05, к2 = а = |321 = 0.025, |312 = 0. Видно, что получаются похожие распределения. Это означает, что для достаточно больших значений коэффициентов диффузии и адвекции уменьшаются
Рис. 3. Пространственно-временные распределения жертвы (a, b) и хищника (c, d) для варианта II; a, c - бездиффузионное приближение; b, d - при k1 = 0.001, к2 = 4г, а = pi2 = Р21 = 0; X = 1.32, n = 4, с = 2.5, р(х) = 1 + 5 sin(2nx), начальное распределение и0(ж) = и*, г>0(ж) = v* + 0.1
Fig. 3. Spatial-temporal distributions of the prey (a, b) and the predator (c, d) for variant II; a, c - the diffusionless approximation; b, d - for ki = 0.001, k2 = 4т, а = P12 = P21 =0; X = 1.32, n = 4, с = 2.5, p(x) = 1 + 5 sin(2rnx), the initial distribution is u0(x) = и*, г>0(ж) = v* + 0.1
Рис. 4. Зависимость от времени численности жертвы (a) и хищника (b) для варианта II; х = 0.75, X = 1.32,
П = 4, с = 2.5, р(х) = 1 + § sin(2rnx), к2 = a = ^, Р12 = 0, Р21 = 0
Fig. 4. Displaying the variation with time of the prey (a) and the predator (b) for variant II; x = 0.75, X = 1.32, n = 4, с = 2.5, p(x) = 1 + § sin(2rnx), k2 = a = , P12 = 0, P21 = 0
Рис. 5. Стационарные распределения хищника v(ж) (сплошная кривая) и жертвы и(х) (штриховая) при «большой диффузии»: a - вариант I, b - II, c - III, d - IV; X = 1.1, n = 4, с = 2.5, p(x) = 1 + 2 sin(2rnx) (штрихпунктир); k1 = 0.05, fc2 = 0.025, a = 0.025, P12 = 0, P21 = 0.025
Fig. 5. Stationary distributions of the predator г>(ж) (solid curve) and the prey u(x) (dashed) for a «large diffusion»: a - variant I, b - II, c - III, d - IV; X = 1.1, n = 4, с = 2.5, p(x) = 1 + 2 sin(2rnx) (dashed dotted line); k1 = 0.05, k2 = 0.025, a = 0.025, P12 = 0, P21 = 0.025
различия между всеми вариантами определения параметров В и С из Таблицы (исчезает «парадоксальность», отмеченная при малых миграционных параметрах).
Интересно сравнить результаты вычислений для большой диффузии, но без учета направленной миграции (а = |312 = |321 = 0) (рис. 6). Видно, что для вариантов I и II распределения хищника и жертвы сонаправлены - экстремумы получаются при близких значениях х, а для вариантов III и IV имеет место разнонаправленность, когда максимуму жертвы и ресурса отвечает минимум в плотности распределения хищника. Вариант I характеризуется наименьшей плотностью хищника, а вариант III - жертвы.
На рис. 7 приведены результаты, демонстрирующие влияние коэффициентов направленной миграции а, Р12, Р21 на формирование стационарных распределений при использовании варианта IV. Увеличение коэффициента а вызывает более плотное заполнение благоприятных по ресурсу частей ареала, выбор знаков и величин коэффициентов позволяет моделировать ситуации от осторожного до безразличного поведения вида по отношению к антагонисту. Разницу в пространственно-временной эволюции видов при «средней диффузии» (к1 = 0.02) демонстри-
Рис. 6. Стационарные распределения хищника г>(ж) (сплошная кривая) и жертвы и(х) (штриховая кривая) при «большой диффузии» и отсутствии направленной миграции: a - вариант I, b - II, c - III, d - IV; X = 1.1, n = 4, с = 2.5, функция ресурса р(х) = 1 + § sin(2rnx) (штрихпунктир), к1 = 0.05, к2 = 0.025, a = р12 = р21 = 0
Fig. 6. Stationary distributions of the predator г>(ж) (solid curve) and the prey u(x) (dashed) for a «large diffusion» and in the absence of directed migration: a - variant I, b - II, c - III, d - IV; X = 1.1, n = 4, с = 2.5, the function of resource is p(x) = 1 + § sin(2rnx) (dashed dotted line), fc1 = 0.05, fc2 = 0.025, a = p12 = p21 = 0
a
Рис. 7. Стационарные распределения хищника г>(ж) (сплошная кривая) и жертвы и(х) (штриховая) для варианта IV: a - Р12 = Р21 =0, b - Р12 = 0, Р21 = 0.025, c - Р12 = 0.025, Р21 = 0, d - P12 = 0.025, P21 = 0.025; функция ресурса р(х) = 1 + 2 sin(2rnx) (штрихпунктир), X = 1.1, n = 4, с = 2.5, к1 = 0.05, к2 = 0.025, a = 0.025
Fig. 7. Stationary distributions of the predator г>(ж) (solid curve) and the prey u(x) (dashed) for variant IV: a - P12 = P21 = 0, b - P12 = 0, P21 = 0.025, c - P12 = 0.025, P21 =0, d - P12 = 0.025, P21 = 0.025; the function of resource is p(x) = 1 + § sin(2rnx) (dashed dotted line), X = 1.1, n = 4, с = 2.5, fc1 = 0.05, fc2 = 0.025, a = 0.025
0
0
0
0
1
0
1
0
1
0
1
Рис. 8. Пространственно-временные распределения жертвы (a, b) и хищника (c, d) для вариантов II (a, c) и IV
(b, d); fcx = 0.02, = a = P21 = ^, P12 = \г; X = 0.95, n = 4, с = 2.5, p(x) = 1 + § sin2n®, начальное распределение u0(x) = 0.8, и0(ж) =0.1
Fig. 8. Spatial-temporal distributions of the prey (a, b) and the predator (c, d) for variants II (a, c) and IV (b, d); fc1 = 0.02, fc2 = a = p21 = ^§L, p12 = ^§L; X = 0.95, n = 4, с = 2.5, p(x) = 1 + § sin2rn®, the initial distribution is uo(x) = 0.8, -Уо(ж) = 0.1
рует рис. 8, где для Х=0.95 слева приведены результаты расчёта на основе варианта II, а справа -IV. Начальное распределение хищника и жертвы было одинаково, а пространственно-временные портреты заметно разнятся. Для обоих вариантов формируются колебательные режимы, но периоды и характеристики распределения по х отличаются.
Заключение
В работах по моделированию пространственно-временного взаимодействия сосуществующих антагонистических видов до настоящего времени мало учитывается многофакторность и неоднородность протекающих процессов. Это связано с трудностями получения надёжных данных о независимом развитии видов, сложности определения миграционных характеристик [17]. При моделировании приходится опираться на полуэмпирические оценки неравномерности распределения ресурсов и видов, делать предположения о характере взаимовлияния популяций и величинах параметров, см. например [7,18-21].
В данной работе предпринята попытка математического исследования системы (1)-(9), относящейся к классу уравнений «диффузия - адвекция - реакция» [5,6]. Для характеристики локальной кинетики используется закон гиперболического роста жертвы и трофическая функция с эффектом Холлинга второго рода. Анализируется область параметров системы, где существенны локальное описание динамики и учёт пространственных эффектов. В качестве бифуркационного выбирается параметр X - отношение коэффициентов естественной убыли хищника и роста жертвы. Проводится анализ устойчивости равновесий, определяются критические значения параметра X.
Установлено, что в бездиффузионном приближении (при отсутствии потоков) модель с постоянными, независящими от х коэффициентами в уравнении (7) (вариант I Таблицы), демонстрирует «парадоксальное» распределение видов (см. рис. 2, a), когда жертва может быть индифферентна к ресурсу при наличии хищника. Рассмотрены три варианта модификации трофической функции с зависимостью её коэффициентов от ресурса (варианты II, III и IV из Таблицы). Показано, что только вариант IV даёт адекватную картину распределения, когда плотности видов коррелируют друг с другом и с функцией ресурса, при этом достичь максимума ресурса жертве мешает именно наличие хищника. Это разница в результатах сохраняется и при малой диффузии.
С увеличением диффузии (и адвекции) результаты на основе моделей I (с постоянными коэффициентами В и С) и IV (с обратной зависимостью от функции ресурса) сближаются, таким образом, вариант IV позволяет описывать поведение системы хищник-жертва как при больших, так и при малых значениях диффузионных и миграционных коэффициентов.
Исследование влияния неоднородности при учёте диффузии и адвекции (таксиса) проводилось с помощью вычислительного эксперимента. Полученные результаты показывают важность предложенной модификации трофической функции для описания пространственно-временных популяционных сценариев.
Список литературы
1. Мюррей Дж. Математическая биология. Т. 2. Пространственные модели и их приложения в биомедицине. М.-Ижевск: НИЦ «Регулярная и хаотическая динамика», Институт компьютерных исследований, 2011. 1104 с.
2. Turchin P.B. Complex Population Dynamics: A Theoretical/Empirical Synthesis (MPB-35). Princeton: Princeton University Press, 2003. 472 p.
3. Базыкин А. Д. Нелинейная динамика взаимодействующих популяций. М.-Ижевск: Институт компьютерных исследований, 2003. 368 c.
4. Rubin A., Riznichenko G. Mathematical Biophysics. New York: Springer, 2014. 273 p. DOI: 10.1007/978-1-4614-8702-9.
5. Cantrell R.S., Cosner C. Spatial Ecology Via Reaction-Diffusion Equations. Chichester: John Wiley and Sons Ltd, 2003. 428 p. DOI: 10.1002/0470871296.
6. Malchow H., Petrovskii S. V., Venturino E. Spatiotemporal Patterns in Ecology and Epidemiology: Theory, Models, and Simulation. New York: Chapman and Hall/CRC, 2008. 469 p.
7. Фрисман Е.Я., Кулаков М.П., Ревуцкая О. Л., Жданова О. Л., Неверова Г. П. Основные направления и обзор современного состояния исследований динамики структурированных и взаимодействующих популяций // Компьютерные исследования и моделирование. 2019. Т. 11, № 1. С. 119-151. DOI: 10.20537/2076-7633-2019-11-1-119-151.
8. Kim K., Choi W. Local dynamics and coexistence of predator-prey model with directional dispersal of predator // Math. Biosci. Eng. 2020. Vol. 17, no. 6. P. 6737-6755.
DOI: 10.3934/mbe.2020351.
9. Courchamp F., Berec J., Gascoigne J.Allee Effects in Ecology and Conservation. Oxford: Oxford University Press, 2008. 256 p. DOI: 10.1093/acprof:oso/9780198570301.001.0001.
10. Свирежев Ю. М.Нелинейные волны, диссипативные структуры и катастрофы в экологии. М.: Наука, 1987. 368 c.
11. Тютюнов Ю.В., Титова Л. И. От Лотки-Вольтерра к Ардити-Гинзбургу: 90 лет эволюции трофических функций // Журнал общей биологии. 2018. Т. 79, № 6. С. 428-448. DOI: 10.1134/S004445961806009X.
12. Епифанов А. В., Цибулин В. Г. Моделирование колебательных сценариев сосуществования конкурирующих популяций // Биофизика. 2016. Т. 61, № 4. С. 823—832.
DOI: 10.1134/S0006350916040072.
13. Епифанов А. В., Цибулин В. Г. О динамике косимметричных систем хищников и жертв // Компьютерные исследования и моделирование. 2017. Т. 9, № 5. С. 799-813.
DOI: 10.20537/2076-7633-2017-9-5-799-813.
14. Holling C. S. Some characteristics of simple types of predation and parasitism // The Canadian Entomologist. 1959. Vol. 91, no. 7. P. 385-398. DOI: 10.4039/Ent91385-7.
15. Budyansky A. V., Frischmuth K., Tsybulin V. G. Cosymmetry approach and mathematical modeling of species coexistence in a heterogeneous habitat // Discrete & Continuous Dynamical Systems -B. 2019. Vol. 24, no. 2. P. 547-561. DOI: 10.3934/dcdsb.2018196.
16. Будянский А. В., Цибулин В. Г. Моделирование многофакторного таксиса в системе «хищник-жертва» // Биофизика. 2019. Т. 64, № 2. С. 343-349. DOI: 10.1134/S0006302919020133.
17. Abrams P. A. The evolution of predator-prey interactions: Theory and evidence // Annu. Rev. Ecol. Syst. 2000. Vol. 31, no. 1. P. 79-105. DOI: 10.1146/annurev.ecolsys.31.1.79.
18. Загребнева А. Д., Говорухин В.Н., Сурков Ф.А. Бифуркации в модели активный хищник -пассивная жертва // Известия вузов. ПНД. 2014. Т. 22, № 3. С. 94-106.
DOI: 10.18500/0869-6632-2014-22-3-94-106.
19. Tyutyunov Y. V., Titova L.I., Senina I. N.Prey-taxis destabilizes homogeneous stationary state in spatial Gause-Kolmogorov-type model for predator-prey system // Ecological Complexity. 2017. Vol. 31. P. 170-180. DOI: 10.1016/j.ecocom.2017.07.001.
20. Mishra P., WrzosekD. The role of indirect prey-taxis and interference among predators in pattern formation // Math. Methods Appl. Sci. 2020. Vol. 43, no. 18. P. 10441-10461.
DOI: 10.1002/mma.6426.
21. Haskell E. C., Bell J. Pattern formation in a predator-mediated coexistence model with prey-taxis // Discrete & Continuous Dynamical Systems - B. 2020. Vol. 25, no. 8. P. 2895-2921. DOI: 10.3934/dcdsb.2020045.
References
1. Murray JD. Mathematical Biology II: Spatial Models and Biomedical Applications. New York: Springer-Verlag; 2003. 814 p. DOI: 10.1007/b98869.
2. Turchin PB. Complex Population Dynamics: A Theoretical/Empirical Synthesis (MPB-35). Princeton: Princeton University Press; 2003. 472 p.
3. Bazykin AD. Nonlinear Dynamics of Interacting Populations. Singapore: World Scientific; 1998. 216 p. DOI: 10.1142/2284.
4. Rubin A, Riznichenko G. Mathematical Biophysics. New York: Springer; 2014. 273 p. DOI: 10.1007/978-1-4614-8702-9.
5. Cantrell RS, Cosner C. Spatial Ecology Via Reaction-Diffusion Equations. Chichester: John Wiley and Sons Ltd; 2003. 428 p. DOI: 10.1002/0470871296.
6. Malchow H, Petrovskii SV, Venturino E. Spatiotemporal Patterns in Ecology and Epidemiology: Theory, Models, and Simulation. New York: Chapman and Hall/CRC; 2008. 469 p.
7. Frisman YY, Kulakov MP, Revutskaya OL, Zhdanova OL, Neverova GP. The key approaches and review of current researches on dynamics of structured and interacting populations. Computer Research and Modeling. 2019;11(1):119-151 (in Russian).
DOI: 10.20537/2076-7633-2019-11-1-119-151.
8. Kim K, Choi W. Local dynamics and coexistence of predator-prey model with directional dispersal of predator. Math. Biosci. Eng. 2020;17(6):6737-6755. DOI: 10.3934/mbe.2020351.
9. Courchamp F, Berec J, Gascoigne J. Allee Effects in Ecology and Conservation. Oxford: Oxford University Press; 2008. 256 p. DOI: 10.1093/acprof:oso/9780198570301.001.0001.
10. Svirezhev YM. Nonlinear Waves, Dissipative Patterns and Catastrophes in Ecology. Moscow: Nauka; 1987. 368 p. (in Russian).
11. Tyutyunov YV, Titova LI. From Lotka-Volterra to Arditi-Ginzburg: 90 years of evolving trophic functions. Biol. Bull. Rev. 2020;10(3):167-185.
DOI: 10.1134/S207908642003007X.
12. Epifanov AV, Tsybulin VG. Modeling of oscillatory scenarios of the coexistence of competing populations. Biophysics. 2016;61(4):696-704. DOI: 10.1134/S0006350916040072.
13. Epifanov AV, Tsybulin VG. Regarding the dynamics of cosymmetric predator - prey systems. Computer Research and Modeling. 2017;9(5):799-813 (in Russian).
DOI: 10.20537/2076-7633-2017-9-5-799-813.
14. Holling CS. Some characteristics of simple types of predation and parasitism. The Canadian Entomologist. 1959;91(7):385-398. DOI: 10.4039/Ent91385-7.
15. Budyansky AV, Frischmuth K, Tsybulin VG. Cosymmetry approach and mathematical modeling
of species coexistence in a heterogeneous habitat. Discrete & Continuous Dynamical Systems -B. 2019;24(2):547-561. DOI: 10.3934/dcdsb.2018196.
16. Budyansky AV, Tsybulin VG. Modeling of multifactor taxis in a predator-prey system. Biophysics. 2019;64(2):256-260. DOI: 10.1134/S0006350919020040.
17. Abrams PA. The evolution of predator-prey interactions: Theory and evidence. Annu. Rev. Ecol. Syst. 2000;31(1):79-105. DOI: 10.1146/annurev.ecolsys.31.1.79.
18. Zagrebneva AD, Govorukhin VN, Surkov FA. Bifurcations in active predator - passive prey model. Izvestiya VUZ. Applied Nonlinear Dynamics. 2014;22(3):94-106 (in Russian).
DOI: 10.18500/0869-6632-2014-22-3-94-106.
19. Tyutyunov YV, Titova LI, Senina IN. Prey-taxis destabilizes homogeneous stationary state in spatial Gause-Kolmogorov-type model for predator-prey system. Ecological Complexity. 2017; 31:170-180. DOI: 10.1016/j.ecocom.2017.07.001.
20. Mishra P, Wrzosek D. The role of indirect prey-taxis and interference among predators in pattern formation. Math. Methods Appl. Sci. 2020;43(18):10441-10461. DOI: 10.1002/mma.6426.
21. Haskell EC, Bell J. Pattern formation in a predator-mediated coexistence model with prey-taxis. Discrete & Continuous Dynamical Systems - B. 2020;25(8):2895-2921.
DOI: 10.3934/dcdsb.2020045.
Цибулин Вячеслав Георгиевич - родился в Ростове-на-Дону (1956), окончил механико-математический факультет Ростовского госуниверситета (1978). Защитил кандидатскую (1990) и докторскую (2011) диссертации. Работает в Институте математики, механики и компьютерных наук имени И. И. Воровича Южного федерального университета. Заведует кафедрой теоретической и компьютерной гидроаэродинамики. Занимается вычислительной гидродинамикой, задачами конвекции, популяционной динамики и др. В соавторстве с В. Н. Говорухиным написал книги «Введение в Maple», «Компьютер в математическом исследовании: Maple, MATLAB и LaTeX».
Россия, 344090 Ростов-на-Дону, ул. Мильчакова, 8а
Институт математики, механики и компьютерных наук имени И. И. Воровича ЮФУ E-mail: [email protected]
Ха Данг Тоан - родился в Ханое (Вьетнам, 1982), окончил факультет математической педагогики Ханойского национального университета (2008). Защитил степень магистра алгебры и теории чисел в Университете естественных наук Национального университета Ханоя. Работает во Вьетнамско-Венгерском индустриальном университете, Ханой, Вьетнам и учится в аспирантуре Института математики, механики и компьютерных наук имени И. И. Воровича Южного федерального университета.
Россия, 344090 Ростов-на-Дону, ул. Мильчакова, 8а
Институт математики, механики и компьютерных наук имени И. И. Воровича ЮФУ
Vietnam-Hungary Industrial University
16, Huu Nghi Str, Son Tay Dist, Hanoi, Vietnam
E-mail: [email protected]
Зеленчук Павел Анатольевич - родился в Ростове-на-Дону (1981), окончил факультет физики Ростовского госуниверситета (2004). Защитил степень магистра математики (2020) Южного федерального университета. Является аспирантом в Институте математики, механики и компьютерных наук имени И. И. Воровича.
Россия, 344090 Ростов-на-Дону, ул. Мильчакова, 8а
Институт математики, механики и компьютерных наук имени И. И. Воровича ЮФУ E-mail: [email protected] ORCID: 0000-0001-6598-8521 AuthorID: 1092911
