Вычислительные технологии
Том 15, № 2, 2010
Исследование апериодического режима в новой гибридной модели динамики биологической популяции*
А. Ю. Переварюха Санкт-Петербургский институт информатики и автоматизации РАН, Россия
e-mail: [email protected]
Рассматривается явление смены стабильного равновесия хаотическим режимом в нелинейной динамической системе в случае, когда не происходит каскада появления топологически неэквивалентных фазовых портретов при изменении управляющих параметров. Разработанная новая непрерывно-дискретная математическая модель вида запас—пополнение основана на наличии пороговых эффектов в раннем онтогенезе анадромных рыб в соответствии с представлениями теории этап-ности развития организмов. Исследуется динамическая система, имеющая четыре нетривиальных стационарных состояния. Режим переходного хаоса реализуется при возникновении локально-несвязных границ областей притяжения регулярных аттракторов.
Ключевые слова: разработка гибридных математических моделей, моделирование динамики популяций, динамический хаос.
Введение
Прогнозирование динамики развития экологических процессов, влияющих на эффективность хозяйственного использования биоресурсов, прсдст^влябт собой весьма труд-норешаемую задачу. Сложность такого прогнозирования часто определяется недостаточным уровнем имеющихся знаний о причинно-следственных связях между основными взаимодействующими факторами в эксплуатируемых биологических системах. При поверхностном анализе проблемы результатом принятых решений становится совершенного неожиданный сценарий развития событий. Ярким примером может служить изучаемая автором ситуация с изменениями в экосистеме Каспия, являющаяся следствием многочисленных попыток "улучшения" состава видов бентоса путем интродукции новых видов или оптимизации режима промысла осетровых рыб. Итогом реализации планов, основанных на подобных прогнозах, является нс только деградация популяции осет~ ровых, приведшая к запрету их промысла, но и неожиданное резкое сокращение численности моллюсков — кормовых объектов осетровых. Вместе с тем бывают ситуации, когда невозможность адекватного прогнозирования является следствием нелинейных эффектов в динамике, связанных со сложным характером зависимости между основными величинами, определяющими развитие биологической системы.
* Работа выполнена при финансовой поддержке Правительства Санкт-Петербурга (грант № 30-04/97).
© ИВТ СО РАН, 2010.
В настоящей статье исследуются нелинейные динамические системы, имеющие биологическую интерпретацию в теории формирования пополнения популяций рыб. Данное научное направление на современном этапе развития в большей мере относится к области математической биологии и к компьютерным наукам, чем к ихтиологии, так как исследует особенности зависимости между численностью нерестового запаса и полученного от него пополнения с целью математической формализации этих величин и создания имитационных моделей экосистем. Реализация разрабатываемых автором моделей основана на применении непрерывно-дискретного представления времени и возможностей современных вычислительных средств для изучения качественного поведения траекторий гибридных систем. Методы вычислительной математики в данной области стали основным инструментом исследований и позволили сделать ряд важнейших открытий, в том числе "детерминированного непериодического потока" Э. Лоренцем [1], что стало началом непосредственного изучения явлений, теоретически ограничивающих возможность предсказания динамики развития нелинейных систем [2].
1. Математическая теория формирования пополнения популяций
Долгое время в ихтиологии считалось, что численность пополнения популяции должна находиться в линейной зависимости от численности нерестившейся в данном году части данной популяции. Впервые проблема появления немонотонной зависимости при возрастании численности нерестового стада была описана в 1954 г. канадским гидробиологом У. Рикером на основе представлений о компенсационной смертности [3]. По классификации Ф. Нива это тип смертности, влияние которой увеличивается при возрастании плотности популяции. Рикер предложил уравнение для расчета численности пополнения R — числа особей, доживших до определенного момента:
R = aS exp(-bS), (1)
где S — величина нерестового запаса; b — коэффициент, отражающий величину, обратную такому количеству выметанной икры, при котором число выжившей молоди максимально, т.е. имеет смысл b < 1; a — безразмерный параметр. На рис. 1 представлена унимодальная кривая, имеющую максимум, одну точку перегиба и единственное нетривиальное пересечение с биссектрисой координатного угла R = S.
После опубликования работы [3] в ихтиологии возникли проблемы, связанные с нелинейной динамикой. Рикер не исследовал свою модель в качестве оператора эволюции динамической системы вида Rn+p = Rn) с построением бифуркационных диаграмм. Внимание исследователей было сосредоточино на проблеме определения параметров. В частности, предлагалось строить кривую с использованием регрессии ln R/S S
ДсШН ы м наблюдений, имеет мало общего с куполообразной кривой Рикера. В.В. Суханов, применявший зависимость (1) при исследовании популяции нерки, в работе [4] писал следующее: "Оказалось, что собственно кривой Рикера результаты расчетов не дали. На плоскости S х R получилось сгущение точек, мало чем напоминающее вышеуказанную кривую. При тех параметрах, при которых популяция испытывала наиболее регулярные колебания, точки на графике образовывали нечто вроде неотчетливого эллипса. Использование метода наименьших квадратов для определения значений па-a b a
Рис. 1. График модели Рикера (1)
чем единица. Полученные результаты говорят о том. что кривая Рикера не является хорошим приближением к эмпирическим данным".
О свойствах, ограничивающих применимость модели (1). в частности о стремлении к нулю количества выжившей молоди при увеличении запаса, было отмечено во многих работах, посвященных реальному моделированию популяций. Наблюдения за воспроизводством популяций рыб [5] показывают, что с увеличением плотности запаса. во-первых, пополнение уменьшается до некоторой постоянной величины, во-вторых, аквариумная популяция стабилизируется, но не в том состоянии, при котором пополнение максимально и которое наиболее выгодно для промысла. Наложение ограничений на правую нисходящую ветвь кривой не только ведет к количественным изменениям рассчитываемых величин, но меняет качественное поведение модели. В задаче исследования эффективности воспроизводства деградирующей популяции у автора возникла другая проблема модель предсказывает нереалистично высокую выживаемость пополнения при существенном уменьшении S:
dR
—— = aexp (—6s)(1 — bS), lim aexp (—bs)(1 — bS) = a.
dS s^o
Оказалось, что моделировать воспроизводство популяции в широком диапазоне зна-S
нением модели Рикера или модели Бивертона—Холта (a,ß константы) i
dN = — (aN (t) + ß )N (t),
были признаны неудачными.
Таким образом, полученные данные определили цель дальнейших исследований расширение методологии математической теории формирования пополнения для описания наблюдаемых в реальности нелинейных эффектов.
2. Развитие динамического хаоса и бифуркации
хаос —> цикл <фп(Я) = <фп+к (R), к = 2\ для модели (1)
Описание свойств моделей запас—пополнение как графиков математических функций оказалось явно недостаточным. О том, что простые одномерные дискретные динамические системы обладают весьма разнообразным качественным поведением, впервые ЗсШВИЛ в "Nature" (1976 г.) биолог Роберт Мэй, обосновавший утверждение, что без выявления диапазона поведения, скрытого в нелинейных детерминированных разностных уравнениях, исследователь может оказаться в затруднении при осмыслении результатов моделирования на ЭВМ [6]. На работу Р. Мэя ссылались математики Дж. Иорк и Э. Ott, изучавшие кризисы странных аттракторов и метаморфозы границ областей притяжения.
Параметр a для динамической системы на основе функции (1) является управляющим. До определенного значения а, не превышающего бифуркационное, система стремится к точечному аттрактору R* = (ln a)/b. Первый метаморфоз поведения системы, представленной в виде итераций Rj+1 = aRj exp (-bRj), происходит, когда производная в неподвижной точке перестает удовлетворять следующему из теоремы Гробмана— Хартмана [7] критерию устойчивости для стационарных точек отображений:
ф'(R) = e-bR - abRe-bR,
_b Ina ln a _b Ina a(1 - ln a) n
ф'(R*) = ae b b - b——ae b b = -i—--= 1 - ln a.
r 4 ' b eln a
Поскольку критерий устойчивости \^'(R*)| < 1, то, приравнивая 1 - ln a = -1, находим
a = e2.
При a = e2 ф'(R*) = -1 и происходит бифуркация удвоения периода. Траектория динамической системы стремится к глобальному аттрактору, состоящему из двух периодических точек: фn(R*) = фn+2(R*). Эти точки являются устойчивыми стационарными точками второй итерации ф2(К). Параметр b не влияет на топологические характеристики фазового портрета.
Очевидно, что нет смысла ожидать образования куполообразной кривой от данных, полученных при наблюдениях за не находящейся в состоянии деградации популяции, воспроизводство которой описывается моделью (1). Популяция развивается как дисси-пативная динамическая система, стремящаяся к регулярному аттрактору. При аппроксимации имеющихся данных получится зависимость, функциональное преобразование на основе которой не будет обладать описанными свойствами: в работе [4] получен эллипс.
a
ния периода, универсальность которого изучена Митчеллом Фейгенбаумом на примере квадратичного отображения xn+1 = uxn(1 - xn) в работе [8]. В результате бесконечной последовательности удвоения периода при значении коэффициента a > 14.3 невозможно выделить замкнутый цикл — реализуется детерминированный хаос, напоминающий стохастический процесс. Основное свойство подобного режима движения траектории заключается в чувствительной зависимости от точности определения начальной точки траектории R0.
Траектория притягивается к странному аттрактору, отличному от конечного объединения гладких подмногообразий фазового пространства [9]. Однако общепринятого строгого определения понятия странного аттрактора не существует, так как последние исследования показали, что возможны странные, но не хаотические аттракторы.
Зависимость поведения траекторий в странных аттракторах от начальных условий ограничивает возможность использования таких моделей для прогнозирования динамики популяций. Вместе с тем известны попытки дать популяционную интерпретацию дискретным разностным уравнениям, приводящим к хаотическому режиму. Модель Тикера в режиме хаоса использовали авторы [10], предположившие, что для горбуши величина параметра а может быть порядка 60. Однако существование какой-либо популяции, динамика которой описывается в виде = аК^ ехр (—ЬЯ^) при а = 60, весьма сомнительно, так как наблюдается огромная, маловероятная для реальных популяци-онных процессов [11], амплитуда хаотических колебаний. Наличие подобной амплитуды хаотических колебаний позволило Бериману и Милштейну обосновать мнение, что популяции, которым свойственны хаотические флуктуации, предрасположены к вымиранию [12]. Известен ряд других ^^-отображений, в которых реализуется сценарий Фейгенбаума, совсем не связанных с биологической проблематикой.
Существует реальная проблема сущностной интерпретации хаоса, фрактальных странных аттракторов и других нелинейных эффектов применительно к биологии. Из-за возможности образования хаотического аттрактора преимущества использования модели Рикера выглядят весьма спорными не с математической, а с биологической точки зрения.
Спектр нелинейных эффектов (1) не ограничивается сценарием перехода к хаосу Фейгенбаума. Для несчетного числа диапазонов при а > 14.2 наблюдаются "окна периодичности", в том числе циклы нечетных периодов. Фейгенбаум рассматривал четные итерации фп(Я,а), п = 2г, то в случае, когда уравнение ф3(Я*) = Я* (при а = 22.255) будет иметь восемь решений, включая тривиальное, и, как показало численное исследование, точки Я*, Я*, Я** будут являться точками устойчивого цикла. При образовании периодического окна происходит явление, названное в [13] "субдукцией" странного аттрактора. Закрытие окна происходит после эффекта внутреннего кризиса странного аттрактора [14]. Внутри окна области притяжения странного и регулярного аттракторов пересекаются и траектория динамической системы через промежуток хаотического движения, называемого "переходным хаосом", стремится к подмножеству фазового пространства, состоящему из циклических точек с периодом 3. Как доказано Ли и Йорком [15], наличие цикла с периодом 3 означает возможность существования континуума непериодических траекторий для дискретной динамической системы. Аналогично после касательных бифуркаций возникают окна других периодов.
Если переход к хаосу через каскад бифуркаций удвоения обсуждался в биологическом контексте, то явления субдукции, кризиса аттрактора при возникновении окон с
2
логической интерпретации.
3. Разработка новой модели формирования пополнения
Для преодоления проблем, возникших с применением известных моделей при исследовании популяции волжской севрюги, автором была разработана следующая модель динамики убыли первоначальной численности поколения на "интервале уязвимости" [0, Т]:
где S — величина нерестового запаса; w(t) отражает уровень размерного развития поколения, влияющий на увеличение пищевых потребностей; g — параметр, учитывающий ограниченность количества доступных кормовых объектов; убывающая функция 0(S) ^ 1 при S ^ то и не влияет на вычисление N(T), если численность запаса достаточно велика; введение в модель быстроубывающей функции ©(S) отражает снижение эффективности воспроизводства при деградации популяции, связанное с уменьшением вероятности встречи особей в местах размножения (важный фактор, влияющий на благополучие осетровых [16]); ( — параметр, учитывающий ограничения темпов развития, не зависящие от численности; с параметр, характеризующий степень выраженности эффекта Олли; a мгновенный коэффициент компенсационной смертности; в — мгновенный коэффициент декомпенсационной смертности. Начальные условия системы уравнений определяются следующим образом: N(0) = AS, w(0) = оде А — средняя плодовитость особей популяции.
Модель основана на том, что темп развития особей находится не в обратно пропорциональной, а в обратной зависимости от численности поколения. При увеличении плотности возникает асимметричное распределение размерной структуры популяции в сторону преобладания особей с меньшими размерами. После гибели части молоди плотность снижается, но удельная биомасса по мере роста особей увеличивается. Это усиливает действие зависящих от плотности факторов.
Графиком предложенной и исследованной с применением численного решения в инструментальной среде моделирования системы уравнений (2) как функции, определя-
T
модальная кривая с уменьшающимся наклоном ниспадающей правой ветви (рис. 2, a = 0.8 • 10"14, с = 2.5 • 10"3). Кривая имеет ненулевую горизонтальную асимптоту, две нетривиальные точки пересечения с биссектрисой координатного угла R = S и характеризуется незнакопостоянным дифференциальным инвариантом Шварца. Такие свойства более соответствуют результатам наблюдений, в частности, приведенной в учебнике экологии Р. Риклефса [17] кривой воспроизводства аквариумной популяции гуппи.
Исследование всех разработанных автором моделей проводилось в инструментальной среде AnyLogic, обладающей достаточным набором численных методов с изменяющимся шагом интегрирования для работы с системами обыкновенных дифференциальных уравнений, включающих первую производную в форме Копти. Библиотеки боль-
7000
6000
н
а 5000
а 4000
о 2000 с
о
С 1000
0 2000 4000 60Ô0 8000 10000 12000 14000 16000 Нерестовый запас S, шт.
Рис. 2. Кривая запас—пополнение тта основе модели (2)
шинства пакетов моделирования содержат практически эквивалентные реализации численных методов, проверенные на тестовых наборах задач. Применяются как одноша-говые методы Рунге—Кутты. так и многошаговые Адамса и Куртиса—Хиршфельдера. Большой вклад в возможность моделирования гибридных систем внесли Ю.Б. Коле-сов и Ю.Б. Сениченков, разработавшие алгоритмы, позволяющие решать последовательные системы и сводящиеся при сопряжении их решения в точках смены поведения к последовательно решаемым системам уравнений различного вида (пакет Му5Чи(1тт).
Поведение траектории динамической системы в виде полугруппы итераций [ф^ }^>о с использованием в качестве оператора эволюции определяющего функциональную итерацию Я^+1 = ) разработанной модели (2) качественно отличается от поведения траектории систем на основе других известных моделей формирования пополнения (Рикера, Шепарда, Бивертона—Холта, Кушинга) возможностью притяжения к двум аттракторам. Соответственно фазовое пространство разделяется на две области притяжения: ^ и Границей между ними является репеллер — неустойчивая особая точка первого пересечения кривой с биссектрисой координатного угла (см. рис. 2). Анализ устойчивости неподвижных точек динамической системы, реализованной в инструментальной среде моделирования, можно проводить с использованием свойства второй итерации ф2(Я). Необходимым и достаточным условием устойчивости неподвижной точки х* одномерного отображения служат неравенства ф2(х) > х при х < х* и ф2 (х) < х при х > х* [18]. Точка Я* второго пересечения с Я^+1 = Я^ является точкой устойчивого состояния равновесия.
Траектории с начальными условиями Я0 ~ Я0, разделенными репел л ером Я*, покидают его окрестность и приближаются к разным аттракторам. Аттрактор области
— точка с координатами (0,0) на фазовой плоскости. Если численность популяции попадает в подмножество фазового пространства П1; то произойдет ее необратимая деградация.
Другое важное отличие предложенной модели от существующих состоит в том, что для динамической системы невозможен сценарий перехода к хаосу через каскад бифуркаций удвоения периода. Условием для реализации сценария Фейгенбаума в унимодальном отображении является отрицательное значение дифференциального инварианта Шварца (шварциана), но знакопостоянство шварциана невозможно при наличии более одной точки перегиба [19].
Отметим, что эффект наличия критически низкой численности Ь учитывает известная модель А.Д. Базыкина, являющаяся одной из модификаций уравнения Ферхюльс-та—Пирла. Введенный Базыкиным в правую часть логистического уравнения сомножитель (М — Ь) действует на всей области значений, которые принимает численность популяции, в том числе и тогда, когда эффект Олли проявляться не дол^кен. При численности популяции, существенно превышающей критическую, использованная функция 0(S) ~ 1.
4. Неунимодальная зависимость пополнения и запаса анадромных рыб
Результаты анализа имеющихся данных (см., например, [20]) о скатившейся в Каспийское море молоди севрюги усложнили решаемую задачу, так как показали наличие очевидной, но существенно более сложной зависимости между количеством пропущен-
700-я-
s¡ tí
« Щ a 600-
н
M R а 500-
a X 5 400 :
O) А I
л л
H в §300-
o о
к с
к о 200 -
и Ц Я :
о К 100 '
у
0 50 100 150 200 250 300 350 Пропущено производителей, тыс. шт.
Рис. 3. Запас пополнение севрюга (среднее по трем соседним точкам)
Нерестовый запас
Рис. 4. Волнообразная зависимость из работы Р. Петермапа [22J
ных производителей и численностью скатившейся молоди. Визуализация данных наблюдений с использованием метода скользящей средней по рекомендациям Рикера [11] (рис. 3) привела к выводу, что на некотором проме^кутке значении S наблюдаются явные принципиальные отклонения от полученной автором теоретической формы кривой запас пополнение (2). а для ряда популяций у кривой наблюдается более одного максимума. Возникла задача моделирования неунимодальной зависимости между запасом и пополнением.
В 2007 г. были опубликованы интересные результаты проводимых в 2002 2004 гг. Е.В. Дмитриевой [21] лабораторных исследований влияния плотности икры на смертность серой жабы Bufo bufo, показавшие, что при различной плотности икры в аквариумах имеется второй максимум выживаемости. В этой же работе исследовалось влияние плотности особей на уровень развития организмов.
О волнообразной восходящей ветви кривой для горбуши р. Карлук писал еще в пионерной работе Рикер [3]. Р. Петерман описал в [22] кривую запас пополнение (рис. 4) популяции горбуши р. Атнарко и сделал вывод, что для популяции возможны два стабильных (отметим. что стабильность C на рисунке сомнительна, так как угол наклона касательной велик) состояния: с низкой и с высокой численностью. По предположению Р. Петермана, такой вид кривой определяется различным характером смертности (компенсационным и декомпенсационньтм) на протяжении периода формирования пополнения в реке и в эстуарии. Однако из предложенной им кривой не ясно, как популяция смогла достичь стабильного состояния. Подобные гипотетические кривые, не являющиеся графиками реальных функций, в [23] приводились для объяснения вспы-тек численности попавших в новую среду видов вселенцев.
В динамической интерпретации зависимости (см. рис. 4) пересечение границы области притяжения B возможно в случае, если допускается интродукция особей ДХо из других ареалов обитания. Для динамической системы зависимость предполагает наличие двух нетривиальных аттракторов C > A > 0 с гладкой границей двух областей притяжения QA и Qc, что не позволяет моделировать наблюдаемое иногда явление!
lim én(R0) = C при R0 е QA.
5. Реализация гибридного времени в динамической системе Яп+1 = 'ф(Яп)
Найти теоретически биологическое обоснование действию факторов, приводящих к появлению сложной волнообразной зависимости запас—пополнение, оказалось возможным в рамках другой биологической теории — этапности развития рыб. С точки зрения теории этапности развития В.В. Васнецова, его последователей и в том числе согласно работам [24, 25], на каждой стадии организм для своего развития требует особых условий. При наличии этих условий он переходит на следующую стадию, причем переход происходит быстро, скачкообразно. Резкие изменения строения совершаются в очень короткий срок, в некоторых случаях в течение меньше 3-4 ч. Рыба меняет характер движения, способ захвата пищи (в связи с изменением ее состава), места своего обитания.
Два важнейших скачкообразных изменения в раннем онтогенезе анадромных рыб происходят при переходе на активное питание и прекращении тактильного контакта с дном, после которого молодь устремляется в потоке воды к морю.
Метод моделирования эффектов, связанных с этапностыо развития, был разработан с применением последовательного описания динамики убыли поколения изменяющимися уравнениями для различных стадий развития молоди. Была реализована нетрадиционная для математической биологии непрерывно-дискретная структура модели, в которой рассматриваются процессы, описываемые дифференциальными уравнениями, и события, меняющие характер протекания процессов.
Особенностью данной структуры является представление модельного времени. Время в гибридных системах рассматривается как множество последовательно расположенных отрезков, где между концом текущего интервала и началом следующего расположена "временная щель", в которой изменяются переменные состояния. Гибридное поведение может являться следствием совместного функционирования непрерывных объектов и дискретных регуляторов или наличия качественных изменений в непрерывном объекте. При достижении особых состояний в пространстве переменных состояния (событий) могут изменяться значения параметров в правых частях ОДУ, форма правой части или число уравнений. События описываются предикатами, выделяющими из всех состояний системы событие, приводящее к смене поведения. Каждому СОСТОЯНИЮ в непрерывно-дискретной модели соответствует такой период развития молоди, в течение которого происходит только ее рост, но не совершается принципиальных изменений в поведении.
В отдельный класс динамических систем гибридные системы выделяют в основном из-за невозможности получить в явном виде и исследовать решение многих нелинейных дифференциальных уравнений [26]. Чаще всего математическим аппаратом, используемым для гибридного автомата, являются дифференциальные уравнения с разрывными коэффициентами в правых частях и переменной структурой [27]. Модель определяется конечным множеством режимов изменения состояния, с каждым из которых связана правая часть системы дифференциальных уравнений первого порядка, соответствующая условиям существования и единственности решения, и множеством переходов между режимами, заданным набором булевских функций, определенных на решениях дифференциальных уравнений.
Переходам поставлены в соответствие условие завершения активности и функция инициализации новых начальных условий. Алгоритмическая реализация гибридного
автомата играет особую роль: наглядно вводит дискретное время, позволяющее в том числе изменять значение параметра а, что приводит к построению ряда значений а, а1, а2, определяющего последовательность изменения действия факторов смертности при переходах между выделенными для моделирования стадиями развития
Модель записывается в виде некоторого набора правил, по которому можно определять, какое дифференциальное уравнение с непрерывными правыми частями следует интегрировать в данный момент времени и какие необходимо выбрать начальные условия.
Первое уравнения системы (2) для описания убыли численности поколения N заменяется уравнением с переменной структурой (изменяющейся правой частью) и отклоняющимся аргументом при переходе к стадии развития Д2:
' —(а'ш(г)м(г) + в(Б)р)М(г), 0 <г < т, ^ = ] -(аМ(т)^(т) + в)М(г), г > т, ■ю(г) < ткь (3)
dt
-(a2w(t)N(t))N(t - я), wki < w(t) < wk,
где [0, т] — длительность периода pазвития D0, определяемая биологическими особенностями вида, подробно исследованными для осетровых рыб О.И. Шмальгаузен [28]; wk
действия факторов смертности и выделяет событие, приводящее к смене режима поведения гибридного автомата при Di ^ D2; я — небольшое по сравнению с временным интервалом расчета убыли численности поколения t G [0,Т] запаздывание.
Необходимость использования уравнения с отклоняющимся аргументом при описании динамики убыли численности поколения при переходе к D2 — старшей (из трех рассматриваемых) стадии развития определяется часто описываемым в экологии явлением: состояние биологической системы в данный момент зависит от процессов, происходивших в течение более или менее длительного интервала времени, предшествующего настоящему моменту. Для явления запаздывания, рассматриваемого еще в моделях В. Вольтерра [29], используются термины "последействие" и "эредитарность". Дифференциальные уравнения с запаздыванием и ранее находили применение в биологических моделях (можно отметить модификацию логистического уравнения Хатчинсоном, модель Венгерски—Каннингема и уравнение Никольсона, известное как "Blowflies equation").
6. Исследование режимов поведения модели с применением (3)
Разработанная гибридная модель представляет интерес для изучения динамического хаоса, так как это пример возникновения хаотического множества, не являющегося аттрактором в одномерном отображении.
Кривая запас—пополнение модели четырежды пересекает биссектрису Я^+1 = Я^ (рис. 5). Нетривиальными стационарными особыми точками динамической системы являются Я1 < Я2* < Я* < Я**. Использование уравнения без перехода ^ ^2 (г,г — я) сохраняет волнообразный характер восходящей ветви кривой, но у нее остаются два пересечения с биссектрисой координатного угла.
В динамике взаимно-однозначного отображения устойчивые и неустойчивые особые точки должны чередоваться. В случае разработанной гибридной модели, как показал
анализ второй итерации ф2(Я) и вычислительные эксперименты, этого не происходит. В области между неустойчивыми первой и третьей особыми точками находится неустойчивая особая точка Я**. Аттракторов, как и в динамической системе на основе (2), по-прежнему два, один из которых тривиальное положение равновесия.
3500
н ч 3000
се 2500
и
к т 2000
0)
и. Ч 1500
о
С
о с 1000
500
0
И )00 2000 3000 4000 50С
Нерестовый запас .V, шт. Рис. о. График пеуттимодалытой зависимости, полученной по (3)
Рис. 6. Переходный хаотический режим, вариант Ко € ^2
При R0 Е [R^Rg] появляются непериодические, никогда точно не повторяющиеся значения и все приблизительные повторения имеют конечную продолжительность. Траектория постепенно заполняет данную область фазового пространства, но не остается здесь навсегда, как в странном аттракторе, а с равной единице вероятностью покидает область по направлению к одному из двух возможных аттракторов. Внутри этого подмножества число непериодических повторяющихся траекторий не является бесконечным. Тип реализующегося ограниченное время хаотического поведения определяется термином "chaotic transient" переходный хаос [30]. Полученная временная диаграмма показана на рис. 6, а, фазовая для плоскости Rn х Rn+1 — на рис. 6, б.
Так как существуют R Е [R^Rg] такие, что ■ф(К) < R1, то траектория системы может устремиться в область притяжения любого из двух аттракторов (рис. 7, a временная диаграмма, б— фазовая). Области Qi и П2 не имеют четких границ, как это наблюдается в модели с применением (2).
Подобное поведение является следствием того, что ^ и П2 не образуют непрерыв-ньтх подпространств в фазовом пространстве и изменение характера границы между областями притяжения аттракторов, так же как изменение типа аттрактора при бифуркации, влияет на режим поведения траектории системы. В случае гладких границ, которыми отличались динамические системы, рассмотренные в разделах 2, 3, при опре-
Рис. 7. Переходный хаотический режим, вариант Ro Е Qi, lim фn(Ro) = 0
делении начальных условий траектория соответствовала Qi или Q2 и малое изменение начальных условий не приводило к выходу фазовой траектории на альтернативный аттрактор. В случае гибридной модели области ^ и Q2 изрезаны и имеют сложную структуру. Область притяжения аттрактора R*4 прерывается отрезками, принадлежащими области притяжения R0, и границы областей притяжения локально-несвязны ("locally disconnected basin boundaries" по классификации фрактальных границ, предложенной авторами [31]). Наличие подобного типа границ приводит к невозможности предсказания относительно того, в непрерывную часть области притяжения какого из двух существующих аттракторов попадет фазовая траектория.
В фазовом пространстве динамической системы с оператором эволюции с применением (3) можно выделить подмножество, не являющееся отдельной неустойчивой стационарной точкой или периодической траекторией, —хаотический репеллер или, как определяется в некоторых работах (например, в [32]), — область переходного хаоса ("basin of transient chaos"), где подобное явление рассматривалось на примере сложного многомерного отображения. Динамические системы с хаотическим репеллером характеризуются чувствительностью к начальным условиям, но это качественно иной вид чустви-тельности, чем наблюдающийся для траектории в странном аттракторе модели (1).
Таким образом, получено одномерное неунимодальное отображение, обладающее хаотическим режимом и характеризующееся неопределенностью асимптотической динамики. В данной динамической системе возможна бифуркация исчезновения стационарных точек RI, R4 (обратная касательная), которая приведет к тому, что Q2 становится областью притяжения возникшего хаотического аттрактора.
Т. Биологические выводы
Принципиальное с биологической точки зрения отличие зависимости, приведенной на рис. 5, от графика из статьи Петермана [22] состоит в том, что для достижения стабильного состояния при высокой численности популяции не требуется внешнего воздействия, так как существуют R G [RtRt] такие, что ^(R0) > Rt- В процессе апериодических колебаний может быть достигнуто состояние, при котором траектория динамической системы покинет диапазон, ограниченный сверху Rt, и окажется в области притяжения устойчивой точки R\. Данное состояние отражает реальную ситуацию для ряда популяций — существование такой численности, при достижении которой популяция превышает порог контроля лимитирующих факторов и происходит стремительное увеличение, "вспышка", ее численности.
В биологической интерпретации свойства динамической системы таковы, что смена устойчивого равновесного состояния хаотическими колебаниями зависит только от численности самой популяции. Если численность под воздействием промысла оказывается ниже определенной, то популяция переходит в состояние апериодических колебаний с дальнейшей возможностью как полного исчезновения, так и, в случае прекращения промыслового давления [33], — восстановления. Небольшое отклонение начальных условий приводит к изменению асимптотического поведения.
Заключение
Модели вида запас—пополнение в ихтиологической практике и, в частности, в осетровом рыбоводном хозяйстве игнорировались на основании незначительности коэффици-
ента корреляции между этими двумя величинами. В настоящее время четыре популяции (белуга, севрюга, русский и персидский осетр) потеряли промысловое значение. Промышленный лов осетровых запрещен.
Интересна история развития процесса, приведшего к деградации промысловых запасов осетровых. Может показаться странным, но данное явление деградации популяций не было вовремя спрогнозировано экспертами [34]. В 1986 г. оптимистично прогнозировалась возможность незначительного увеличения уловов волжской севрюги. Как известно, к 1992 г. уловы осетровых снизились в 2.5 раза, а за следующие восемь лет — еще в 10 раз. Экологи стали писать о деградации запасов в середине 1990-х. Авторы прогноза, непосредственно изучавшие динамику промысловых запасов осетровых Волго-Каспийского региона и участвовавшие в определении режима промыла, в работе [35] в 2000 г. писали: "Несмотря на первоначальное увеличение уловов осетровых, в дальнейшем наступили стагнация и спад численности". Таким образом, неучет описанной в настоящей работе нелинейности зависимости эффективности воспроизводства от численности не позволил экспертам вовремя принять решение о приостановке промысла.
Особую роль в динамике процесса сыграло искусственное воспроизводство. В конце 1970-х по расчетам планировалось довести выпуск молоди до 90 млн шт., что было достигнуто к 1984 г. Согласно прогнозам, такие объемы выпуска могли позволить довести вылов осетровых до 30 тыс. т. Однако существенного увеличения уловов, как отмечалось ранее, не произошло. Нелинейные модели формирования пополнения предсказывают существенное уменьшение эффективности воспроизводства при перенасыщении экологической ниши молодью одного вида. Массовый одномоментный выпуск приводит к неестественной концентрации популяции и быстрому истощению кормовой базы.
Анализ разработанных моделей показывает, что цель увеличения стабильных, получаемых длительное время уловов в два раза (как это планировалось для осетровых рыб Каспия) является нереальной. После достижения определенного значения выпуска молоди динамическая система перестает реагировать на сколь угодно большой рост этой величины. В интерпретации модели Рикера искусственный выпуск молоди может перевести траекторию системы из стационарного состояния в хаотический режим при дополнительном бифуркационном параметре, влияющем на критерий устойчивости аттрактора:
н, = 1п а/(1- АМ) , <щр= _ 1п _ АМ _
Ь аМ
и данное свойство становится еще одной проблемой при использовании модели (1).
Автор благодарит Г.И. Гутенёву (КаспНИРХ) за помощь при изучении данных наблюдений за естественным воспроизводством осетровых рыб Каспия.
Список литературы
[1] Lorenz E.N. Deterministic nonperiodic flow // J. Atmosph. Sei. 1963. Vol. 20, No. 2. P. 130-141.
[2] Короновский A.A., Трубецков Д.И. Нелинейная динамика в действии. Саратов: Изд-во ГосУНЦ "Колледж", 2002. 324 с.
Ricker W. Stock and recruitment //J. Fish. Res. Board. Can. 1954. Vol. 11, No. 5. R 559-623.
Суханов В.В. Исследование модели популяции нерки Oncorhynchus nerka в условиях изменчивой кормовой базы // Вопр. ихтиологии. 1973. Т. 13, вып. 4. С. 627-632.
Бигон Н.М. Экология. М.: Мир, 1989. Т. 1. 291 с.
May R.M. Simple mathematical model with very complicated dynamics // Nature. 1976. Vol. 261, No. 5560. P. 459-467.
Малинецкий Г.Г., Потапов А.Б. Нелинейная динамика и хаос. М.: КомКнига, 2006. 240 с.
Feigenbaum M.J. Universal behavior in nonlinear systems // Phys. D. 1983. Vol. 7, No. 1-3. P. 16-39.
Арнольд В.П., Афраймович B.C., Ильяшенко Ю.С., Шильников Л.П. Теория бифуркаций // Современные проблемы математики. Фундамент, направл. 1985. Т. 5. С. 5-220.
Животовский Л.А., Храмцов В.В. Модель динамики численности горбуши Oncorhynchus gorbusha // Вопр. ихтиологии. 1996. Т. 36, № 3. С. 369-385.
Ricker W.E. Computation and Interpretation of Biological Statistics of Fish Population. Ottawa: Supply and Services Canada, 1975. 408 p.
Berryman A., Milstein J. Are ecological systems chaotic and if not, why not? // Trends Ecol. Evol. 1989. Vol. 4. P. 26-32.
Grebogi C., Ott E., Yorke J. Crises, sudden changes in chaotic attractors and transient chaos // Phys. D. 1983. Vol. 7. P. 181-200.
Grebogi C., Ott E., Yorke J. Chaotic attractors in crisis // Phys. Rev. Lett. 1982. Vol. 48, No. 22. P. 1507-1510.
Li Т., Yorke J. Period three implies chaos // Amer. Math. Monthly. 1975. Vol. 82, No. 10. P. 985-990.
Переварюха А.Ю. Анализ воспроизводства популяций рыб на основе динамических систем // Экологические системы и приборы. 2008. № 1. С. 40-44.
Ricklefs R. The Economy of Nature. Portland: Chiron Press, 1976. 302 p.
Магницкий H.A., Сидоров С.В. Новые методы хаотической динамики. М.: Едиториал УРСС, 2004. 320 с.
Шарковский А.Н., Коляда С.Ф., Сивек А.Г. Динамика одномерных отображений. Киев: Наук, думка, 1989. 216 с.
Переварюха Т.Ю. Современное состояние и некоторые вопросы восстановления биоразнообразия осетровых рыб // Вестник Астраханского гос. техн. ун-та. 2008. № 3. С. 33-38.
Дмитриева Е.В. Влияние плотности икры на темпы развития и смертность серой жабы (Bufo bufo) в лабораторных условиях // Зоолог, журн. 2007. Т. 86, № 2. С. 229-235.
Peterman R.M. A simple mechanism that causes collapsing stability regions in exploited salmonid population //J. Fish. Res. Board Can. 1977. Vol. 34. P. 1130-1142.
Свирежев Ю.М., Логофет Д.О. Устойчивость биологических сообществ. М.: Наука, 1978. 352 с.
Васнецов В.В. Этапы развития костистых рыб // Очерки по общим вопросам ихтиологии. М.: Изд-во АН СССР, 1953. С. 207-217.
[25] Еремеева Е.Ф., Смирнов А.И. Теория этапности развития и ее значение в рыбоводстве // Теоретические основы рыбоводства. М.: Наука, 1965. С. 129-138.
[26] Сениченков Ю.Б. Численное моделирование гибридных систем. СПб.: Изд-во Политехи, ун-та, 2004. 206 с.
[27] Колесов Ю.Б., Сениченков Ю.Б. Моделирование систем. Динамические и гибридные системы. СПб.: БХВ, 2006. 224 с.
[28] Шмальгаузен О.И. Закономерности предличиночных и личиночных периодов развития осетровых // Биологические основы осетроводства. М.: Наука, 1983. С. 196-199.
[29] Volterra V. Theory of Functionals and of Integral and Integro-Differential Equation. N.Y.: Dover PubL, 1959. 304 p.
[30] Grebogi C., Ott E., Yorke J. Chaos, strange attractors and fractal basin boundaries in nonlinear dynamics // Science. 1987. Vol. 238, No. 4827. P. 632-638.
[31] MacDonald S., Grebogi C., Ott E., Yorke J. Fractal basin boundaries // Phys. D. 1985. Vol. 17, No. 2. P. 125-153.
[32] de Souza S.L., Caldas L., Viana R.L. Sudden changes in chaotic attractors and transient basins in a model for rattling in gearboxes // Chaos, Solitons and Fractals. 2004. Vol. 21. P. 763-772.
[33] Christie W.J. Changes in the fish species composition of the Great Lakes //J. Fish. Res. Board Can. 1974. Vol. 31, No. 5. P. 827-853.
[34] Власенко А.Д. Естественное воспроизводство осетра Волго-Каспия в современных условиях и его роль в формировании запаса этого вида // Динамика численности промысловых рыб. М.: Наука, 1986. С. 177-189.
[35] Ходоревская Р.П., Довгопол Г.Ф.. Журавлёва О.Л. Динамика промысловых запасов осетровых Волго-Каспийского региона // Вопр. рыболовства. 2000. Т. 1, вып. 2-3. С. 160-162.
Поступила в редакцию 17 февраля 2009 г., с доработки — 3 декабря 2009 г.