МЕХАНИЗМЫ ФОРМИРОВАНИЯ ПАЧЕЧНОЙ ДИНАМИКИ В СИСТЕМЕ МИГРАЦИОННО СВЯЗАННЫХ СООБЩЕСТВ ТИПА ХИЩНИК-ЖЕРТВА Текст научной статьи по специальности «Математика»

Прикладные задачи

нелинейной теории колебаний и волн

Известия высших учебных заведений. Прикладная нелинейная динамика. 2023. Т. 31, № 2 Izvestiya Vysshikh Uchebnykh Zavedeniy. Applied Nonlinear Dynamics. 2023;31(2)

Научная статья УДК 517.925.42, 574.34

DOI: 10.18500/0869-6632-003030 EDN: FGDHHM

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

Е.В. Курилова, М. П. Кулаковш, Е.Я. Фрисман

Институт комплексного анализа региональных проблем ДВО РАН, Биробиджан, Россия E-mail: [email protected], [email protected], [email protected] Поступила в редакцию 22.09.2022, принята к публикации 26.12.2022, опубликована онлайн 27.02.2023, опубликована 31.03.2023

Аннотация. Цель работы — изучение периодических режимов динамики миграционно связанных неидентичных сообществ типа хищник-жертва, возникающих при частичной синхронизации колебаний их численностей. Комбинация этих колебаний приводит к образованию режимов динамики, объединяющих в себе как быстрые всплески численностей (пачечная динамика), так и медленные релаксационные колебания (тоническая динамика), которые характеризуются различным соотношением синхронной и несинхронной динамики в определенные периоды времени. В работе особое внимание уделено описанию сценариев перехода между разными типами пачечной активности. Эти типы отличаются между собой не столько размерами, формой и числом быстрых всплесков численностей в пачке, сколько очередностью появления этих всплесков относительно релаксационного цикла. Методы. При исследовании предложенной модели динамики миграционно связанных сообществ используются методы бифуркационного анализа динамических систем, а также геометрические методы, основанные на разделении полной системы на быстрые и медленные уравнения (две подсистемы). Результаты. Показано, что динамика первой подсистемы с релаксационным циклом определяет динамику второй — с пачечной динамикой — посредством гладкой зависимости режима от численности хищника и не гладкой зависимости от численности жертв. Были сконструированы инвариантные многообразия, на которых реализуются участки с тонической (медленное многообразие) и пачечной (быстрое многообразие) активностью рассматриваемой системы. Заключение. Описаны сценарии возникновения пачек с различной формой, которые определяются внешним видом быстрого инвариантного многообразия и расположения его частей относительно релаксационного цикла. Показано, что переходы между разными типами пачек сопровождаются сменой периода колебаний, степенью синхронизации, и в итоге система переходит к квазипериодической динамике, когда оба сообщества не синхронны между собой.

Ключевые слова: хищник-жертва, миграция, синхронизация, пачечная и тоническая динамика, быстро-медленная система.

Благодарности. Работа выполнена в рамках государственного задания Института комплексного анализа региональных проблем ДВО РАН.

Для цитирования: Курилова Е.В., Кулаков М.П., Фрисман Е. Я.Механизмы формирования пачечной динамики в системе миграционно связанных сообществ типа хищник-жертва//Известия вузов. ПНД. 2023. T. 31, № 2. С. 143-169. DOI: 10.18500/0869-6632-003030. EDN: FGDHHM

Статья опубликована на условиях Creative Commons Attribution License (CC-BY 4.0).

©Курилова Е. В., Кулаков М.П., Фрисман Е. Я., 2023

Article

DOI: 10.18500/0869-6632-003030

Mechanisms leading to bursting oscillations in the system of predator-prey communities coupled by migrations

E. V. Kurilova, M. P. Kulakovm, E. Ya. Frisman

Institute for Complex Analysis of Regional Problems, Far Eastern Branch, Russian Academy of Sciences, Birobidzhan, Russia E-mail: [email protected], [email protected], [email protected] Received 22.09.2022, accepted 26.12.2022, available online 27.02.2023, published 31.03.2023

Abstract. The purpose is to study the periodic regimes of the dynamics for two non-identical predator-prey communities coupled by migrations, associated with the partial synchronization of fluctuations in the abundance of communities. The combination of fluctuations in neighboring sites leads to the regimes that include both fast bursts (bursting oscillations) and slow oscillations (tonic spiking). These types of activity are characterized by a different ratio of synchronous and non-synchronous dynamics of communities in certain periods of time. In this paper, we describe scenarios of the transition between different types of burst activity. These types of dynamics differ from each other not so much in size, shape, and number of spikes in a burst, but in the order of these bursts relative to the slow-fast cycle. Methods. To study the proposed model, we use the bifurcation analysis methods of dynamic systems, as well as geometric methods based on the division of the full system into fast and slow equations (subsystems). Results. We showed that the dynamics of the first subsystem with a slow-fast limit cycle directly determines the dynamics of the second one with burst activity through a smooth dependence of regime on the number of predators and a non-smooth dependence on the number of prey. We constructed the invariant manifolds on which there are parts of dynamics with tonic (slow manifold) and burst (fast manifold) activity of the full system. Conclusion. We described the scenario for bursting with different waveforms, which are determined by the appearance of the fast invariant manifold and the location of its parts relative to the slow-fast cycle. The transitions between different types of burst are accompanied by a change in the oscillation period, the degree of synchronization, and, as a result, the dynamics becomes quasi-periodic when both communities are not synchronous with each other.

Keywords: predator-prey, migration, synchronization, tonic spiking, bursting, slow-fast system.

Acknowledgements. This work was carried out within the framework of the state targets of the Institute for Complex Analysis of Regional Problem of the Far Eastern Branch of the Russian Academy of Sciences.

For citation: Kurilova EV, Kulakov MP, Frisman EYa. Mechanisms leading to bursting oscillations in the system of predator-prey communities coupled by migrations. Izvestiya VUZ. Applied Nonlinear Dynamics. 2023;31(2):143-169. DOI: 10.18500/0869-6632-003030

This is an open access article distributed under the terms of Creative Commons Attribution License (CC-BY 4.0).

Введение

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

Изучение динамических моделей биологических систем, независимо от области их применения, основывается, как правило, на аппарате обыкновенных дифференциальных уравнений. При исследовании моделей нескольких взаимодействующих сообществ особый интерес вызывают нелинейные явления, возникающие при синхронизации динамики на разных территориях [4-6]

и связанные, например, с механизмами формирования пространственно-временной неоднородности, сложного динамического поведения и т.п. [7-12].

К одному из популярных направлений нелинейной динамики относятся исследования нейронной активности, которые описываются быстро-медленными системами. В таких системах часто возникают так называемые взрывные колебания [13-16]. За последние десятилетия были обнаружены и описаны разные виды таких колебаний, а также подробно изучены бифуркационные механизмы, приводящие к ним [17-21].

Настоящая работа входит в число таких исследований и посвящена изучению динамики численности двух смежных неидентичных миграционно связанных сообществ, каждое из которых является автоколебательной подсистемой и функционирует по принципу «хищник-жертва» с откликом Холлинга II типа [22,23]. Рассматриваемая система является модификацией уравнений Базыкина [24,25] — одной из базовых моделей популяционной биологии, в зарубежной литературе именуемой моделью Розенцвейга-Макартура [26,27].

Ранее выполненное нами исследование синхронизации регулярных колебаний, возникающих в рассматриваемой системе, выявило несколько особенностей [28-30]. Во-первых, в случае сильной связи возможна полная синхронизация циклов на разных территориях, даже при значительной разнице между сообществами [29]. Однако это приводит к единственно возможному типу динамики — предельному циклу. Во-вторых, снижение силы связи приводит к очень быстрой десинхронизации, при которой каждое сообщество испытывает колебания численности с собственным ритмом. В результате, в случае слабой связи синхронизация возможна лишь для идентичных сообществ, а мало отличающиеся сообщества неспособны к полной синхронизации. В целом, перечисленные результаты хорошо согласуются с результатами других авторов. Однако дальнейшее исследование показало, что слабосвязанные неидентичные сообщества способны, по крайней мере, к частичной синхронизации именно в случае большой разницы между значениями смертности хищников [28]. Кроме того в результате исследования были обнаружены режимы, сочетающие в себе как участки с быстрыми взрывными колебаниями (пачечной активности), так и с медленной тонической динамикой. Возникновение таких режимов происходит при изменении различий между рассматриваемыми сообществами в результате их частичной синхронизации. Это приводит к образованию сложно устроенной фазовой траектории, сочетающей участки синхронной и несинхронной динамики в разные промежутки времени. Степень синхронизации значительно влияет на изменение рисунка динамического поведения, что в свою очередь сказывается на изменении вида пачечной активности [30].

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

1. Модель динамики двух неидентичных миграционно связанных сообществ

Уравнения динамики численности двух соседних неидентичных миграционно связанных сообществ «хищник-жертва», описываемых уравнениями Базыкина [24,25] с функциональным откликом Холлинга II типа, имеют вид:

(1Хх

= АхХх

К -Хх ВХхУх

м

^ = -СП +

м

ЛХ2

К 1 + НХ\' БХ¥1 + М (У2 -У1)

= А2Х2

1 +АХ! К -Х2 ВХ2У2

(1)

§ = -С* +

К 1 + Н Х2

Ш + М к -*)

где Х\ и Х2 — численности популяций жертвы, У\ и У2 — численности популяций хищника в первом и втором сообществах, соответственно; К — устойчивая равновесная численность популяции жертв в каждом местообитании в отсутствие хищника; В — удельная скорость потребления популяцией хищника популяции жертвы при единичной плотности обеих популяций; С — естественная смертность хищника, Б/В — коэффициент переработки потребленной хищником биомассы жертвы в собственную биомассу; Н — коэффициент насыщения хищника; М — коэффициент миграции хищника. Неидентичность сообществ выражается в существенной разности максимальных скоростей воспроизводства популяций жертвы, обозначенных как А\ и А2, соответственно (А1 = А2). Это отражает ситуацию, когда на сопредельных территориях обитают два принципиально разных вида жертв, но имеющих одинаковую пищевую ценность для хищника.

После замены переменных: Хг = Схг/8, Уг = А1у1/Ь (г = 1,2), модель (1) преобразуется в следующую систему с пятью параметрами:

йХ1 Х1У1

— = Х1(1 - ОХ1) - ,

т 1 + пх1

йу 1

~Ж = -с т + 1 + ЛХ1

1 Х1 1 1

+ С1Ш — У2 2

йх2 /л ч Х2У2

— =Х2(1 — ах2) - 1+Х2,

{^У2- ш)

(2)

— = -с + ей 1 + Ь,Х2

Х

+ С2т\ —у1 1

{^-т - у2)

где Х% и уг — относительные численности жертв и хищника, к = НС/Б — коэффициент насыщения хищника, а = С/(К Б) — коэффициент самолимитирования жертвы, а = С/А^ — относительная скорость снижения (убыли) численности хищников (смертности) и тСг = М/Аг — коэффициент миграции хищников ( г = 1, 2). В статье [28] более подробно описана постановка задачи и приведено аналитическое исследование системы (1).

Надо заметить, что использованная в системе (2) замена переменных приводит к тому, что два неидентичных сообщества с разными скоростями роста численности жертв (А1 = А2) становятся подобными сообществам, в которых смертность хищников различна (с 1 = С2). Однако в этом случае связь между сообществами оказывается несимметричной (с 1т = С2т), то есть не совпадают темпы оттока особей, покидающих одну территорию (эмигранты), и темпы их притока на другую (иммигранты), что можно интерпретировать как гибель хищников в процессе миграции. Однако из-за разницы в смертности хищника (с 1 = С2) эта гибель будет различной при движении из первого или второго сообщества, что прямо указывает на наиболее «выгодное» направление для эмиграции: территория, где ниже смертность хищника или выше рождаемость жертв, что эквивалентно, согласно использованной замене переменных.

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

а = 0.1, к = 0.5, С1 = 0.002, С2 = 0.5, т = 0.005,

где значения а, к приводят к возникновению предельного цикла в каждой из подсистем, соответствующих разным сообществам. Значение С1 = 0.002 соответствует низкой смертности хищников на первой территории, а С2 = 0.5 — высокой смертности хищников на второй. Согласно приведенной замене переменных из (1) в (2), сообщество с низкой смертностью хищника имеет высокую скорость воспроизводства жертв, и наоборот, высокая смертность хищника присуща сообществу с медленной скоростью роста числа жертв. Такая ситуация возможна, когда рядом находятся местообитания жертв с принципиально разным репродуктивным потенциалом. Например, жертвами

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

В ходе исследования будем варьировать а и с\ при фиксированных Ъ = 0.5, с2 = 0.5, т = 0.005, следя за устойчивостью особых точек и изучая тип динамического поведения.

Аналитическое исследование показало, что система (2) имеет несколько особых точек:

• тривиальная Ао(0; 0; 0; 0) с нулевыми численностями всех популяций;

• полутривиальные А(1/а; 0; 0; 0), А2(0; 0; 1/а; 0), А3(1/а; 0; 1/а; 0) с нулевой численностью хищников;

„ (1 + 2т (1 + т)Ро С2т$2\ „ ств2 1 + 2т (1 + т)02\

• полутривиальные Щ —--;--; 0;-и В2 0;-~2Г; —--;-—

V Р1 С102 ' ^ С2Р2 Р1 С1Р2 )

(Р1 = 1 — Ъ + т — 2Ът, Р2 = Р1 — а — 2ат), соответствующие нулевым численностям одной популяции жертвы в разных сообществах;

• три нетривиальные еАх^, у^^х^, У2} = 0,1, 2) с ненулевыми численностями всех

популяций.

Были вычислены координаты всех точек, исследованы условия их существования и описаны простейшие бифуркации [28, 30]. Построение бифуркационных диаграмм выполнено с использованием открытой системы исследования динамических систем Ма1СоП [31], а также собственного набора программ. Для выделения области квазипериодической динамики использован метод карт показателей Ляпунова, а также метод карт динамических режимов, построенных для соответствующего отображения Пуанкаре. Результат этого исследования резюмирован на рис. 1.

На рис. 1, a линиями отмечены основные локальные бифуркации особых точек системы (1), а цветом выделено несколько областей, соответствующих разным типам динамики. Линия ТС соответствует транскритической бифуркации между точками В1 и Е1. Между линиями Н и Н- находится область существования асимптотически устойчивых предельных циклов (белая область), возникающих вокруг отмеченных в скобках особых точек. При пересечении линии Н- циклы теряют устойчивость и формируются релаксационные предельные циклы, известные также как быстро-медленный цикл, канард или уточное решение [21,32-34]. На рис. 1, d-e показаны примеры перехода предельных циклов ¿1 и 12 в релаксационные циклы С1 и С2, которые отличаются сильно возросшим размахом колебаний, а также структурой, показанной во втором столбце. На фазовых портретах видно, что в отличие от предельных циклов ¿1 и ¿2, циклы С1 и С2 «прилипают» к отмеченной ветви параболы (М^д) и одной из осей координат (Оу2), а также содержат участки с быстрым (две или более стрелок) и медленным (одна стрелка) движением фазовых точек.

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

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

(a) (b) (c)

H(Bi) TC(Bi,Ei) H(E0,Ei) I: TC(BbEl) H(E0),H+(Ei) II:

0 cx 0.015 ' 0 cx 0.004 ' 0.01 cx 0.015

Рис. 1. а — Бифуркационная диаграмма системы (2) и b — ее увеличенный фрагмент. Над линией Н устойчива точка В1 или Е1. Серая область — отсутствие колебаний (асимптотически устойчива точка А0); белая С — образование устойчивых предельных циклов; зеленая S — возникновение пачек; пурпурная Q — квазипериодическая динамика. Линия Н — бифуркация Андронова-Хопфа, TC — транскритическая бифуркация, HD — гомоклиническая бифуркация, SB — рождение пачечной динамики, KB и TB — возникновение квазипериодической динамики. c — карта динамических режимов в области II, где цветом закодирован период цикла отображения Пуанкаре. d-f — примеры предельных циклов в системе (2) (цвет онлайн)

Fig. 1. a — Bifurcation diagram of the system (2) and b — its enlarged fragment. Point B1 or E1 is stable above the line H. The gray domain is the absence of oscillations (the point Ao is asymptotically stable); white domain С is the region of stable limit cycles; green domain S corresponds to burst dynamics; magenta domain Q is quasi-periodic dynamics. Lines H, TC and HB are Andronov-Hopf, transcritical and homoclinic bifurcations respectively; SB is the birth of burst dynamics; KB and TB are the emergence of quasi-periodic dynamics. c — map of dynamic regimes in region II where color coded is the cycle period of the Poincare mapping. d-f — examples of limit cycles in system (2) (color online)

динамики численности имеют ромбовидную, треугольную (рис. 2, а), усеченную трапециевидную (рис. 2, Ь) или более сложные формы (рис. 2, с), и отличаются количеством витков на быстром пачечном многообразии. Число витков растет по мере приближения к области квазипериодической динамики ( (рис. 2, й-/).

В пурпурной области ( (см. рис. 1, а) модельная траектория плотно покрывает такие множества, как тор (при пересечении линии ТВ) (см. рис. 2, /) или бутылка Клейна (линии КВ) (см. рис. 2, й, е), а периоды колебаний сообществ на разных территориях связаны иррациональным отношением. Такой тип динамики ранее был рассмотрен через анализ сечений Пуанкаре [30]. В данном случае на специально подобранной секущей плоскости траектория формирует замкнутые инвариантные кривые, которые в случае динамики на торе или бутылке Клейна, имеют разное взаимное положение. Обнаружено, что в первом случае кривые вложены друг в друга, а во втором — расположены достаточно далеко друг от друга.

Далее по мере роста С1 наблюдается только квазипериодический тип динамики, и лишь при С1 ^ С2 происходит полная взаимная синхронизация. Хаотические режимы при принятых значениях параметров в данной модели не обнаружены.

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

10

(a) о=0.1, с1=0.003

■ 6

- V1

- \ Х2 \ .

- \ .

- 1^2 Jy1 .

0

2600

15

(b) a=0.08, с1=0.004

p'^wJiAftWliF^^jJ_WllM^iiitmmll^^t-.__I^HiiiiiwunmilMII

4700

10

(с) a=0.1, с1=0.0038

2200

t

(d) a=0.152, c1=0.00877

\

1. , 1 , , 1

У2 У1

3000 4800

5000

(e) a=0.13, c1=0.00877

U 1 IkYJ Г2) ii 1 tan, i 1 Ik

itf | , rim»1 1111 |тиг......iw limps Ittilllllllllllllliilli^Ninilllllll^itaHlll'lllilMilUllinllllllU1

5000 4700

У 2 5000

f) a=0.1497, c1=0.0143

0

ill

5000 4900

штш

t

5000

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

Fig. 2. Examples of bursting in system (2) from region S (fig. 1) which differ in the waveform of burst (color online)

0

t

t

8

0

t

t

8

0

Скорости изменения переменных у1 и у2 определены параметрами С1 и с2, причем с учетом выбранных значений параметров 0 < С1 ^с2< 1, следовательно можно принять, что 0 < С1 ^ 1 — малый параметр. Это позволяет рассматривать данную модель, как быстро-медленную систему. В том смысле, что первые два уравнения составляют медленную подсистему, в которой наблюдается только медленный цикл, а последние два — быструю, в которой возможны колебания с периодом много большим, чем у медленной подсистемы. Воспользуемся для исследования системы (2) методологией изучения быстро-медленных систем.

При исследовании таких динамических систем часто используются геометрические методы, основанные на разделении полной системы (2) на быстрые и медленные уравнения (подсистемы). При этом вариация медленной переменной, являющейся бифуркационным или сканирующим параметром, в предельном случае позволяет упростить поиск бифуркаций состояний равновесия, приводящих к рождению предельных циклов быстрой (двумерной) подсистемы, рассматриваемой в сингулярном пределе С1 = 0. Также изменение сканирующего параметра при малых С1 =0 позволяет построить медленное критическое многообразие, являющееся скелетом для различного типа активностей в полной быстро-медленной модели (2). Очевидно, что быстрая подсистема модели (2) обладает, по крайней мере, двумя такими многообразиями: многообразием покоя и тоническим многообразием, на котором возникают быстрые всплески численностей (спайки), показанные на рис. 2. Причем синергия быстрых и медленных взаимодействий может приводить к весьма нетипичным бифуркационным явлениям, возможным только в полной модели [17,35,36].

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

2. Бифуркации в медленной и быстрой подсистеме

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

йх 1 / л \ — = Х1(1 - аХ1) -

йу±

сИ

йх2 = йу2

а а

Х1У1

-ТГ = СП - ш(1 + т) +

0.

1 + к Х1 Х1 1

= /1(хъ Уl),

1 т +-У2

1 + к Х1

= Clgl(хl, yl,y2),

(3)

В системе (3) присутствует малый параметр сГ (0 < сГ ^ с2 < 1), который входит во второе уравнение дважды. Наличие малого Сг приводит к тому, что влияние «замороженной» переменной у2 на систему значительно слабее, чем уГ и хГ.

Для второй подсистемы аналогично можно записать:

' йхг = й/ут = 0

м а '

Х Х

< ~Ж = Х2(1 -аХ2) - ТГлХ;

= ( Х , ),

(4)

йУ2 . С2Х2у2 (с2 \ , ,

— = -С2У2 + 1 , +с2т\— У1 -у2\ = д2(х2,у2, У1).

й

1 + к Х

Система (4) существенно зависит от параметра уГ, поскольку с^т/сГ ^ сГ.

Назовем систему (3) ведущей или сильной, а (4) — ведомой или слабой. Покажем, что сильная система (3) определяет динамику (4), а также совместную динамику полной системы (2).

Для начала рассмотрим бифуркации в подсистеме (3) при изменении параметра у2. Состояния равновесия данной подсистемы лежат в плоскости (х\, у\), на пересечении двух нульклин Л(Ж1, ш) =0 и (л(х1, ух,у2) = 0, где у2 — параметр.

Как известно [17,33], для систем вида (3) в предельном случае с1 =0 уравнение /1(х1, у1) = 0 задает в фазовом пространстве критическое многообразие или медленную кривую: Мед = = { (х1, у1) Е Я1 /1 (х1, у1) = 0} и {х1 = 0, у1 ^ 1}. Заметим, что она становится поверхностью для полной системы (2). Эта кусочно-заданная кривая содержит устойчивую М^д (верхняя часть параболы на рис. 3, a) и неустойчивую М^ ветви (нижняя часть параболы и часть прямой х1 = 0). На рис. 3, a эти ветви показаны сплошной и пунктирной линиями, и выступают каркасом релаксационного цикла в ведущей системе (3) при достаточно малых с1 = 0. Как показано на рис. 3, а одна часть цикла лежит на устойчивой ветке многообразия М^д, другая на .

Точка (ъ) на рис. 3, а и рис. 5 соответствует центру «масс» предельного цикла (среднее значение координат) и находится из стандартной формулы [36]:

т

(ъ) = 1 У ф(1)(И, 0

Рис. 3. a — Нульклины системы (3) с нанесенным релаксационным циклом системы (2) — канард С1. Стрелочки здесь и на рис. 3, с — направления движения вдоль цикла, их количество — скорость движения. b и c — Бифуркационная диаграмма системы (3) и (4) при вариации параметра у2 и у-^, соответственно. Серой полосой выделен диапазон значений переменных цикла С\, попадающих в разные части диаграммы b и c. d — Динамика полной системы (2), соответствующая циклу С\ (цвет онлайн)

Fig. 3. a — Nullcline of system (3) with a slow-fast limit cycle of system (2) (canard Ci) lying on top of nullcline. b and c — Bifurcation diagram of system (3) and (4) with variation of the parameter y2 and y1, respectively. The gray bars highlight the range of the limit cycle C1 falling into different parts of the diagram b and c. d — Dynamics of the full system (2) corresponding to cycle C1 (color online)

где z = ф(£) — уравнение, описывающее изменение координат предельного цикла с периодом Т. Для системы (3) z = (х\,уi), для (4) z = (х2, у2). При вычислении этого интеграла рассматриваются численно полученные координаты точек предельного цикла на одном витке, выполняется суммирование их координат, а результат делится на число точек, лежащих на цикле.

На рис. 3, a показан канард, возникающий в системе (3) при фиксированном значении «замороженной» переменной у2. Более того, этот цикл с точностью до малых порядка 0(с\) совпадает с периодическим решением полной системы (2) при указанных значениях параметров. Число стрелок на рис. 3, a указывает на скорость движения вдоль канарда таким образом, что в горизонтальном направлении цикл содержит медленное движение, а в вертикальном — быстрое.

При малом сi=0 вариация параметра у2 двигает вторую нульклину, определяемую уравнением д\(х\,у\, у2) = 0, относительно неподвижной первой, не зависящей от значений у2. Вместе с ней сдвигается состояние равновесия. При сi = 0 нетривиальная стационарная точка z = (Х\, у\) совпадает с вершиной параболы f\(x\,yi) = 0, а при с\ > 0 смещается внутрь цикла. На рис. 3, b показана зависимость координат (Х\) состояний равновесия системы (3) от у2. При некотором значении у2 в системе (3) происходит бифуркация Андронова-Хопфа (точка H), и из стационарной точки рождается устойчивый предельный цикл. С ростом параметра он увеличивает свой размер (амплитуду колебаний). Чтобы это показать, на рис. 3, b приведена зависимость крайних значений численности жертв (max(xi) и min(xi)) медленной подсистемы (3), а также среднего значения координат цикла {х\) при изменении у2. В точке, обозначенной как «canard explosion», цикл стремительно увеличивается в размере и «прилипает» (стремится) к соответствующим ветвям медленной кривой Meq, а центр цикла значительно отдаляется от стационарной точки. В результате, фазовый портрет подсистемы (3) при допустимых значениях у2 > 0 представляет собой канард, который практически не меняется даже при значительных вариациях «замороженной» переменной у2.

В отличие от сильного сообщества (3), слабое сообщество (4) существенно зависит от параметра у\. Рассмотрим локальные бифуркации системы (4) в зависимости от у\. Координаты стационарных точек Х2 и У2 несложно вычислить как решение системы

/2(^2, У2) = 0,

92(Х2, У2,Уl) = 0.

На рис. 3, c показана зависимость стационарных точек от параметра у\. Эта зависимость представлена прямой M0 = {(у\, у2) Е R2\у2 = с2ту\/(с\(\ + ш)),х2 = 0} и кривой M}q = = {(У1,У2) Е R2\(C2mhy 1 -ci((1 + m)h- 1)у2)2 -сю2т(а + h)yi + с2((1 + т)(1 + h) -1)У2 = 0, х2 > 0}. Они пересекаются в точке LP, соответствующей транскритической бифуркации, в которой совпадают координаты стационарных точек и происходит обмен устойчивостью. Справа от нее устойчивы точки, расположенные на прямой, но с нулевой стационарной численностью жертв слабого сообщества (х2 = 0). Слева устойчивы точки на параболе с ненулевой стационарной численностью жертв (х2 > 0). Во второй точке пересечения совпадают лишь численности хищника, но не жертв. Поэтому качественных изменений там не происходит. Покажем, что множества M0 и M^, так же как Meq, можно рассматривать в качестве приближения части инвариантных многообразий полной системы (2).

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

Во-первых, если вариация у\ лежит правее точки ЬР, то динамика фазовой переменной у2 строго следует за всеми изменениями ух (см. рис. 3, с). Как уже было показано, ведущая система может генерировать лишь периодические колебания, поэтому динамика численностей хищников на обоих участках оказывается синфазной (см. рис. 3, й). Численность жертв х2 в ведомой системе в асимптотическом случае оказывается нулевой. Данный предельный цикл (канард Сх) формируется вокруг полутривиальной стационарной точки Вх. На рис. 3, а показано, что при проецировании на плоскость (хх, у\) цикл стремится к ветвям многообразия Мед, а при проецировании на плоскость (ух, у2) — лежит на участке прямой М0д (см. рис. 3, с).

Во-вторых, параметр у\ может изменяться так, что фазовые точки расположены слева от точки ЬР на диаграмме. В этом случае цикл в ведущей системе провоцирует в ведомой системе противофазные к ней колебания с ненулевой численностью жертв х2. Пример такой динамики показан на рис. 4, а, где аттрактор полной системы (2) лежит только на участке параболы М^. Справа показано, как во времени изменяются значения фазовых переменных у\ и у2, а также расположение фазовой кривой на одной из ветвей медленной поверхности при проецировании в трехмерное подпространство (хх,х2, у2). В полной системе (2) этот цикл формируется вокруг нетривиальной точки Е\.

В-третьих, вариация у\ может попадать как на левую, так и на правую от точки ЬР часть бифуркационной диаграммы. В результате цикл ведущей системы будет генерировать в ведомой системе чередующиеся участки синфазной и противофазной динамики. Пример такого поведения показан на рис. 4, й-/. В период синфазной динамики фазовая кривая лежит на поверхности плоскости х2 = 0, а в период противофазной — вне этой плоскости (см. рис. 4, й), а также частично на медленной поверхности Ми на неизвестном многообразии, которое предстоит найти (см. рис. 4,/). В этом случае, очевидно, существует петля гистерезиса такая, что увеличение или уменьшение значения у\ приводит к движению фазовых точек вдоль разных стационарных

1

У2

0 -0.5

1

У2 0

У1 1

(d)

H ф C+S[

•SN LP \ Mlq

к 1— |\

a=0.164, с1=0.0035

\ X1 -1— т- (b) [Y

У2

У1 L

7000 t 11000 x

a=0.163, ^=0.00285 y1

-1- V1 \ (e) \

У1 У2

% NX2I CT4 K. -

16000

20000

-0.5 0

У1 1

2

/

/Х2

S1

Рис. 4. a, d — Бифуркационные диаграммы системы (4) при изменении параметра у\, на которые нанесены фазовые точки периодического решения системы (2). b, e — Примеры динамики системы (2). c, f — Расположение фазовых точек относительно критических многообразий (цвет онлайн)

Fig. 4. a, d — Bifurcation diagrams of system (4) with a change in the parameter y\ on which the phase points of the periodic solution of system (2) are plotted. b, e — Examples of dynamics for system (2). c, f — Location of phase points with respect to critical manifolds (color online)

2

5

0

0

2

2

5

t

состояний ведомой системы (ветвей), пересекающихся в точке ЬР (см. рис. 4, й). Когда значение у\ снижается, то происходит движение вдоль прямой М^, в том числе левее точки ЬР. Затем может возникнуть небольшой всплеск численностей жертв и хищника на втором участке, очевидно не связанный с бифуркациями в ведомой системе, поскольку он наблюдается вдалеке от точки Н. Траектория при этом выходит на параболу М^, а переменная у\ переходит к росту, пока не достигнет точки ЬР. Затем она опять возвращается на прямую М°. Движение фазовой траектории в этом случае происходит по часовой стрелке вдоль треугольника, образованного двумя ветвями многообразия.

В последних двух случаях характер динамики в противофазной части цикла существенно меняется (возникает пачечная активность) в зависимости от того, насколько близко к точке Н попадает вариация , а также от происходящих в ведомой системе (4) других качественных изменений. В примерах, приведенных выше, предельный цикл ведомой системы (4) рождается (исчезает) в двух точках Н, лежащих на устойчивой ветке стационарных состояний . Однако при достаточно малых с\ размах колебаний в ведущей системе (3) самостоятельно никогда не достигает этой точки. Поэтому при достаточно малых с\ (параметров из области С на рис. 1, а) пачечная активность на основе цикла ведомой системы может быть лишь частью переходного процесса, если начальное значение у 1(0) расположено между двумя точками Н бифуркации Андронова-Хопфа (см. рис. 4, а или й).

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

На приведенных диаграммах с помощью кривых шах(у2) и шт(У2) (см. рис. 5) можно проследить, как меняется размер предельного цикла ведомой системы (4) при разных значениях параметров а и с\. Указанный подход составляет основу так называемой техники продолжения по параметру. Ветвь кривой, обозначенная через (У2) на рис. 5, показывает зависимость средних координат предельного цикла. Заметим, что при уменьшении у\ средняя ветвь (У2) приближается к нижней шт(У2).

Таким образом, на верхней ветке стационарных состояний М^ ведомой системы (4) образуется предельный цикл, который можно рассматривать в качестве каркаса пачечной динамики полной системы (2). При изменении величины у\ меняется размер и расположение предельного цикла ведомой системы (4). Несложно понять, что в фазовом пространстве полной системы (2) каждый из этих циклов лежит на гладкой поверхности параболической формы, исходящей из точки Н. Эту поверхность можно рассматривать как приближение части инвариантного многообразия системы (2), на котором располагаются быстрые всплески численности хищника и жертв (быстрое многообразие). По крайней мере, той его части, на которой лежат затухающие быстрые периодические движения. Обозначим его через Мро. Наличие этого многообразия позволяет сразу объяснить механизм рождения пачек треугольной формы, показанных на рис. 5, а. В случае таких пачек, сначала траектория медленно движется вниз вдоль прямой МО (многообразие покоя). Затем, достигнув минимально возможного значения у\, траектория «перескакивает» на быстрое многообразие Мро (левее точки Н). Этот скачок сопровождается стремительным ростом переменной х\ . На этом многообразии траектория совершает несколько витков, соответствующих только затухающим колебаниям фазовых переменных. Обозначим эту часть динамики как Б. Достигнув точки Н, траектория покидает быстрое многообразие и оказывается на многообразии покоя . При переходе через точку ЬР траектория переходит на прямую МО, и цикл замыкается.

4 3

У2 1

0

4'

3

У2 2-

(a)

4 3

2

У2 1

0

max(y2X^

SN < miny) ■ HB' lb ' H LP VI ей.

9 1'чвнр 'Tel1 I 'in

II . шаШ 1 ' d

f'i i I ■ Tig' i i iFTKi i> Tjgi ■

-2 -1 0 y1 1

Рис. 5. Слева — бифуркационные диаграммы системы (4) при изменении параметра yi, на которые нанесены фазовые точки периодического решения системы (2), в центре — примеры динамики системы (2), справа — вид аттрактора и его расположение относительно многообразия покоя (цвет онлайн)

Fig. 5. On the left — the bifurcation diagrams of system (4) with a change in the parameter yi and phase points of the periodic solution of system (2). In the middle — examples of the dynamics of system (2). On the right — the view of the attractors and its location relative to the slow manifold (color online)

0

Во втором и третьем примерах на рис. 5, Ь-с возникают пачки усеченной ромбовидной формы, что указывает на немного иной механизм формирования быстрых колебаний и существование еще одной ветки быстрого многообразия. Кроме того, переход с быстрого многообразия Мро на многообразие покоя Мпроисходит значительно правее точки Н во всех случаях. Во втором примере на быстром многообразии можно выделить два участка динамики и , соединенные канардом С (см. рис. 5, Ь). Участок состоит из затухающих быстрых колебаний при максимальных значениях х\ ^ 0, — из расходящихся колебаний при минимальных значениях х\ ~ 0 (см. рис. 5, Ь).

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

При низких или даже близких к нулю значениях х\ колебаний нет или рождаются расходящиеся колебания 52, а величина у\ только падает (см. рис. 5, Ь). Из этого можно предположить, что быстрое многообразие Мро не единственное или, по крайней мере, оно имеет дополнительную ветку, на которой реализуются быстрые расходящиеся колебания.

Данное наблюдение было сделано ранее, в работах [28,30], однако там не было доказано наличие второй ветки быстрого многообразия. Для проверки данного утверждения выделим еще одну подсистему, в которой величина х\ является постоянной величиной:

сгд1(х1,у1, у2),

(5)

д2(х2,У2,У 1).

В системе (5) также присутствует малый параметр с\ (0 < с\ ^ с2 < 1). Однако, из-за нелинейной зависимости производной у\ от значения «замороженной» переменной х\, влияние х на динамику и бифуркации системы (5) должно существенно отличаться от рассмотренных зависимостей поведения систем (3) и (4) от переменной .

Проследим за типичными бифуркациями системы (5) при значениях параметров, соответствующих пачечной активности в системе (2) с пачкой сложной формы, состоящей из усеченного ромба и треугольника, при а = 0.11, с1 = 0.0035 (остальные параметры аналогичны использованным выше). В данном случае в полной системе (2) присутствуют затухающие колебания на многообразии Мро, а также серия затухающих и расходящихся колебаний на искомом инвариантном множестве Мро.

Система (5) имеет единственное нетривиальное состояние равновесия Е*(у\,х2, У2), которое можно найти как решение системы уравнений:

91(х1, У!, У2) = 0, ¡2(Х2, У 2) = 0,

92(т ,х2, т) = 0.

Характер устойчивости состояний равновесия системы (5) показан на рис. 6, а, где сплошная линия соответствует зависимости нетривиальной стационарной численности т от значения параметра х\, а пунктир — тривиальному состоянию равновесия (0, 0, 0). В точке Н+ происходит суперкритическая, а в Н- — субкритическая бифуркация Андронова-Хопфа. Слева от Н+ кривая у1(х1) приближается к у\ =0 с положительной стороны, справа от Н- она приближается к т = 0 из отрицательной области. Данная зависимость имеет разрыв в точке х* = (к + 2тк — 1 —т)/((1 — к)(к + 2тк — 1)), которую вместе с точками Н+ и Н- можно рассматривать как некоторые пороговые значения х\, определяющие тип динамики. Кроме того, система (5) имеет полутривиальное состояние равновесия (0,1/а, 0), координата у\ = 1/а которой лежит выше данных линий и не показана на рис. 6, а. Сепаратрисы, исходящие из тривиальной и полутривиальной точки, порождают некоторые области притяжения циклов в системе (5).

Несложно показать, что правее точки х* производная йу\/(И > 0 для любых х2 > 0 и У2 > 0, и не важно, монотонно ли они изменяются или колеблются. В результате, если х\ > х*, то переменная у \(Ь) может только неограниченно расти при у 1(0) > 0 (вернее, для начальных условий выше сплошной на рис. 6, а). Но поскольку в системе (5) присутствует малый

ахх = 0

М = ,

(У\ ( ,л , ч , хш слт

^ = С1 ^+т) + тгкк + Т и =

х2 х2 2

— = х2(1 — ах2) — = к(х2, у2),

2 2 х2 2

^ = —с 2У2 + ГГ2х2 + С2т[~Ш—У2

г-У2)

Рис. 6. a — Бифуркационная диаграмма системы (5) при а = 0.11, с1 = 0.0035 и эскизы фазовых портретов. b — Переходная динамика системы (5) при движении к устойчивому предельному циклу Ls и c — при при движении от неустойчивого предельного цикла Lu. d — Приближение быстрых многообразий системы (2) динамикой (циклами) системы (5) (цвет онлайн)

Fig. 6. a — Bifurcation diagram at a = 0.11, ci = 0.0035 and phase portraits sketches of system (5). b — Transient dynamics of system (5) when moving towards a stable limit cycle Ls and c — from an unstable limit cycle Lu. d — Approximation of fast manifolds of system (2) by dynamics (cycles) of system (5) (color online)

параметр с\, то значения у\ изменяются настолько медленно, что, прежде чем уйти из окрестности неустойчивого цикла Ьи, фазовая кривая совершает большое число витков в его окрестности. Расстояние между витками медленно увеличивается по мере роста у\. Со временем траектория стягивается к лучу, и все фазовые переменные начинают неограниченно и монотонно возрастать (рис. 6, с). Можно показать, что с точностью до малых порядка 0(с±) этот луч совпадает с медленным многообразием или М0 системы (4), поскольку эти кривые совпадают с медленными кривыми системы (5) при с 1=0 и, очевидно, сохраняются при малом возмущении 0 < с 1 < 1.

На рис. 6, с показана динамика системы (5) при х\ = 3.5, где стрелки показывают направления движения траектории. Из рисунка видно, что траектория наматывается на поверхность параболической формы Мро, вид которой не зависит от начальных условий. В данном случае не важно, откуда начинать интегрирование (снаружи или внутри этой фигуры) или насколько далеко от цикла Ьи. Траектория быстро притягивается к этой поверхности, а затем медленно по ней движется, пока не стянется к лучу М\ .

При х\ < х* производная йу\/(И < 0 для любых х2 > 0 и у2 > 0. Поэтому, если х\ < х*, то переменная у\ (¿) только падает до своего стационарного значения, определяемого зависимостью У\(х\). Правее точки Н+ рождается устойчивый предельный цикл Ь3, и переменные начинают испытывать колебания. По мере достаточно медленного снижения величины у\ (из-за малого параметра с{) фазовая траектория наматывается на поверхность, имеющую форму однополостного гиперболоида, а не параболоида, как можно было ожидать (рис. 6, Ь). Переход траектории на эту поверхность происходит следующим образом. При малых 0 < х\ < х* существует некоторое критическое значение переменной у* такое, что при у\ > у* переменные у\, х2 и у2 достаточно монотонно падают, а траектория системы (5) движется по М0. При достижении и преодолении этого значения, то есть при у\ < у*, в системе (5) происходит скачок величин х2 и у2, что приводит к возникновению колебаний. По мере падения значений у\ в системе (5) сначала возникают затухающие колебания х2 и у2, которые при любых начальных условиях не достигают кривой М^д, а затем они сменяются расходящимися. Двигаясь в сторону устойчивого предельного цикла Ь3, траектория совершает большое число витков вдоль поверхности Мр0, расстояние между которыми сокращается по мере приближения к циклу Ь3. Более того, численные расчеты показывают, что часть динамики полной системы (2), соответствующая затухающим и расходящимся колебаниям, достаточно точно совпадает с двумя приведенными на рис. 6, й решениями системы (5). Например, если начальная точка систем (5) и (2) лежит на медленном многообразии М0 при у\(0) > у* и х\ =0, то их динамика с точностью до малых совпадет в части расходящихся колебаний на Мр0 до момента перехода траектории системы (2) на многообразие Мро, в то время как траектория системы (5) продолжит двигаться к циклу Ь3. Аналогично для начальных условий, взятых на многообразии Мро, динамика (5) и (2) совпадает до малых даже при переходе на ветки медленного многообразия М}д и М0.

Таким образом, инвариантные многообразия системы (5) можно рассматривать как более точное приближение быстрого многообразия системы (2), образованное участками траектории системы (5) при движении к предельному циклу Ь3 или движении от цикла Ьи. Они представляют собой две двумерные «пересекающиеся» поверхности: параболоид Мро и однополостный гиперболоид Мр0.

Для приближенного построения параболической поверхности Мро (рис. 7, Ь) достаточно рассмотреть несколько траекторий с начальными точками из разных частей неустойчивого цикла Ьи и выразить координаты фазовых точек через цилиндрические координаты. Затем отсортировать полученные координаты фазовых точек по росту (или падению) величины у\ и произвести интерполяцию полученных точек, которую будем использовать в качестве параметризации искомой поверхности. На рис. 6, с приведена одна из используемых траекторий. Каждая дополнительная траектория будет заполнять пробелы в местах, где увеличивается расстояние между витками. Параметризация гиперболоида Мр0 (рис. 7, а) выполняется похожим образом. Однако в качестве стартовых точек подбираются такие значения у 1(0) и у2(0) на медленной кривой М0 (х2(0) = 0), при которых траектория системы (5) максимально приближается к медленной кривой М^д в месте перехода от затухающих колебаний к расходящимся. Одновременно с этим необходимо, чтобы переход от монотонной динамики к колебаниям (скачок с кривой М0 на поверхность Мр0) происходил при максимально возможном значении у\. Тем самым оценивается величина у*.

На рис. 7, а-с показан вид полученных поверхностей. Несложно заметить, что параболоид Мро пересекает медленную кривую М^д в своей вершине (точка Н). В то время как для поверхности Мр0 кривая М^д выступает своего рода осью «симметрии» и не пересекается с ней. Вместе с тем однополостный гиперболоид Мр0, по всей видимости, касается кривой М0 в точке у*, соответствующей переходу от монотонной динамики к колебаниям в системе (5). Кроме того, на рис. 7, а-с продемонстрировано, что участки динамики системы (2) лежат на данных поверхностях. В этом случае траектория полной системы (2) разделена на две части. В первую

Рис. 7. a-c — Две ветки быстрого многообразия Мрэ0 и Мро, на которых лежат участки траектории системы (2), a — с расходящимися колебаниями S2, для которых xi ~ 0, b — с затухающими колебаниями Si, для которых xi ^ 0. d — Пример поведения системы (2) с выделенными участками динамики. e-f — Расположение траектории относительно медленных многообразий M0q и Меq (цвет онлайн)

Fig. 7. a-c — Two branches of fast manifold Мрэ0 and Мро on which lie parts of the trajectory of system (2) a — with divergent oscillations S2 (xi ~ 0) and b — with damped vibrations Si (xi ^ 0). d — The example of the behavior of system (2) with selected parts of dynamics. e-f — Location of the trajectory with respect to slow manifolds M0q and Меq (color online)

часть попадают участки динамики, для которых значения х\(£) ^ 0, во вторую часть — значения х\(£) » 0. В результате первая часть траектории лежит на Мро и М0д (см. рис. 7, а). Вторая — на МР0, Мед и М0д (см. рис. 7, Ь).

Заметим, что при проектировании в другие подпространства быстрые многообразия Мро и МРо не пересекаются так, как показано на рис. 7, с, и имеют немного другой вид. В подпространстве (х1, у\, у2) Мро и МРо — двумерные кривые, которые лежат на цилиндрической поверхности с основанием М^д и М]!д (хр > 0 и у2 > 0) (рис. 7, е), образующим каркас релаксационного цикла в системе (3) (см. рис. 3, а). Соответственно, быстрые колебания также лежат на данном цилиндре. В подпространстве (хр ,х2, у2) МРо — сильно сплюснутый параболоид, а Мро полностью сливается с плоскостью хр =0 (рис. 7,/). В обоих случаях между двумя быстрыми многообразиями расположено многообразие покоя, по которым траектория совершает медленное движение по поверхности цилиндра с основанием М^д и Ме/д (см. рис. 7, е).

3. Сценарий образования пачек

Наличие разных типов активности в подсистемах (3), (4) и (5) (монотонный рост или падение численностей, а также затухающие или расходящиеся колебания) позволяет описать сценарий формирования пачечных режимов колебаний в полной модели (2). Этот сценарий построен на анализе взаимного расположения быстрых и медленных многообразий, а также на наблюдениях за фазовой траекторией, которая способна периодически (или нерегулярно) перескакивать между разными аттракторами медленной и быстрой подсистемы. Сначала фазовая точка движется только по релаксационному циклу (канарду) С, соответствующему фазовой траектории первой (медленной) подсистемы (3). Затем траектория притягивается к одному из периодических решений быстрой подсистемы (4) или (5) (при разных х\). В этом случае решения системы (2) при сi = 0 сходятся к аттрактору быстрой подсистемы, которым является ее устойчивое состояние равновесия или предельный цикл. Когда аттрактор быстрой подсистемы является экспоненциально устойчивым, то он зависит гладко от у\. Однако гладкость сохраняется лишь при х\ = х*. Изменение у\ при некоторых фиксированных х\ приводит к образованию гладкого притягивающего инвариантного многообразия, которое состоит из медленной и быстрой части. Перечислим их.

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

• Кривые, образованные проекцией кривых М3д и М^ на двумерные поверхности М0 и М^д, которые являются состояниями равновесия слабого сообщества (4) и (5) при сi = 0.

• Трехмерная параболическая поверхность Мро и однополостный гиперболоид Мр0 (см. рис. 5, a и рис. 7), которые образованы циклами слабого сообщества (5).

Первые две части многообразия соответствуют медленной тонической динамике или состоянию покоя (rest state) обоих сообществ. Как было показано выше, из-за наличия более одной ветки возможно несколько вариантов расположения такого аттрактора (канарда) в фазовом пространстве. Например, канард может лежать только на ветке М0 (С), тогда динамика хищников на обеих территориях будет синхронной, но с нулевой численностью жертв во втором сообществе. Либо он может лежать только на М^дд (С2), тогда динамика будет противофазной с ненулевой численностью жертв. Либо канард может попеременно (регулярно) располагаться на разных ветках, тогда динамика сообществ будет содержать периоды синхронной и несинхронной динамики, совпадающие с периодами околонулевой (х\ ~ 0) и ненулевой (х\ ^ 0) численности жертв во втором сообществе (С).

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

В первом случае (х ^ 0), если у\ ниже определенной величины (левее точки H на рис. 7,b), траекторию перебрасывает на параболическую поверхность Мро, где она совершает несколько витков. Содержательно такое динамическое поведение может соответствовать усилению конкурентных взаимоотношений между хищниками слабого (местные особи) и сильного сообщества (иммигранты). В результате сильная конкуренция приводит к быстрым всплескам численности в слабом сообществе, которое проявляется в виде затухающих колебаний численностей хищников и жертв на этой территории. Эти затухающие колебания продолжаются до тех пор, пока число хищников сильного сообщества растет до определенной величины (точка H на рис. 7, b). Достигнув этой величины, траектория мягко возвращается в состояние покоя.

Во втором случае, когда количество жертв сильного сообщества х\ близко к нулю (х ~ 0), число хищников у\ падает. Если при этом у\ достигает определенного значения у*, то траекторию резко перебрасывает на гиперболоид Мр0, на который она наматывается. Аналогично прошлому случаю, конкуренция иммигрантов из сильного сообщества с местными особями инициирует быстрые расходящиеся колебания численностей в слабом сообществе. Эти колебания продолжаются до достижения минимально возможной численности хищника у\, при которой происходит стремительное восстановление жертв х\ сильного сообщества и медленное восстановление хищников. В самом крайнем случае это происходит в месте «пересечения» гиперболоида Мр0 с параболоидом Мро в подпространстве (у\,х2, у2), то есть при у\ ^ 0. После этого траекторию жестко перебрасывает с гиперболоида на параболоид, и колебания сменяются затухающими на фоне роста числа хищника у\ сильного сообщества.

При вариации параметров системы (2) меняется как размер релаксационного предельного цикла, так и расположение точек бифуркации в подсистемах (3), (4) и (5), которые выступают «краями» быстрого многообразия. В результате разным типам пачечной активности с разной формой пачки соответствует различное взаимное расположении канарда и поверхностей Мро и Мр0. На рис. 8 схематически показан механизм (сценарий) формирования пачек разной формы. Для этого две ветки быстрого многообразия Мро и Мр0 показаны вложенными друг в друга и соединены предельным релаксационным циклом С, а точки, где они соединяются, выделены звездочкой.

На рис. 8, а показано, что если асимптотически устойчив только релаксационный цикл, то поверхности Мро и Мр0 существуют, но расположены в стороне от него. В результате для параметров из области С на рис. 1, а, где в предельном случае в системе (2) возможен только релаксационный цикл, быстрые колебания возникают при некоторых возмущениях начальных условий, но как часть переходной динамики. Например, если стартовую точку взять в окрестности

(d)

Mpo

'PO

^Si

i CJ

a=0.16, Cj=0.004

a=0.07, c1=0.003

a=0.085, c1=0.0035 a=0.08, c1=0.0045

0

C C

Si Si

27000 t 30000 21000 t 25000 23000 t 25000 28000 t 30000

Рис. 8. Первая строка — сценарий формирования пачек разной формы в системе (2), связанный с различным расположением канарда С и быстрых многообразий Мро и Мр>о, где Si для xi ^ 0, S2 для xi ~ 0. Вторая строка — примеры динамики, соответствующие этим сценариям: a — предельный релаксационный цикл (канард), b — треугольная пачка, c — усеченная ромбовидная пачка и d пачка, состоящая из усеченного ромба и треугольника (цвет онлайн)

Fig. 8. The first line is the scenario for the formation of various waveforms bursting in the system (2), associated with different positions of the canard С and fast manifolds. The second line is examples of dynamics corresponding to these scenarios: a — slow-fast limit cycle (canard), b — triangular-wave bursting, c — truncated rhomboid-wave bursting and d — bursting consisting of a truncated rhombus and a triangle (color online)

5

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

В случае треугольной пачки канард С касается поверхности Мро, и траектория начинает самостоятельно двигаться по нему (рис. 8, Ь). Это происходит лишь после полного восстановления численности жертв хр в сильном сообществе, то есть при хр ^ 0. В случае ромбовидной и усеченной ромбовидной пачки канард дополнительно касается Мро (рис. 8, с-й), что происходит задолго до восстановления численности жертв хр, то есть при хр ~ 0. Но касание канарда с гиперболоидом Мро возможно в разных местах. Если траектория «пересекает» параболоид Мро, а затем достигает гиперболоид Мро, возникают усеченные ромбовидные пачки (см. рис. 8, с). Если же траектория переходит на гиперболоид не пересекая параболоид, то возникает более сложная пачка, в которой затухающие колебания сменяются расходящимися, а затем вновь затухающими (см. рис. 8, й). Естественно, никакого пересечения траектории и быстрого многообразия в четырехмерном пространстве не происходит. Это «пересечение» соответствует скачку численностей хищников ур и жертв хр, что хорошо заметно при проецировании в другие подпространства (см. рис. 7, е и/).

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

Надо отметить, что сценарий на рис. 8 описывает случай полностью периодической динамики, когда после каждого полного оборота фазовая точка строго возвращается в ту же точку фазового пространства (где численности достигают строго тех же значений), например, в точку, отмеченную звездочкой. В этом случае время между появлением всплесков численностей (расстояние между пачками) оказывается постоянным на каждом полном обороте цикла. В случае квазипериодической динамики (область Q на рис. 1) траектория на каждом своем витке возвращается в малую окрестность исходной точки. В результате точки соединения канарда с одной из веток быстрого многообразия образуют замкнутую инвариантную кривую. Несложно понять, что образованная в результате этого поверхность оказывается неориентируемой, а траектория движется по бутылке Клейна, если канард соединен с «внутренней» веткой Мро быстрого многообразия (левее точки Н). Если канард соединен с «внешней» веткой Мро (правее Н), то поверхность будет ориентируемой, а траектория будет лежать на торе. В случае, если параметры попадают в один из языков синхронизации (на рис. 1, Ь), динамика вновь оказывается периодической, однако траектория возвращается в исходную точку, лишь совершив несколько полных оборотов на торе или бутылке Клейна.

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

Заключение

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

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

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

Как указывалось ранее, исследование быстро-медленных систем является популярным направлением нелинейной динамики. В ходе исследования таких систем обнаруживаются разные режимы динамики, исследование которых основывается на выделении, как правило, лишь двух типичных времен — быстрого и медленного [ 18,20]. Траектории в этом случае формируются на основе единственного быстрого и медленного релаксационного цикла (канарда). Так, например, происходит в трехмерной модификации модели Розенцвейга-Макартура, описывающей динамику систем типа жертва-хищник-суперхищник [9]. Качественно такая динамика соответствует сценарию «Складка/Хопф» [17], который является наиболее простым для рассматриваемой системы (2). В медленной подсистеме (слабом сообществе) при этом наблюдаются две бифуркации: седло-узловая (типа складка) и мягкая (суперкритическая) Андронова-Хопфа, в зависимости от значения «замороженной» переменной (численность хищника сильного сообщества). Этими двумя бифуркациями начинается и заканчивается фаза пачечной динамики. Но в большинстве случаев (в большей части области S на рис. 1, a) пачечные режимы представляют собой быстрые колебания другого типа — это комбинация расходящихся и затухающих пачечных режимов. Два этих типа колебаний возникают при разных значениях дополнительной «замороженной» переменной (численность жертв сильного сообщества) и связаны с медленным многообразием в разных его частях. В этом случае расходящиеся колебания наблюдаются во время восстановления жертв на второй территории и, как правило, предшествуют восстановлению жертв первого сильного сообщества. После чего они сменяются затухающими колебаниями и два сообщества сосуществуют. В этом случае в системе, очевидно, существует более двух типичных времен. Подобный сценарий возникает, например, в некоторых трехмерных моделях нейронной активности (модель Хиндмарша-Роуза) [37,38]. Для них потребовалось выделять три разные временные шкалы, для каждой из которых была построена своя медленная и быстрая подсистема, а также описано три слоя медленных и быстрых многообразий, вложенных друг в друга [34,36,38].

В данной работе можно предположить, что затухающие колебания в рассматриваемой модели связанных сообществ — результат мягкой бифуркации Хопфа в системе с медленным временем,

а расходящиеся — жесткой (субкритической) бифуркации Хопфа в системе с супермедленным временем. Два этих цикла заканчиваются седло-узловой бифуркацией предельных циклов либо гомоклинической бифуркацией. В общем-то, описанные качественные перестройки укладываются в сценарии формирования пачечной активности, описанные Е. М. Инжикевичем [17]. Однако, из-за четырехмерной природы изучаемой модели и негладкой зависимости от «замороженной» переменной (численность жертв сильного сообщества), возникают новые комбинации этих сценариев, связанные с «невозможными» для традиционных трехмерных систем способами соединения траекторий на быстрых и медленных многообразиях (через четвертое измерение). В частности такие, которые «протыкают» быстрое многообразие и порождают траектории на бутылке Клейна, или комбинации, допускающие «пересечение» быстрых многообразий. Такая особенность динамики рассматриваемого сообщества, без сомнения, требует дальнейшего подробного изучения.

Список литературы

1. Фрисман Е.Я., Кулаков М.П., Ревуцкая О. Л., Жданова О. Л., Неверова Г. П. Основные направления и обзор современного состояния исследований динамики структурированных и взаимодействующих популяций // Компьютерные исследования и моделирование. 2019. Т. 11, № 1. С. 119-151. DOI: 10.20537/2076-7633-2019-11-1-119-151.

2. Mukhopadhyay B., Bhattacharyya R. Role of predator switching in an eco-epidemiological model with disease in the prey // Ecological Modelling. 2009. Vol. 220, no. 7. P. 931-939. DOI: 10.1016/j.ecolmodel.2009.01.016.

3. Saifuddin M., Biswas S., Samanta S., Sarkar S., Chattopadhyay /.Complex dynamics of an eco-epidemiological model with different competition coefficients and weak Allee in the predator // Chaos, Solitons & Fractals. 2016. Vol. 91. P. 270-285. DOI: 10.1016/j.chaos.2016.06.009.

4. Jansen V.A.A. The dynamics of two diffusively coupled predator-prey populations // Theoretical Population Biology. 2001. Vol. 59, no. 2. P. 119-131. DOI: 10.1006/tpbi.2000.1506.

5. Liu Y The Dynamical Behavior of a Two Patch Predator-Prey Model. Theses, Dissertations, & Master Projects. Williamsburg: College of William and Mary, 2010. 46 p.

6. Saha S., Bairagi N., Dana S. K. Chimera states in ecological network under weighted mean-field dispersal of species // Frontiers in Applied Mathematics and Statistics. 2019. Vol. 5. P. 15. DOI: 10.3389/fams.2019.00015.

7. Shen Y, Hou Z, Xin H. Transition to burst synchronization in coupled neuron networks // Physical Review E. 2008. Vol. 77, no. 3. P. 031920. DOI: 10.1103/PhysRevE.77.031920.

8. Баханова Ю.В., Казаков А. О., Коротков А. Г. Спиральный хаос в моделях типа Лотки-Вольтерры // Журнал Средневолжского математического общества. 2017. Т. 19, № 2. С. 13-24. DOI: 10.15507/2079-6900.19.201701.013-024.

9. Bakhanova Y. V., Kazakov A. O., Korotkov A. G, Levanova T.A., Osipov G. V. Spiral attractors as the root of a new type of «bursting activity» in the Rosenzweig-MacArthur model // The European Physical Journal Special Topics. 2018. Vol. 227, no. 7-9. P. 959-970. DOI: 10.1140/epjst/e2018-800025-6.

10. Huang T., Zhang H. Bifurcation, chaos and pattern formation in a space- and time-discrete predator-prey system// Chaos, Solitons & Fractals. 2016. Vol. 91. P. 92-107. DOI: 10.1016/j.chaos. 2016.05.009.

11. Banerjee M., Mukherjee N., Volpert V. Prey-predator model with a nonlocal bistable dynamics of prey // Mathematics. 2018. Vol. 6, no. 3. P. 41. DOI: 10.3390/math6030041.

12. Yao Y, Song T., Li Z. Bifurcations of a predator-prey system with cooperative hunting and Holling III functional response // Nonlinear Dynamics. 2022. Vol. 110, no. 1. P. 915-932. DOI: 10.1007/s11071-022-07653-7.

13. Smirnov D. Revealing direction of coupling between neuronal oscillators from time series: Phase

dynamics modeling versus partial directed coherence // Chaos: An Interdisciplinary Journal of Nonlinear Science. 2007. Vol. 17, no. 1. P. 013111. DOI: 10.1063/1.2430639.

14. Dahasert N., Ozturk I., Kilic R. Experimental realizations of the HR neuron model with programmable hardware and synchronization applications // Nonlinear Dynamics. 2012. Vol. 70, no. 4. P. 2343-2358. DOI: 10.1007/s11071-012-0618-5.

15. Wang L., Liu S., Zeng Y. Diversity of firing patterns in a two-compartment model neuron: Using internal time delay as an independent variable // Neural Network World. 2013. Vol. 23, no. 3. P. 243-254. DOI: 10.14311/NNW.2013.23.015.

16. Santos M.S., Protachevicz P. R., Iarosz K. C., Caldas I. L., Viana R. L., Borges F. S., Ren H.-P., Szezech Jr. J.D., Batista A. M., Grebogi C. Spike-burst chimera states in an adaptive exponential integrate-and-fire neuronal network // Chaos: An Interdisciplinary Journal of Nonlinear Science. 2019. Vol. 29, no. 4. P. 043106. DOI: 10.1063/1.5087129.

17. Izhikevich E. M.Neural excitability, spiking and bursting // International Journal of Bifurcation and Chaos. 2000. Vol. 10, no. 6. P. 1171-1266. DOI: 10.1142/S0218127400000840.

18. Han X., Jiang B., Bi Q. Symmetric bursting of focus-focus type in the controlled Lorenz system with two time scales // Physics Letters A. 2009. Vol. 373, no. 40. P. 3643-3649. DOI: 10.1016/j.physleta.2009.08.020.

19. Gu H., Pan B., Xu / Bifurcation scenarios of neural firing patterns across two separated chaotic regions as indicated by theoretical and biological experimental models // Abstract and Applied Analysis. 2013. Vol. S141. P. 374674. DOI: 10.1155/2013/374674.

20. Chen J., Li X., Hou J., Zuo D. Bursting oscillation and bifurcation mechanism in fractional-order Brusselator with two different time scales // Journal of Vibroengineering. 2017. Vol. 19, no. 2. P. 1453-1464. DOI: 10.21595/jve.2017.18109.

21. Макеева А. А., Дмитричев А. С., Некоркин В. И. Циклы-утки и торы-утки в слабонеоднородном ансамбле нейронов ФитцХью-Нагумо с возбуждающими связями // Известия вузов. ПНД. 2020. Т. 28, № 5. С. 524-546. DOI: 10.18500/0869-6632-2020-28-5-524-546.

22. 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.

23. Holling C. S. The functional response of predators to prey density and its role in mimicry and population regulation // The Memoirs of the Entomological Society of Canada. 1965. Vol. 97, no. S45. P. 5-60. DOI: 10.4039/entm9745fv.

24. Базыкин А. Д.Математическая биофизика взаимодействующих популяций. М.: Наука, 1985. 181 с.

25. Bazykin A. D. Nonlinear Dynamics of Interacting Populations. World Scientific Series on Nonlinear Science Series A: Vol. 11. New-Jersey, London, Hong Kong: World Scientific, 1998. 216 p. DOI: 10.1142/2284.

26. Rosenzweig M.L., MacArthur R. H. Graphical representation and stability conditions of predator-prey interactions // The American Naturalist. 1963. Vol. 97, no. 895. P. 209-223. DOI: 10.1086/ 282272.

27. Rinaldi S., Muratori S. Slow-fast limit cycles in predator-prey models // Ecological Modelling. 1992. Vol. 61, no. 3-4. P. 287-308. DOI: 10.1016/0304-3800(92)90023-8.

28. Кулаков М. П., Курилова Е. В., Фрисман Е. Я.Синхронизация, тоническая и пачечная динамика в модели двух сообществ «хищник-жертва», связанных миграциями хищника // Математическая биология и биоинформатика. 2019. Т. 14, № 2. С. 588-611. DOI: 10.17537/2019.14.588.

29. Курилова Е. В., Кулаков М.П., Фрисман Е.Я. Последствия синхронизации колебаний численностей в двух взаимодействующих сообществах типа «хищник-жертва» при насыщении хищника и лимитировании численности жертвы // Информатика и системы управления. 2015. № 3(45). С. 24-34.

30. Курилова Е. В., Кулаков М. П. Квазипериодические режимы динамики в модели миграционно связанных сообществ «Хищник-жертва» // Региональные проблемы. 2020. Т. 23, № 2. С. 3-11. DOI: 10.31433/2618-9593-2020-23-2-3-11.

31. Dhooge A., Govaerts W., Kuznetsov Y.A., Meijer H. G. E., Sautois B. New features of the software MatCont for bifurcation analysis of dynamical systems // Mathematical and Computer Modelling of Dynamical Systems. 2008. Vol. 14, no. 2. P. 147-175. DOI: 10.1080/13873950701742754.

32. Benoit E., Callot J. L., Diener F., Diener M.Chasse au canard // Collectanea Mathematica. 1981. Vol. 31-32. P. 37-119.

33. Арнольд В. И., Афраймович В. С., Ильяшенко Ю. С., Шильников Л. П.Теория бифуркаций // Динамические системы - 5. Итоги науки и техники. Современные проблемы математики. Фундаментальные направления. Т. 5. М.: ВИНИТИ, 1986. С. 5-218.

34. Ersoz E. K., Desroches M., Mirasso C. R., Rodrigues S. Anticipation via canards in excitable systems // Chaos: An Interdisciplinary Journal of Nonlinear Science. 2019. Vol. 29, no. 1. P. 013111. DOI: 10.1063/1.5050018.

35. Shilnikov A., Cymbalyuk G. Homoclinic bifurcations of periodic orbits on a route from tonic spiking to bursting in neuron models // Regular and Chaotic Dynamics. 2004. Vol. 9, no. 3. P. 281-297. DOI: 10.1070/RD2004v009n03ABEH000281.

36. Коломиец М. Л., Шильников А. Л. Методы качественной теории для модели Хиндмарш-Роуз // Нелинейная динамика. 2010. Т. 6, № 1. С. 23-52. DOI: 10.20537/nd1001003.

37. Hindmarsh J.L., Rose R. M. A model of neuronal bursting using three coupled first order differential equations // Proc. R. Soc. Lond. B. 1984. Vol. 221, no. 1222. P. 87-102. DOI: 10.1098/ rspb.1984.0024.

38. Linaro D., Champneys A., Desroches M., Storace M. Codimension-two homoclinic bifurcations underlying spike adding in the Hindmarsh-Rose burster // SIAM Journal on Applied Dynamical Systems. 2012. Vol. 11, no. 3. P. 939-962. DOI: 10.1137/110848931.

References

1. 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-1119-151.

2. Mukhopadhyay B, Bhattacharyya R. Role of predator switching in an eco-epidemiological model with disease in the prey. Ecological Modelling. 2009;220(7):931-939. DOI: 10.1016/j.ecolmodel. 2009.01.016.

3. Saifuddin M, Biswas S, Samanta S, Sarkar S, Chattopadhyay J. Complex dynamics of an eco-epidemiological model with different competition coefficients and weak Allee in the predator. Chaos, Solitons & Fractals. 2016;91:270-285. DOI: 10.1016/j.chaos.2016.06.009.

4. Jansen VAA. The dynamics of two diffusively coupled predator-prey populations. Theoretical Population Biology. 2001;59(2):119-131. DOI: 10.1006/tpbi.2000.1506.

5. Liu Y. The Dynamical Behavior of a Two Patch Predator-Prey Model. Theses, Dissertations, & Master Projects. Williamsburg: College of William and Mary; 2010. 46 p.

6. Saha S, Bairagi N, Dana SK. Chimera states in ecological network under weighted mean-field dispersal of species. Frontiers in Applied Mathematics and Statistics. 2019;5:15. DOI: 10.3389/ fams.2019.00015.

7. Shen Y, Hou Z, Xin H. Transition to burst synchronization in coupled neuron networks. Physical Review E. 2008;77(3):031920. DOI: 10.1103/PhysRevE.77.031920.

8. Bakhanova YV, Kazakov AO, Korotkov AG. Spiral chaos in Lotka-Volterra like models. Middle

Volga Mathematical Society Journal. 2017;19(2):13-24 (in Russian). DOI: 10.15507/20796900.19.201701.013-024.

9. Bakhanova YV, Kazakov AO, Korotkov AG, Levanova TA, Osipov GV. Spiral attractors as the root of a new type of «bursting activity» in the Rosenzweig-MacArthur model. The European Physical Journal Special Topics. 2018;227(7-9):959-970. DOI: 10.1140/epjst/e2018-800025-6.

10. Huang T, Zhang H. Bifurcation, chaos and pattern formation in a space- and time-discrete predator-prey system. Chaos, Solitons & Fractals. 2016;91:92-107. DOI: 10.1016/j.chaos.2016.05.009.

11. Banerjee M, Mukherjee N, Volpert V. Prey-predator model with a nonlocal bistable dynamics of prey. Mathematics. 2018;6(3):41. DOI: 10.3390/math6030041.

12. Yao Y, Song T, Li Z. Bifurcations of a predator-prey system with cooperative hunting and Holling III functional response. Nonlinear Dynamics. 2022;110(1):915-932. DOI: 10.1007/s11071-022-07653-7.

13. Smirnov D. Revealing direction of coupling between neuronal oscillators from time series: Phase dynamics modeling versus partial directed coherence. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2007;17(1):013111. DOI: 10.1063/1.2430639.

14. Dahasert N, Oztürk I, Kilic R. Experimental realizations of the HR neuron model with programmable hardware and synchronization applications. Nonlinear Dynamics. 2012;70(4):2343-2358. DOI: 10.1007/s11071-012-0618-5.

15. Wang L, Liu S, Zeng Y. Diversity of firing patterns in a two-compartment model neuron: Using internal time delay as an independent variable. Neural Network World. 2013;23(3):243-254. DOI: 10.14311/NNW.2013.23.015.

16. Santos MS, Protachevicz PR, Iarosz KC, Caldas IL, Viana RL, Borges FS, Ren HP, Szezech Jr JD, Batista AM, Grebogi C. Spike-burst chimera states in an adaptive exponential integrate-and-fire neuronal network. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2019;29(4):043106. DOI: 10.1063/1.5087129.

17. Izhikevich EM. Neural excitability, spiking and bursting. International Journal of Bifurcation and Chaos. 2000;10(6):1171-1266. DOI: 10.1142/S0218127400000840.

18. Han X, Jiang B, Bi Q. Symmetric bursting of focus-focus type in the controlled Lorenz system with two time scales. Physics Letters A. 2009;373(40):3643-3649. DOI: 10.1016/j.physleta.2009.08.020.

19. Gu H, Pan B, Xu J. Bifurcation scenarios of neural firing patterns across two separated chaotic regions as indicated by theoretical and biological experimental models. Abstract and Applied Analysis. 2013;S141:374674. DOI: 10.1155/2013/374674.

20. Chen J, Li X, Hou J, Zuo D. Bursting oscillation and bifurcation mechanism in fractional-order Brusselator with two different time scales. Journal of Vibroengineering. 2017;19(2):1453-1464. DOI: 10.21595/jve.2017.18109.

21. Makeeva AA, Dmitrichev AS, Nekorkin VI. Cycles-canards and torus-canards in a weakly inhomogeneous ensemble of FitzHugh-Nagumo neurons with excitatory synaptic couplings. Izvestiya VUZ. Applied Nonlinear Dynamics. 2020;28(5):524-546 (in Russian). DOI: 10.18500/ 0869-6632-2020-28-5-524-546.

22. Holling CS. Some characteristics of simple types of predation and parasitism. The Canadian Entomologist. 1959;91(7):385-398. DOI: 10.4039/Ent91385-7.

23. Holling CS. The functional response of predators to prey density and its role in mimicry and population regulation. The Memoirs of the Entomological Society of Canada. 1965;97(S45):5-60. DOI: 10.4039/entm9745fv.

24. Bazykin AD. Mathematical Biophysics of Interacting Populations. Moscow: Nauka; 1985. 181 p. (in Russian).

25. Bazykin AD. Nonlinear Dynamics of Interacting Populations. World Scientific Series on Nonlinear

Science Series A: Vol. 11. New-Jersey, London, Hong Kong: World Scientific; 1998. 216 p. DOI: 10.1142/2284.

26. Rosenzweig ML, MacArthur RH. Graphical representation and stability conditions of predator-prey interactions. The American Naturalist. 1963;97(895):209-223. DOI: 10.1086/282272.

27. Rinaldi S, Muratori S. Slow-fast limit cycles in predator-prey models. Ecological Modelling. 1992;61(3-4):287-308. DOI: 10.1016/0304-3800(92)90023-8.

28. Kulakov MP, Kurilova EV, Frisman EY. Synchronization and bursting activity in the model for two predator-prey systems coupled by predator migration. Mathematical Biology and Bioinformatics. 2019;14(2):588-611 (in Russian). DOI: 10.17537/2019.14.588.

29. Kurilova EV, Kulakov MP, Frisman EY. Effects of synchronization by fluctuations in numbers of two predator-prey communities at saturation predator growth and limitation of the victim number. Information Science and Control Systems. 2015;(3(45)):24-34 (in Russian).

30. Kurilova EV, Kulakov MP. Quasi-periodic dynamics in a model of predator-prey communities coupled by migration. Regional Problems. 2020;23(2):3-11 (in Russian). DOI: 10.31433/26189593-2020-23-2-3-11.

31. Dhooge A, Govaerts W, Kuznetsov YA, Meijer HGE, Sautois B. New features of the software MatCont for bifurcation analysis of dynamical systems. Mathematical and Computer Modelling of Dynamical Systems. 2008;14(2):147-175. DOI: 10.1080/13873950701742754.

32. Benoît E, Callot JL, Diener F, Diener M. Chasse au canard. Collectanea Mathematica. 1981; 31-32:37-119 (in French).

33. Arnold VI, Afraimovich VS, Ilyashenko YS, Shilnikov LP. Bifurcation theory and catastrophe theory. In: Dynamical Systems V. Vol. 5 of Encyclopaedia of Mathematical Sciences. Berlin: Springer-Verlag; 1994. P. 1-205.

34. Ersoz EK, Desroches M, Mirasso CR, Rodrigues S. Anticipation via canards in excitable systems. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2019;29(1):013111. DOI: 10.1063/ 1.5050018.

35. Shilnikov A, Cymbalyuk G. Homoclinic bifurcations of periodic orbits on a route from tonic spiking to bursting in neuron models. Regular and Chaotic Dynamics. 2004;9(3):281-297. DOI: 10.1070/RD2004v009n03ABEH000281.

36. Kolomiets ML, Shilnikov AL. Qualitative methods for case study of the Hindmarch-Rose model. Russian Journal of Nonlinear Dynamics. 2010;6(1):23-52 (in Russian). DOI: 10.20537/nd1001003.

37. Hindmarsh JL, Rose RM. A model of neuronal bursting using three coupled first order differential equations. Proc. R. Soc. Lond. B. 1984;221(1222):87-102. DOI: 10.1098/rspb.1984.0024.

38. Linaro D, Champneys A, Desroches M, Storace M. Codimension-two homoclinic bifurcations underlying spike adding in the Hindmarsh-Rose burster. SIAM Journal on Applied Dynamical Systems. 2012;11(3):939-962. DOI: 10.1137/110848931.

È

s^Jb

Курилова Екатерина Викторовна — родилась в г. Артем Приморского края (1983). Окончила Биробиджанский государственный педагогический институт по специальности «Математика» (2005). Работает младшим научным сотрудником лаборатории математического моделирования популяционных и экологических систем Института комплексного анализа региональных проблем ДВО РАН. Область научных интересов — популяционная биология, нелинейная динамика, синхронизация. Автор 18 работ.

Россия, 679016 Биробиджан, ул. Шолом-Алейхема, 4 Институт комплексного анализа региональных проблем ДВО РАН E-mail: [email protected] AuthorID (eLibrary.Ru): 626712

Кулаков Матвей Павлович — родился в Биробиджане Еврейской автономной области (1982). Окончил Биробиджанский государственный педагогический институт (2004). Защитил диссертацию на соискание ученой степени кандидата физико-математических наук по специальности «Биофизика» (2018). Старший научный сотрудник лаборатории математического моделирования популяционных и экологических систем Института комплексного анализа региональных проблем ДВО РАН. Научные интересы связаны с использованием методов нелинейной динамики, качественной теории дифференциальных уравнений, бифуркаций, хаоса и синхронизации для изучения особенностей функционирования пространственно распределенных живых систем: популяций животных и биологических сообществ. Автор более 50 работ.

Россия, 679016 Биробиджан, ул. Шолом-Алейхема, 4

Институт комплексного анализа региональных проблем ДВО РАН

E-mail: [email protected]

ORCID: 0000-0002-7060-2731

AuthorID (eLibrary.Ru): 170285

Фрисман Ефим Яковлевич — родился в 1948 году в Сталинабаде (ныне Душанбе). Окончил Новосибирский государственный университет (1971). Защитил диссертацию на соискание ученой степени кандидата биологических наук по специальности «Генетика» (1982) и доктора биологических наук по специальности «Биофизика» (1989). Член-корреспондент Российской академии наук по специальности «Общая биология» (2011). После окончания университета работал в Институте автоматики и процессов управления ДВО РАН: занимал должности от стажера-исследователя до заведующего лабораторией. С 2002 по 2017 — директор Института комплексного анализа региональных проблем ДВО РАН. Ныне — научный руководитель этого института. Автор более 200 публикаций по общим вопросам математического моделирования и по конкретным моделям экологических и популяционных систем. Научные интересы связаны с моделированием динамики популяционных и экологических систем, математической популяционной генетикой и задачами оптимального управления.

Россия, 679016 Биробиджан, ул. Шолом-Алейхема, 4

Институт комплексного анализа региональных проблем ДВО РАН

E-mail: [email protected]

ORCID: 0000-0003-1629-2610

AuthorID (eLibrary.Ru): 1986