Известия РАН. Серия биологическая, 2019, № 4, стр. 341-352

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

Е. Я. Фрисман 1, К. В. Шлюфман 1*, Г. П. Неверова 12

1 Институт комплексного анализа региональных проблем ДВО РАН
679016 г. Биробиджан, Еврейская автономная область, ул. Шолом-Алейхема, 4, Россия

2 Институт автоматики и процессов управления ДВО РАН
690041 г. Владивосток, ул. Радио, 5, Россия

* E-mail: shlufman@mail.ru

Поступила в редакцию 21.05.2018
После доработки 18.10.2018
Принята к публикации 20.10.2018

Полный текст (PDF)

Аннотация

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

Многие виды животных с выраженным сезонным характером процессов размножения демонстрируют четкие колебания численности популяций. Динамику таких популяций удобно описывать и исследовать с помощью достаточно простых математических моделей, отражающих большой спектр возможных динамических режимов. Для популяций лососевых видов рыб соответствующая модель – хорошо известное уравнение Рикера (Рикер, 1979). Один из базовых режимов этой модели, существующий в достаточно широкой области биологически содержательных значений популяционных параметров, – цикл с периодом два года. Такой режим динамики известен для тихоокеанской горбуши Oncorhynchus gorbuscha, у которой есть две репродуктивно изолированные линии четных и нечетных лет, порождающие двухлетний цикл с весьма различающимися значениями численностей в фазах цикла (Каев, 2002). При этом следует учесть, что горбуша – составная часть большого сообщества биологических видов, вступающих с ней в различные экологические взаимодействия (Каев, 2002; Павлов и др., 2015). Так, горбуша – основной сезонный кормовой ресурс для многих обитателей бассейнов нерестовых водоемов, включая беспозвоночных, питающихся рыбой, погибшей после нереста, и поставляющих корм для вылупившейся молоди. Кроме того, с горбушей могут активно конкурировать другие лососевые за места на нерестилищах, а для ряда речных рыб отложенные икринки – дополнительный корм.

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

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

Простейшая математическая модель, порождающая такие режимы, – модель Рикера, но с периодически меняющимся параметром. Следует отметить, что, несмотря на значительное число работ, имеющих отношение к уравнению Рикера с периодическим параметром (Ашихмина и др., 2004; Kon, 2005; Sacker, 2007; Sacker, Bremen, 2010), полного и подробного исследования возникающих и сосуществующих динамических режимов в этих моделях до сих пор не проводилось. В данной работе приводятся результаты, позволяющие заполнить некоторые из имеющихся пробелов.

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

МАТЕРИАЛЫ И МЕТОДЫ

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

Отметим, что применение простых аналитических моделей к описанию биологических объектов и их исследование – важная часть объяснительной математической теории в экологии (Розенберг, 2005). Мы хорошо понимаем ограничения модели Рикера: это модель с дискретным временем, разработанная для видов с сезонным характером размножения и унимодальной зависимостью величины пополнения от численности родительской часть популяции. Однако есть два обстоятельства, которые оправдывают использование этой модели:

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

во-вторых, выявленные закономерности могут быть обобщены и экстраполированы на другие более распространенные ситуации, поскольку модели, подобные уравнению Рикера, можно рассматривать в качестве некоторых интегралов достаточно сложных моделей с непрерывным временем, описывающих динамику биологических сообществ, состоящих из нескольких видов (Ризниченко, Рубин, 1993; Базыкин, 2003).

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

Уравнение Рикера с периодическим мальтузианским параметром, позволяющее учитывать циклические воздействия как экзогенной, так и эндогенной природы на численность популяции, – модификация модели запасопополнения (Рикер, 1979):

(1)
${{X}_{{n + 1}}} = {{X}_{n}}{\kern 1pt} {{\alpha }_{n}}\exp \left( { - \beta {\kern 1pt} {{X}_{n}}} \right),$

где переменная Xn интерпретируется как значение численности рассматриваемой популяции в момент времени n = 0, 1, …; коэффициент αn – периодический мальтузианский параметр (т.е. αn + k = αn для любого n ≥ 0, где k – период), который соответствует репродуктивному потенциалу популяции и определяется биологическими особенностями вида. Отметим, что уравнение (1) неавтономное, так как параметр αn зависит от момента времени n. Дальнейшее исследование модели (1) будет ориентировано на случай, когда период параметра αn составляет два года, т.е. k = 2.

Параметр β (β > 0) характеризует “емкость” экологической ниши популяции. Нетрудно показать, что с помощью замены (масштабирования) переменной без потери общности этот параметр может быть исключен и уравнение (1) может быть представлено в виде

(2)
${{x}_{{n + 1}}} = {{\alpha }_{n}}{{x}_{n}}{\kern 1pt} \exp \left( { - {{x}_{n}}} \right),$
где безразмерная переменная xn = βXn интерпретируется как “относительная” численность (отношение текущего значения численности популяции к уровню, равному 1/β и обеспечивающему максимальное пополнение, равное αn/(eβ)).

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

$\alpha = \frac{{{{\alpha }_{n}} + {{\alpha }_{{n + 1}}}}}{2},\,\,\,\,\rho = \frac{{{{\alpha }_{n}} - {{\alpha }_{{n + 1}}}}}{2}.$

Тогда уравнение (2) может быть записано в виде

(3)
${{x}_{{n + 1}}} = {{x}_{n}}{\kern 1pt} \left( {\alpha + {{{\left( { - 1} \right)}}^{n}}{\kern 1pt} \rho } \right)\exp \left( { - {{x}_{n}}} \right).$

Параметр α в уравнении (3) соответствует среднему значению мальтузианского параметра, относительно которого происходят колебания, а ρ – амплитуда колебаний. Так как мальтузианский параметр имеет смысл только при положительных значениях, то значение α ± ρ должно быть строго >0. Отсюда возникает естественное условие |ρ| < α.

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

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Положения равновесия и их устойчивость. Уравнение Рикера c мальтузианским параметром периода 2 – уравнение (3) – имеет тривиальное стационарное решение

(4)
$\bar {x} = 0,$

при α > 0 и |ρ| < α.

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

Исследование стационарных решений с нулевой численностью на устойчивость позволяет заключить, что условие вырождения популяции (устойчивость нулевой численности) определяется как средним значением, так и размахом (амплитудой) колебаний репродуктивного потенциала: в случае когда |ρ| = (α2 – 1)1/2 (точнее, (α2 – 1)1/2 < |ρ| < α), популяция вырождается. Отметим, что если в отсутствие колебаний репродуктивного потенциала (ρ = 0) популяция вырождается только при α < 1, то в случае колебаний возможно вырождение и при больших α, но это происходит лишь при достаточно высоких размахах колебаний и область вырождения закономерно сужается с ростом α. Вместе с тем при наличии колебаний репродуктивного потенциала популяция может существовать и не вырождаться даже в случае, когда меньшее из значений репродуктивного потенциала, реализуемое в фазе минимума кормовых запасов и равное α – |ρ|, оказывается меньше единицы. Для этого требуется, чтобы была достаточная компенсация в фазе максимума: (α + |ρ|)(α − |ρ|) > 1.

Уравнение (3) может иметь и нетривиальное стационарное решение

(5)
$\bar {x} = \ln (\alpha ),$

но оно существует только при ρ = 0. Это означает, что при постоянном репродуктивном потенциале, т.е. в случае отсутствия колебаний кормовых запасов, существует ненулевое равновесное значение численности рассматриваемой популяции. Однако появление колебаний мальтузианского параметра исключает возможность такого равновесия.

Вопрос об устойчивости решения (5) подробно рассмотрен в литературе (Рикер, 1979; Скалецкая и др., 1979; Шапиро, Луппов, 1983), где показано, что оно устойчиво при 1 < α < e2.

2-циклы и их устойчивость. В уравнении Рикера (3) с периодическим мальтузианским параметром при ρ ≠ 0 с ростом значения α при переходе через границу α2 – ρ2 = 1 на смену потерявшему устойчивость тривиальному решению (4) приходит устойчивый 2-цикл. Содержательно это означает, что двухлетние колебания кормовых ресурсов приводят к закономерному возникновению двухлетних колебаний численности зависимой популяции.

Значения численностей популяций в фазах этого цикла – элементы цикла – могут быть найдены из уравнения, где правая часть представляет собой результат 2 раза подряд итерированного уравнения (3):

(6)
$\bar {x} = {{f}^{{\left( 2 \right)}}}\left( {\bar {x},\alpha ,\rho } \right).$

Для уравнения (3) характерно согласование фазы переменной x уравнения (3) с фазой периодического мальтузианского параметра. Этому согласованию отвечает деление уравнения (6) на два подобных,

(7)
$x = \left( {{{\alpha }^{2}} - {{\rho }^{2}}} \right)x\exp \left[ { - x - \left( {\alpha - \rho } \right){\kern 1pt} x\exp \left( { - x} \right)} \right],$
(8)
$x = \left( {{{\alpha }^{2}} - {{\rho }^{2}}} \right)x\exp \left[ { - x - \left( {\alpha + \rho } \right){\kern 1pt} x\exp \left( { - x} \right)} \right].$

каждое из которых определяет только один элемент 2-цикла. Одно из них описывает четные элементы 2-цикла, а другое – нечетные.

Уравнения (7) и (8) определяют неявно заданные зависимости между x, ρ и α. Численное решение этих уравнений при фиксированном значении ρ позволяет построить график зависимости переменной x уравнения (3) от значений α. На рис. 1 представлены графические решения уравнений (7) и (8) в пространстве (α, x) для ρ = 0, 0.1 и 1.

Рис. 1.

Устойчивые решения уравнения (7) при ρ = 0 – (1), неустойчивые решения уравнения (7) при ρ = 0 – (2), устойчивые решения уравнения (7) – (3), неустойчивые решения уравнения (7) – (4), устойчивые решения уравнения (8) – (5), неустойчивые решения уравнения (8) – (6).

При ρ = 0 уравнения (7) и (8) совпадают и их график соответствует бифуркационной диаграмме классической модели Рикера (рис. 1а). В зависимости от числа решений уравнения (6) и наблюдаемого динамического режима можно выделить следующие диапазоны значений α:

при 0 < α < (ρ2 + 1)1/2 = 1 существует единственное неотрицательное решение уравнения (6) – тривиальное равновесие $\bar {x} = 0,$ которое устойчиво. Если репродуктивный потенциал <1, то популяция вырождается;

при α = (ρ2 + 1)1/2 = 1 происходит транскритическая бифуркация (Ван и др., 2005), в результате которой нулевое решение теряет устойчивость и появляется еще одно неотрицательное решение уравнения (6) – нетривиальное равновесие $\bar {x} = \ln (\alpha ),$ которое устойчиво в диапазоне 1 < α < < αbif(ρ) = e2. Другими словами, если репродуктивный потенциал больше 1, то возможно ненулевое равновесное значение численности. При постоянных внешних условиях и ненулевой начальной численности это равновесие непременно достигается, если репродуктивный потенциал не очень высок: 1 < α < αbif;

при α = αbif (ρ) = e2 нетривиальное равновесие $\bar {x} = \ln (\alpha )$ теряет устойчивость и происходит бифуркация – рождение 2-цикла. В результате при α > αbif(ρ) = e2 существуют четыре решения уравнения (6): тривиальное равновесие $\bar {x} = 0,$ нетривиальное равновесие $\bar {x} = \ln (\alpha )$ и элементы 2-цикла, которые располагаются на “вилке” по обе стороны от неустойчивого ненулевого равновесия. Таким образом, при высоком репродуктивном потенциале наблюдаются внутрипопуляционные колебания численности.

Описанная эволюция динамических режимов и числа решений уравнений (6) при ρ = 0 соответствует результатам исследований классической модели Рикера (Скалецкая и др., 1979).

Даже небольшое увеличение значения ρ – малое по амплитуде колебание репродуктивного потенциала – вносит существенные изменения в картину динамического поведения уравнения (3). Как видно на рис. 1б, 1в, при ρ > 0 из точки с координатами ((ρ2 + 1)1/2, 0) расходятся две кривые, соответствующие положительным решениям уравнений (7) и (8). Это свидетельствует о том, что в данном случае при потере устойчивости тривиального решения в положительной области фазового пространства появляется 2-цикл (в дальнейшем будем называть его “первый” 2-цикл). Другими словами, если ρ 0, тогда с ростом значения α при переходе через границу α2 – ρ2 = 1 двухлетние колебания кормовых ресурсов приводят к закономерному возникновению двухлетних колебаний численности зависимой популяции.

Дальнейшее увеличение значения α приводит к впечатляющему изменению амплитуды цикла. Вначале амплитуда небольшая и медленно растет (кривые на рис. 1б, в медленно расходятся) с ростом α. Однако при приближении α к e2 и переходе через это значение амплитуда резко возрастает. Такое поведение амплитуды позволяет заключить, что здесь возникают явления, аналогичные параметрическому резонансу, которые и приводят к формированию 2-цикла с большой амплитудой. Дело в том, что, как известно, даже в отсутствие колебаний кормовых ресурсов при значениях репродуктивного потенциала больше e2 возникают двухгодичные колебания численности, вызванные внутрипопуляционными факторами, а именно плотностно-зависимой регуляцией. Амплитуда этих колебаний быстро растет с ростом репродуктивного потенциала. Соответственно, колебания кормовых ресурсов при небольших значениях α приводят к вынужденным колебаниям численности, а при приближении и переходе через e2 вызывают раскачку внутрипопуляционных колебаний, размах которых (относительно среднего значения) может многократно превышать амплитуду колебаний кормовых ресурсов.

На рис. 1б, в видно, что уравнения (7) и (8) в диапазоне (ρ2 + 1)1/2 < α < αbif имеют в общей сложности два положительных решения и, следовательно, в этом диапазоне модель (3) имеет один 2-цикл, соответствующий единственному виду двухлетних колебаний, синхронных колебаниям репродуктивного потенциала. Однако при α > αbif (ρ) ситуация существенно меняется. Уравнения (7) и (8) имеют теперь шесть положительных решений. Два из них являются элементами синхронного первого 2-цикла: соответствующие значения переменной x лежат на крайне верхней и крайне нижней ветвях графиков. Четыре других решения соответствуют двум новым 2-циклам уравнения (3), которые рождаются в результате касательной бифуркации при α = αbif (ρ) (Кузнецов, 2001).

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

Исследования модели (3) на устойчивость и, соответственно, границы областей устойчивости тривиального решения (4) и 2-циклов (6) были проведены ранее (Шлюфман и др., 2016). Параметрический портрет уравнения Рикера с мальтузианским параметром периода 2, формируемый линиями транскритической, касательной и удвоения периода бифуркаций, представлен на рис. 2.

Рис. 2.

Параметрические области различных устойчивых решений уравнения (3). I – устойчиво тривиальное решение; II – устойчив единственный 2-цикл (первый 2-цикл); III – бистабильность: устойчивы 2 разных 2-цикла (первый и второй); IV – бистабильность: сосуществуют устойчивый второй 2-цикл и устойчивые динамические режимы, возникшие в результате бифуркаций первого 2-цикла. TC – линия транскритической бифуркации, PD – линии бифуркации удвоения периода, SN – линия касательной бифуркации; для рис. 2 и 4.

На параметрическом портрете имеется область III, в которой есть два разных устойчивых 2-цикла (рис. 2). При одних и тех же значениях популяционных параметров в зависимости от начального значения численности популяция выходит на разные циклические режимы, но имеющие один и тот же период два года. Результаты дополнительного анализа показали, что если устойчивый первый 2-цикл, возникший в результате касательной бифуркации, синхронен колебаниям мальтузианского параметра, вызванным колебаниями кормовой базы, то устойчивый второй 2‑цикл, появившийся в результате транскритической бифуркации, асинхронен этим колебаниям и поддерживается внутренними плотностно-зависимыми факторами, направленными на регуляцию роста численности. Амплитуда второго 2‑цикла меньше амплитуды первого, поскольку внешние колебания раскачивают первый цикл, но подавляют “второй”, однако при достаточно больших α (или малых ρ) полностью подавить его не могут.

Бассейны притяжения двухлетних циклов. Выясним, при каких начальных условиях (значениях численности популяции x0) популяция выходит на синхронный или асинхронный режим. Другими словами, выясним, как именно фазовое пространство (фазовая ось Ox0) разбивается на бассейны притяжения двух сосуществующих устойчивых решений – двух разных 2-циклов. Чтобы получить представление о таком делении фазовой оси исследуемой модели (3), были построены карты бассейнов притяжения устойчивых решений в плоскости (x, ρ) при фиксированных значениях α.

Построение карт бассейнов притяжения выполнялось следующим способом. При выбранном значении α с мелким шагом перебирались начальные состояния (x0, ρ0) и находились соответствующие решения уравнения (3). Каждое решение вычислялось на протяжении 5000 итераций, а по результатам последних 500 шагов выяснялось, какой именно двухлетний цикл оказался в установившемся режиме – синхронный или асинхронный. На карте бассейнов притяжения точка (x0, ρ0) окрашивалась в заданный цвет в соответствии с полученным режимом.

На рис. 3 представлена карта бассейнов притяжения при α = 0.95. Бассейны притяжения двух сосуществующих 2-циклов уравнения (3) для фиксированного значения ρ, не превышающего бифуркационное значение 0.57, представляют собой последовательность чередующихся интервалов на прямой ρ0= ρ, исходящей из точки (0, ρ0) в плоскости (x0, ρ0). Это означает, что при одной и той же амплитуде ρ колебаний мальтузианского параметра при разных начальных значениях численности популяции x0 решения уравнения (3) выходят на разные двухлетние колебания. При попадании начального значения численности в светло-серую зону эти колебания синхронны внешним воздействиям и имеют большую амплитуду, при попадании в темно-серую зону – установившиеся колебания численности асинхронны внешним воздействиям и имеют меньшую амплитуду.

Рис. 3.

Карта бассейнов притяжения уравнения (3) в фазовом пространстве (x0, ρ) при α = 9.5.

С увеличением размаха колебаний кормовых ресурсов интервалы, соответствующие областям притяжения асинхронного второго 2-цикла, сужаются, а бассейн притяжения синхронного первого 2-цикла расширяется. С дальнейшим ростом значений ρ наблюдается исчезновение областей притяжения асинхронного цикла и остается только один устойчивый синхронный цикл. Это позволяет заключить, что в большой части фазового пространства ведущий фактор – внешнее воздействие: под него подстраиваются популяционные колебания, и оно определяет их фазу.

Бифуркации удвоения периода. В классическом уравнении Рикера двухлетние колебания с ростом репродуктивного потенциала теряют устойчивость и в результате бифуркации удвоения периода появляются устойчивые четырехлетние колебания, которые с ростом значений параметра также бифурцируют по сценарию Фейгенбаума. В итоге наблюдается цепочка изменений динамических режимов 2 → 4 → 8 → 16 → 32 → … → хаос, при этом в области хаоса возникают окна периодичности (в частности, трехлетние колебания), которые также бифурцируют по сценарию удвоения периода (Скалецкая и др., 1979). В связи с этим возникает естественный вопрос: как реализуется каскад бифуркаций удвоения периода в уравнении Рикера при периодическом изменении мальтузианского параметра, с учетом того что в этой модели сосуществуют два устойчивых 2-цикла? Для анализа происходящих бифуркаций рассмотрим графическое представление элементов циклов уравнения Рикера с периодическим мальтузианским параметром в зависимости от параметра α. Элементы как устойчивых, так и неустойчивых циклов могут быть найдены из уравнения вида

(9)
$\bar {x} = {{f}^{{\left( k \right)}}}\left( {\bar {x}} \right),$
где f($\bar {x}$) – правая часть исследуемого уравнения, k – число выполненных итераций рекуррентным уравнением, $\bar {x}$ – элемент цикла. С помощью уравнения (9) можно находить элементы циклов, периоды которых одновременно и не превышают k, и делят k нацело. Для получения элементов циклов большей длины необходимо увеличение числа выполняемых итераций k.

В отсутствие колебаний кормовых ресурсов (ρ = 0) уравнение (9) соответствует бифуркационной диаграмме классической модели Рикера (Скалецкая и др., 1979). С ростом значений α усложнение динамики реализуется по сценарию Фейгенбаума. При α = αbif = e2 устойчивое нетривиальное равновесие теряет устойчивость, происходит бифуркация удвоения периода и рождается устойчивый 2-цикл. При α ≈ 12.49 2-цикл теряет устойчивость, в результате чего от каждого из элементов 2-цикла появляются еще по 2 элемента, т.е. рождается 4-цикл, и т.д. Следует отметить, что в результате каскада бифуркаций удвоения периода 4-цикл эволюционирует до 128-цикла в узком диапазоне значений α (12.49 < α < 14.79).

При наличии колебаний кормовых ресурсов (ρ > 0) потеря устойчивости нулевого равновесия в результате транскритической бифуркации сопровождается появлением устойчивого 2-цикла (первый 2-цикл). Как отмечено выше, в диапазоне (ρ2 + 1)1/2 < α < αbif это единственный 2-цикл, вызванный возбуждением колебаний периодически меняющимися внешними воздействиями. Однако при α > αbif в уравнении (3) появляются еще два 2-цикла, один из которых устойчивый (второй 2-цикл) и асинхронный колебаниям кормовых ресурсов, а другой неустойчивый (третий 2-цикл).

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

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

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

Для более детального представления о возникающих и сосуществующих динамических режимах была построена карта динамических режимов. Для этого использовался метод сканирования, позволяющий обнаруживать как глобально, так и локально асимптотически устойчивые решения в условиях мультистабильности. В диапазоне [0, ln(α + ρ)] последовательно с малым шагом перебирались начальные состояния x0уравнения (3) для каждой фиксированной пары значений параметров α и ρ в узлах прямоугольной равномерной сетки, покрывающей область {(α, ρ)|0 < α < 17; 0 < ρ < 6.2}. Далее итерационным вычислением находилось решение, из которого исключался переходный режим, вызванный начальным состоянием x0. Оставшаяся часть решения была использована для определения асимптотически устойчивого периодического решения, в бассейн притяжения которого попало выбранное начальное состояние x0. Таким образом, для каждой пары значений (α, ρ) накапливались сведения об имеющихся периодических и хаотических решениях уравнения (3), на которые оно выходит при разных начальных условиях x0. Карта динамических режимов, построенная по результатам сканирования, представлена на рис. 4. Для простоты восприятия мелкие области устойчивости циклов с периодами 16 и более (продолжающих серию бифуркаций до появления хаоса) были объединены в одну область.

Рис. 4.

Карта устойчивых динамических режимов уравнения (3) в параметрическом пространстве (α, ρ). н. д. – нерегулярная динамика; ki – область устойчивости k-цикла из i-серии бифуркации удвоения периода; ki/lj – область мультистабильности двух локально устойчивых периодических решений: k-цикла из i-серии бифуркации и l-цикла из j – серии бифуркации; точки A, B, C и E, расположенные на оси Oα – точки бифуркаций классического уравнения Рикера.

На рис. 4 в области 1 устойчиво нулевое стационарное решение. При выходе из нее (в сторону больших значений параметра α) на линии транскритической бифуркации тривиальное решение и первый 2-цикл обмениваются устойчивостью, т.е. тривиальное решение теряет устойчивость, а первый 2-цикл ее приобретает. Другими словами, при достаточном репродуктивном потенциале (и при ненулевой начальной численности) в сообществе присутствует рассматриваемый вид (тот самый, для которого независимо меняющийся вид служит кормом) и численность этого вида меняется синхронно с изменением объема пищевых ресурсов. Область устойчивости этого 2-цикла состоит из частей 21 и 21/22. В части 21 первый 2‑цикл – единственно возможное устойчивое решение – возможны только синхронные колебания. В области 21/22 первый синхронный 2-цикл делит фазовое пространство уравнения (3) со вторым асинхронным устойчивым 2-циклом, возникшим в результате касательной бифуркации. Таким образом, в области 21/22 в зависимости от начальных условий возможны два типа двухгодичных колебаний: либо синхронные с изменением запасов корма, либо асинхронные, причем амплитуды синхронных больше амплитуд асинхронных, поскольку первые резонируют с колебаниями корма, а вторые ими частично гасятся.

Дальнейшая эволюция первого 2-цикла происходит по сценарию Фейгенбаума. В результате удвоения периода синхронный 2-цикл теряет устойчивость и на смену ему приходит устойчивый 4-цикл, как мы выяснили, тоже синхронный. Область его устойчивости обозначена через 41, 41/22 и 41/42. В части 41 бассейн притяжения 4‑цикла заполняет всю область фазового пространства, соответствующую положительным значениям переменной x, а в частях 41/22 и 41/42 делит ее с устойчивыми асинхронными циклами 22 и 42 соответственно.

Отметим, что область устойчивости второго асинхронного 2-цикла проходит полосой, которая заключена между линией касательной бифуркации и линией бифуркации удвоения периода, исходящей из точки B. Данная полоса на карте состоит из частей 21/22, 41/22, 81/22 и н.д./22, где н.д. – нерегулярная динамика. В каждой из них, кроме последней, второй асинхронный 2-цикл сосуществует с одним из синхронных циклов первой серии бифуркаций удвоения периода. Отметим, что область 21/22 (устойчивости синхронного и асинхронного 2-циклов) включает в себя интервал (A, B) ≈ (e2 ≈ 7.38, 12.49) оси абсцисс Oα, на котором устойчив 2-цикл классического уравнения Рикера (Скалецкая и др., 1979).

Второй 2-цикл, как и первый, эволюционирует в хаос по сценарию Фейгенбаума. Возникающие в результате этой эволюции периодические решения в тексте упоминаются как “вторые”. Их вытянутые области устойчивости представляют собой последовательность примыкающих одна к другой полос. Накладываясь одна на другую, области устойчивости первой и второй последовательностей бифуркаций образуют различные сочетания бистабильности. Поясним это на примере второго 4-цикла. Его область устойчивости простирается вдоль области устойчивости второго 2‑цикла и примыкает к ней по линии бифуркации удвоения периода. На карте ее части обозначены через 41/42, 81/42 и н.д./42. Заметим, что, подобно 2‑циклам, области устойчивости обоих 4-циклов включают в себя интервал (B, C) (12.49, 14.24) оси абсцисс Oα, на котором устойчив 4-цикл классического уравнения Рикера (Скалецкая и др., 1979).

Анализ представленной на рис. 4 карты динамических режимов позволяет сделать определенные выводы. В рассматриваемой области параметрического пространства бистабильность имеется там, где находится область устойчивости асинхронных циклов второй последовательности бифуркации. Бистабильность проявляется в двух качественно разных сочетаниях. Первое сочетание – это два устойчивых периодических решения, одно из которых относится к первой последовательности бифуркаций, а второе – ко второй. Второе сочетание – это устойчивое периодическое решение второй последовательности бифуркаций и непериодическое решение, возникшее по сценарию Фейгенбаума в конце первой серии бифуркаций.

На рис. 5 приведены карты бассейнов притяжений устойчивых режимов при α = 12 и α = 13.5. При α = 12 (рис. 5а) карта содержит бассейны трех динамических режимов: первого 2-цикла, первого 4-цикла (результат бифуркации удвоения периода первого 2-цикла) и второго 2-цикла. При α = 13.5 (рис. 5б) число динамических режимов на карте существенно увеличивается. Здесь наблюдаются все режимы из серии бифуркаций удвоения периода первого 2-цикла, также второй 2-цикл и возникающий из него 4-цикл. Для упрощения восприятия рисунка бассейны притяжения циклов длиной >8 объединены в одну область 16-цикл и т.д.. Из рис. 5а видно, что с увеличением амплитуды колебаний мальтузианского параметра у первой серии бифуркации удвоения периода происходит увеличение периода цикла до наступления хаоса. Во “втором” 2-цикле с ростом ρ наблюдается обратный порядок бифуркации: происходит уменьшение длины цикла от 4 до 2. Как и ранее, бассейн притяжения режимов, рожденных в результате бифуркаций второго 2-цикла, с увеличением абсолютного значения ρ заметно сужается.

Рис. 5.

Карта бассейнов притяжения уравнения (3) при α =12 (а) и 13.5 (б). ki – область устойчивости k-цикла из i-серии бифуркации удвоения периода.

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

Хорошо известно, что в классической модели Рикера в зоне хаоса неоднократно появляются “окна” регулярной динамики – возникают и бифурцируют по сценарию Фейгенбаума циклы нечетной длины (3, 5, и т.д.). Достаточно давно (Скалецкая и др., 1979) численно было установлено, что 3-цикл бифурцирует в интервалах (15.985, 16.035) и (22.25, 24.5), а 5-цикл – в (18.51, 18.59). Для исследования эволюции названных циклов в модели Рикера с периодически изменяющимся мальтузианским параметром нами была построена карта динамических режимов в пространстве (α, ρ) для значений α = 0–25 и ρ = 0–25 (рис. 6).

Рис. 6.

Карта динамических режимов, демонстрирующая эволюцию 3-цикла. Числа соответствуют длинам наблюдаемых циклов, н.д. – нерегулярная динамика, 0 – тривиальное решение.

На рис. 6 в области мультистабильности параметрического пространства (α, ρ) часть полосы устойчивости 2-цикла, появившегося в результате касательной бифуркации, и возникших из него по сценарию Фейгенбаума 4-циклов, 8-циклов и т.д., вырезана, поэтому на карте отчетливо видны узкие полосы области устойчивости 6- и 12-циклов, а также 10- и 20-циклов. Эти области берут свое начало от оси Oα из уже упомянутых в предыдущем абзаце очень узких интервалов. С ростом амплитуды колебаний мальтузианского параметра ρ рассматриваемые области изгибаются, как бы повторяя линию изгиба границы первой серии бифуркации 2-цикла. Затем они утолщаются и заканчиваются резким сужением до формы “усов”, расходящихся в разные стороны. Области существования устойчивых циклов такой формы, расположенные внутри области хаоса, в литературе именуются креветкоподобными (shrimp-like) (Gallas, 1994; Gomez et al., 2014).

На карте динамических режимов имеется еще несколько достаточно крупных областей устойчивости 6-цикла. Одна из них содержит в себе интервал (22.25, 24.5) оси Oα. Две другие располагаются симметрично относительно оси Oα при больших значениях ρ.

Визуальный анализ карты позволяет заключить, что при ρ > 0 бифуркация удвоения периода первого 2-цикла происходит при меньших значениях α, чем у классического уравнения Рикера (при ρ = 0). При этом для второго 2-цикла, возникшего в результате касательной бифуркации, справедливо обратное утверждение: с ростом амплитуды этот цикл возникает при бóльших значениях, чем α = e2.

Отметим, что циклы нечетной длины возможны только при ρ = 0, а при ρ > 0 существуют только циклы четной длины (кратные периоду внешнего воздействия). Когда запасы кормовых ресурсов колеблются с периодом 2, циклы нечетной длины не наблюдаются. Вместо них в уравнении (3) возникают циклы с удвоенным периодом, которые бифурцируют по сценарию Фейгенбаума: 6, 12, …, 10, 20, …, 14, 28, …, и т.д.

Возникающие в уравнении (3) колебания имеют особенности, отличающие эти колебания от циклов классического уравнения Рикера. В частности, из всех 4-циклов, изображенных на рис. 7а–7в, только в одном случае (рис. 7а) график имеет хорошо узнаваемую пилообразную форму, свойственную 4-циклам классического уравнения Рикера. В случае (рис. 7б) решение отличается тем, что имеет локальный рост значения переменной x от меньшего к большему в интервале длиной в один период цикла, т.е. располагается на четырех значениях. В случае (рис. 7в) 4-цикл имеет два ярко выраженных локальных максимума, различающихся почти в 2 раза, а локальные минимумы расположены в узком диапазоне значений x. Другими словами, это означает, что при периодических изменениях мальтузианского параметра с большой амплитудой (ρ = 10 при α = 10.5) система может демонстрировать разные приращения переменной за одну итерацию практически из одной и той же области малых значений x.

Рис. 7.

Графики фрагментов решений уравнения (3). a – 4-цикл при ρ = 2.7 и α = 10.5, б – 4-цикл при ρ = 7.37 и α = = 10.5, в – 4-цикл при ρ = 10 и α = 10.5, г – 8-цикл при ρ = 5.68 и α = 10.5, д – 8-цикл при ρ = 6.5 и α = 12, е – 8-цикл при ρ = 8.7 и α = 10.5, ж – 10-цикл при ρ = 0 и α = 18.59, з – 10-цикл при ρ = 5.7 и α = 12. n – номер итерации.

Для 8- и 10-циклов уравнения (3) (рис. 7г–7е, 7з), как и для 4-циклов, характерен рост на протяжении трех итераций подряд. Для сравнения на рис. 7ж представлен фрагмент 10-цикла классического уравнения Рикера при α = 18.59.

ЗАКЛЮЧЕНИЕ

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

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

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

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

  1. Ашихмина Е.В., Израильский Ю.Г., Фрисман Е.Я. Динамическое поведение модели Рикера при циклическом изменении одного из параметров // Вестн. ДВО РАН. 2004. № 5. С. 19–28.

  2. Базыкин А.Д. Нелинейная динамика взаимодействующих популяций. Москва; Ижевск: Ин-т компьютерных исследований, 2003. 367 с.

  3. Ван Д., Ли Ч., Чоу Ш.-Н. Нормальные формы и бифуркации векторных полей на плоскости. Перев. с англ. / Под ред. Ильяшенко Ю.С. М.: МЦНМО, 2005. 416 с.

  4. Каев А.М. Временная структура миграционного потока горбуши Oncorhynchus gorbuscha в Охотское море // Изв. ТИНРО. 2002. № 1–3. С. 904–920.

  5. Кузнецов С.П. Динамический хаос. М.: Физматлит, 2001. 296 с.

  6. Павлов Д.С., Кириллова Е.А., Кириллов П.И., Нездолий В.К. Покатная миграция, поведение и распределение молоди рыб в низовьях реки Озерной (юго-западная Камчатка) // Изв. РАН. Сер. биол. 2015. № 1. С. 52–62.

  7. Ризниченко Г.Ю., Рубин А.Б. Математические модели биологических продукционных процессов. М.: Изд. МГУ, 1993. 301 с.

  8. Рикер У.Е. Методы оценки и интерпретации биологических показателей популяций рыб. М.: Пищ. пром-сть, 1979. 408 с.

  9. Розенберг Г.С. Модели потенциальной эффективности популяций и экологических систем // Вестн. Нижегород. ун-та. Сер. Биология. 2005. № 1. С. 163–180.

  10. Скалецкая Е.И., Фрисман Е.Я., Шапиро А.П. Дискретные модели динамики численности популяции и оптимизации промысла. М.: Наука, 1979. 168 с.

  11. Шапиро А.П., Луппов С.П. Рекуррентные уравнения в теории популяционной биологии. М.: Наука, 1983. 132 с.

  12. Шлюфман К.В., Неверова Г.П., Фрисман Е.Я. 2-циклы уравнения Рикера с периодически изменяющимся мальтузианским параметром: устойчивость и мультистабильность // Нелинейная динамика. 2016. Т. 12. № 4. С. 553–565.

  13. Шлюфман К.В., Неверова Г.П., Фрисман Е.Я. Динамические режимы модели Рикера с периодически изменяющимся мальтузианским параметром: устойчивость и мультистабильность // Нелинейная динамика. 2017. Т. 13. № 3. С. 363–380.

  14. Gallas J. Dissecting shrimps: Results for some one-dimensional physical models // Phys. A. 1994. V. 202(1–2). P. 196–223.

  15. Gomez F., Stoop R.L., Stoop R. Universal dynamical properties preclude standard clustering in a large class of biochemical data // Bioinformatics. 2014. V. 30(17). P. 2486–2493.

  16. Kon R. Attenuant cycles of population models with periodic carrying capacity // J. Diff. Eq. Appl. 2005. V. 11(4–5). P. 423–430.

  17. Sacker R.J. A note on periodic Ricker maps // J. Diff. Eq. Appl. 2007. V. 13(1). P. 89–92.

  18. Sacker R.J., Bremen H.F. A conjecture on the stability of the periodic solution of Ricker’s equation with periodic parameters // Appl. Math. Comp. 2010. V. 217. P. 1213–1219.

Дополнительные материалы отсутствуют.