CC BY
47
УДК 629.075
моделирование влияния скачкообразных
изменений в онтогенезе
на динамику системы запас-пополнение
А. Ю. Переварюха1,
аспирант
Санкт-Петербургский институт информатики и автоматизации РАН
Представлена новая нелинейная математическая модель системы запас-пополнение, реализованная в форме гибридного автомата в инструментальной среде разработки моделей на основе дифференциальных уравнений с разрывной правой частью. Дается математическая интерпретация наблюдаемому явлению деформации купола кривой воспроизводства для популяции осетра. Исследование модели показывает наличие естественных циклических колебаний численности популяции осетра в условиях отсутствия чрезмерного промысла.
Ключевые слова — гибридные динамические системы, прикладная нелинейная динамика, моделирование популяционных процессов.
Введение
Экология формулирует не строгие законы, а обобщенные принципы, называемые часто по именам ученых, впервые описавших те или иные эффекты: Олли, Либиха, «конкурентного исключения» Гаузе, «минимума видов» Ремане. Соответственно, математическое моделирование биологических систем не имеет в основе какого-то фундаментального и не оспариваемого никем положения, подобного F = Gm1m2/r2. Представьте, что при одних начальных условиях сила гравитации была бы обратно пропорциональна квадрату расстояния, а при других — только корню квадратному. В подобной лабильности — сложность математической биологии, но это же делает ее интереснейшей сферой науки, где до сих пор можно предлагать новые модели и формулировать теории.
Теория и модели формирования пополнения
Подход, основанный на роли компенсационной смертности как важнейшего регулирующего
1 Научный руководитель — доктор техн. наук, ве-
дущий научный сотрудник Санкт-Петербургского института информатики и автоматизации РАН В. В. Михайлов.
фактора, развивается в рамках теории формирования пополнения. Модели «запас-пополнение» формализуют процесс формирования численности поколения, начиная с ранних стадий развития. В предыдущей работе [1] был проведен анализ методами нелинейной динамики известной модели У. Рикера. Предложенная им в 1954 г. теория — пример того, какие проблемы могут возникать в математической биологии, ограничивающейся статистикой. В дальнейших работах Рикер развивал способы оценки коэффициентов по имеющимся данным наблюдений, не установив возможности возникновения топологически неэквивалентных фазовых портретов (в том числе образование странного аттрактора) при изменении параметров этой модели. Наш подход заключается в применении методов нелинейной динамики для изучения качественного поведения динамических моделей экологических процессов.
Автором была разработана новая модель убыли первоначальной численности поколения на интервале уязвимости £е[0, Т] в виде системы дифференциальных уравнений, учитывающая влияние уровня развития особей на смертность. Дополнительно введена функция 8(5), отражающая снижение эффективности воспроизводства при деградации популяции, которая связана с уменьшением вероятности встречи особей в сезон размножения (эффект Олли):
— = —(аю^)Щ^ + Є(Я)Р)#(0 dt
dw
dt
— = —^, к < 1, 0(S) =
s = N
Ь=Т , 4=0
1 — ехр(—cS)
= Юо, ^=0 = А£
(1)
где ю(ф) отражает изменение пищевых потребностей по мере развития особей; 5 — величина нерестового запаса; убывающая функция 85) ^ 1 и «не вмешивается» в расчеты, когда численность запаса достаточно велика; g — ограничивающий фактор, учитывающий количество доступных кормовых организмов; с — параметр, характеризующий степень выраженности эффекта Олли; X — средняя плодовитость особей; а, р — константы. Графиком предложенной и исследован-
5 ■ 103, шт.
б)
5 • 10, шт.
■ Рис. 1. Кривая модели (1) (а), Рикера (б), Бивертона—Холта (в)
ной с применением численного решения системы уравнений (1) как функции, вычисляющей численность выживших особей поколения R к моменту времени Т новой модели, является куполообразная кривая с уменьшающимся наклоном ниспадающей правой ветви (рис. 1, а). Кривая имеет две нетривиальные точки пересечения с биссектрисой координатного угла R = Б и, соответственно, две области притяжения, что определяет качественное отличие поведения траектории динамической системы от модели У. Рикера (рис. 1, б) или модели Р. Бивертона и С. Холта (рис. 1, в).
Последующее изучение деградирующих популяций позволило автору отметить некоторые интересные, но сложные и неожиданные эффекты в формировании пополнения, редко наблюдаемые у находящихся в стабильном состоянии популяций. Данные о воспроизводстве стабильной популяции будут группироваться на графике в окрестности точки равновесия. Появление куполообразной последовательности точек на графике должно служить индикатором нестабильности популяционных процессов.
Анализ характера зависимости в формировании пополнения осетра
Целью проводимых исследований является системный анализ практики хозяйственного использования биоресурсов Каспия, приведшего к деградации промысловых запасов многих популяций. Разработка моделей непосредственно возникла в результате изучения воспроизводства различных анадромных рыб и не является попыткой умозрительной популяционной интерпретации математического аппарата. Задачей работы стало дальнейшее приближение кривой модели формирования пополнения и наблюдаемого в реальности описанного ниже эффекта.
Анализ имеющихся данных [2-4] о воспроизводстве волжских популяций осетровых показал наличие очевидной, но более сложной (ср. рис. 1, а) зависимости между количеством пропущенных производителей и численностью скатившейся молоди. Визуализация данных наблюдений с использованием метода скользящей средней (которым советует пользоваться У. Рикер [5]) привела к выводу, что на некотором промежутке значений Б наблюдаются явные принципиальные отклонения от полученной в предыдущей работе теоретической формы кривой запас-пополнение, а для некоторых популяций у кривой наблюдается более одного максимума. Зависимость между нерестовым запасом и количеством скатившейся в море молоди для «осетра» характеризуется деформацией купола (рис. 2). Этот эффект свиде-
И /1КП
И : ® 400 :
И 4UU : 5 ч^о : * *
^ 300 : •
о 250 :
Ч : О 900 ' •
Я \ <д 150 *
и -2 100 : •
4 ±ии - 5 50- •
* *
О °1 1111 1111 1111 1111 . 1 1 1
О 500 1000 1500 2000 2500 3000
Пропущено производителей осетра, тыс. экз.
■ Рис. 2. Данные о воспроизводстве осетра в Волге (среднее по трем точкам)
тельствует о более сложном и распределенном во времени механизме действия компенсационной смертности. Общий характер графика является следствием длительного воздействия интенсивного промысла на многовозрастную популяцию.
Необходимо отметить, что под названием «осетр» промысловая статистика и наблюдения ЦНИОРХ суммировали без разделения данные о двух обитающих в Каспийском море видах — русском и менее многочисленном персидском осетре, фигурировавшем в отчетах часто под определением «летненерестящийся».
Проблема существования волнообразных кривых воспроизводства выявлена не автором, так называется один из подразделов работы У. Рикера «Stock and recruitment» [6], где приведена гипотетическая кривая воспроизводства горбуши реки Мак-Клинтон, но тогда такая форма кривой не получила необходимой математической интерпретации. По-видимому, эффект характерен для анадромных рыб, для которых свойственно наличие периода ската в пресной воде и солоноводного периода, в котором дополнительно действует компенсационный фактор.
Выявленный сложный характер зависимости заставил существенно переосмыслить лежащие в основе модели (1) предположения о механизмах, влияющих на скорость убыли первоначальной численности поколения, и привел к идее дифференцированного описания изменений численности в зависимости от стадии развития. Изучение биологических особенностей вида позволило сделать выводы о существенном влиянии пороговых эффектов, имеющих теоретическое биологическое обоснование в трудах В. В. Васнецова и в работе [7]. С учетом этих представлений коэффициент а нельзя считать ни постоянным, ни монотонно изменяющей значение функцией a(t) = = aw(t), и правая часть первого уравнения не является гладкой.
Гибридные системы и методы их моделирования
В реальных системах может происходить взаимодействие длительных и мгновенных процессов. При создании модели такой системы возникает необходимость приостанавливать длительные процессы и выполнять заданную последовательность дискретных действий. Еще Вито Воль-терра в 1874 г. решал задачу определения траектории полета снаряда в комбинированном гравитационном поле Земли и Луны — «ограниченную задачу трех тел» — методом деления времени на короткие интервалы. А Анри Пуанкаре [8] вводил понятие «купюры» и писал об элементах, которые необходимо будет считать принадлежащими двум разным непрерывностям. В современных условиях развития вычислительных систем целесообразно воспользоваться гибридным представлением времени:
х = {{<Gap _рге^, [0, Ti ], Gap_posti }
{Gup_ pren, \Tn-i, Tn ] GaP_ Postn }n }
где Gap_pre — «временная щель» для вычисления согласованных начальных условий и проверки предиката на левом конце промежутка очередного длительного поведения; Gap_post — «временная щель» для вычисления новых начальных условий на правом конце текущего промежутка т.; T. — время срабатывания перехода или точка, в которой становится истинным предикат события, приводящего к смене поведения. Гибридное время состоит из пронумерованной и упорядоченной последовательности «кадров», в которых непрерывная составляющая времени сменяется дискретными отсчетами. Скорость мгновенных процессов полагается бесконечной. Гибридное поведение может являться следствием совместного функционирования непрерывных объектов и дискретных регуляторов или наличия мгновенных качественных изменений в непрерывном объекте. В таких случаях до появления современных вычислительных мощностей применялся метод «припасовывания», когда нелинейные функции аппроксимировались кусочно-линейными. Иногда в моделируемой системе может происходить изменение количества составляющих ее элементов, т. е. размерность совокупного вектора состояний всей системы будет зависимой от времени величиной.
Гибридные системы выделяют в отдельный класс динамических систем в достаточной степени из-за невозможности получить в явном виде и исследовать решение многих нелинейных дифференциальных уравнений [9].
Модификация модели запас-пополнение с гибридным временем
Основным математическим аппаратом, используемым для моделирования гибридного поведения, являются дифференциальные уравнения с разрывными коэффициентами в правых частях и переменной структуры [10]. Модель записывается в среде автоматизации моделирования с набором правил для определения, какое дифференциальное уравнение с непрерывной правой частью следует интегрировать в данный момент времени и какие необходимо выбрать начальные условия.
Представим модель в виде системы с гибридным временем и кусочно-непрерывными правыми частями. Возможность реализации гибридного автомата в визуальной среде моделирования играет особую роль: наглядно вводит дискретное время, позволяющее, в том числе, изменять значение параметра а, что приводит к построению последовательности а, а1, а2. Гибридный автомат, в сущности, есть расширение идеи дискретных карт состояний ^1а1еЛагЬ>), где узлам графической формы сопоставлены процессы, описываемые дифференциальными уравнениями.
Первое уравнение из (1) для текущей численности поколения N заменяется уравнением с дважды изменяющейся правой частью:
dN
dt
-(aw(t)N(t) + 9(S)P)N(t), 0 < t < т -(a1N(x)/w(t)+P)N(t), t > т, w(t) < wk1 ,(2) -(a2w(t)N(t))N(t), wk1 < w(t) < wk
где 0 < t < т — период, определяемый биологическими особенностями вида, подробно исследованными Шмальгаузен [11]; интерпретируется как уровень развития, при достижении которого изменяется характер действия факторов смертности: смена типа питания или выход из-под пресса хищничества.
Вторая разработанная модель определяется конечным множеством состояний, с каждым из которых связана правая часть системы дифференциальных уравнений первого порядка, и множеством переходов между состояниями. Каждому переходу поставлены в соответствие условие завершения активности и функция инициализации новых начальных условий. Под состоянием в непрерывно-дискретных моделях понимаются интервалы, соответствующие системе уравнений с неизменной структурой правой части. Дополнительное основание выбора гибридного представления (2) связано с особенностью решаемых автором задач: при смешанном характере воспроизводства необходимо прерывать непрерывный расчет модельных уравнений для учета искусственного выпуска молоди в определенное время.
Примененная структура модели имеет биологическое обоснование с точки зрения теории этап-ности развития рыб, разработанной В. В. Васнецовым и его последователями (например, [12]). Суть этапного подхода в том, что на каждой стадии организм для своего развития требует иных условий. При наличии необходимых условий организм переходит на следующую стадию, причем переход происходит быстро, скачкообразно. Резкие изменения строения совершаются в очень короткий срок, в некоторых случаях меньше 3-4 ч; при этом изменяются все системы органов почти одновременно. Изменяется характер движения, изменяется способ захвата пищи в связи с изменением ее состава. Рыбы при этом меняют места своего обитания, переходя в новую по своему характеру стадию. Каждому состоянию в непрерывно-дискретной модели соответствует такой период развития молоди, в течение которого происходит только рост, но не совершается принципиальных изменений в поведении рыбы.
Два важнейших скачкообразных изменения в раннем онтогенезе анадромных рыб происходят при переходе на активное питание и при прекращении тактильного контакта с дном, после которого молодь устремляется в потоке воды к морю.
Таким образом, в разработанной модели (2) объединены представления теории этапности развития организмов и теории формирования пополнения.
Исследование разрабатываемых автором моделей рассчитано на применение современных инструментальных сред разработки имитационных моделей, в данном случае AnyLogic5. Использование вычислительных средств позволило сделать ряд важнейших открытий в нелинейной динамике. Эдвард Лоренц исследовал на ЭВМ не имеющую аналитического решения систему уравнений, что привело его к обнаружению знаменитого странного аттрактора. Теория универсальности поведения нелинейных систем сформулирована Митчеллом Фейгенбаумом в ходе исследования итерационного процесса с применением программируемого микрокалькулятора [13].
Продукт AnyLogic предлагает достаточный выбор численных методов с изменяющимся шагом интегрирования для работы с системами обыкновенных дифференциальных уравнений (но только в форме Коши) и дифференциальных уравнений с отклоняющимся аргументом. Это средство, позволяющее использовать популярный язык программирования Java, но имеющее и ряд недостатков, особенно в отношении визуализации результатов численных расчетов, может служить до некоторой степени малобюджетной альтернативой системе MATLAB. Возможностями по построению гибридных моделей обладает
и другой отечественный продукт, разработанный в Санкт-Петербурге, — MvStudium. Численное решение гибридных уравнений представляет самостоятельную задачу вычислительной математики, которую решили разработчики пакета MvStudium Ю. Б. Колесов и Ю. Б. Сениченков.
Библиотеки большинства современных пакетов моделирования содержат практически одни и те же реализации численных методов, проверенные на тестовых наборах. Применяются как одношаговые методы Эйлера, Рунге—Кутты, так и многошаговые Адамса, Куртиса—Хиршфель-дера. Профессиональные пакеты моделирования автоматически выбирают лучший метод и позволяют решать последовательные системы уравнений, сводящиеся к последовательно решаемым системам уравнений различного вида и сопрягать их решения в точках смены поведения [10].
Анализ фазового портрета
В предложенном гибком математическом аппарате итоговый вид «двоякоизогнутой» (как называл такой эффект Рикер) восходящей ветви кривой (рис. 3) зависит от степени влияния компенсационных и декомпенсационных факторов на предыдущих стадиях. Фазовое пространство динамической системы в виде полугруппы итераций {^}р0 с использованием в качестве оператора эволюции, определяющего функциональную итерацию Rj+1 = у(^) модели (2), отличается наличием двух областей притяжения и О2, границей между которыми служит репеллер — неустойчивая особая точка первого пересечения кривой с биссектрисой координатного угла. Соответственно, для траектории существует возможность притяжения к двум аттракторам. Аттрактор области — точка с координатами (0, 0) на плоскости R х Б (R = ЩТ)). Если численность популяции попадает в подмножество фазового пространства Ор то произойдет необратимая деградация популяции.
Исследование динамической системы с моделью, определяющей кривую на рис. 3, показало, что аттрактором для О2 является цикл периода 2, отвечающий режиму периодических автоколебаний: у%й*) = уа+2(Й*). На рис. 4 представлена диаграмма Ламерея (с характерным прямоугольником вокруг биссектрисы координатного угла), полученная в среде AnyLogic.
Объяснение волнообразной формы кривой с биологической точки зрения заключается в том, что высокий процент выживаемости икры в благоприятных условиях приводит к недостаточной обеспеченности молоди кормовой базой, к замедлению скорости роста, увеличению продолжительности периода уязвимости. В результате за-
5 • 102, шт.
■ Рис. 3. Кривая запас-пополнение с деформацией купола модели (2)
5 • 102, шт.
■ Рис. 4. Цикл в динамической системе с применением (2)
медляется выход поколения из-под пресса хищников. Массовый и совершенный в одном месте рыбоводными заводами выпуск молоди приводит к неестественной для среды ее концентрации, что увеличивает внутривидовую конкуренцию и плот-ностно-зависимую смертность.
Заключение
Исследование системы запас-пополнение, динамика которой определяется зависимостью в форме кривой (см. рис. 3), позволило сделать важный теоретический вывод. Популяции осетра в «неуправляемом» режиме существования свойственны естественные длиннопериодические (соотносящиеся со средним возрастом полового созревания данного вида) циклические колебания численности, обусловленные взаимозависимостью величин нерестового запаса и пополнения.
Разработанные модели позволяют объяснить следующее явление: достижение запланированных масштабов выпуска молоди осетровых восемью рыбоводными заводами в 80-е гг. (в среднем 79 млн шт. в год) не привело к ожидаемому вслед за этим увеличению вылова. По мере увеличения выпуска молоди (в 1984 г. выпуск превысил 104 млн шт.), достигнув своего максимума в объ-
еме 27,3 тыс. т (без учета вылова Исламской Республикой Иран) в 1977 г., стали сокращаться уловы осетровых. Планы стабильного получения улова в 30 тыс. т за счет выпуска 90 млн шт. 3-граммовой молоди так и остались несбыточными. В настоящее время стоимость выращивания 3-граммовой особи оценивается в 4,5 руб.
Модель (2) предсказывает эффект замедления темпов созревания и роста при масштабном искусственном воспроизводстве. Соответственно, увеличение плотности и внутривидовой конкуренции среди младших возрастных групп приводит к нарушению естественной возрастной структуры популяции в сторону непропорционального преобладания неполовозрелых особей. Коэффициенты промыслового возврата для искусственной молоди, изначально принятые как 3 %, впоследствии неоднократно пересчитывались и дискутировались. В 1989 г. они определялись для осетра 1,2 %, севрюги 1 %, белуги 0,1 %, в 1998 г. — 0,7; 0,83 и 0,07 % соответственно.
Говоря об актуальности работ в данной области, автор считает нужным отметить, что наивно ожидать от создания хороших имитационных
Литература
1. Переварюха А. Ю. Нелинейная динамическая модель системы запас-пополнение // Информационно-управляющие системы. 2008. № 2. С. 9-14.
2. Власенко А. Д. Оценка пополнения запасов волжского осетра за счет естественного воспроизводства // Осетровое хозяйство внутренних водоемов СССР: Тез. и реф. II Всесоюзного совещ. 26 февраля-2 марта 1979 г. Астрахань, 1979. С. 38-40.
3. Лагунова В. С. Воспроизводство осетра на Нижней Волге // Осетровые на рубеже 21 века: Тез. докл. Междунар. конф. 11-15 сентября 2000 г. Астрахань, 2000. С. 68-69.
4. Вещев П. В., Гутенева Г. И. Современное состояние эффективности естественного воспроизводства осетровых рыб в различных нерестовых зонах Нижней Волги // Проблемы изучения, сохранения и восстановления биологических ресурсов в XXI веке: Материалы Междунар. науч.-практ. конф. 16-18 октября 2007 г. / КаспНИРХ. Астрахань, 2007. С. 25-28.
5. Рикер У. Е. Методы оценки и интерпретации биологических показателей популяций рыб. М.: Пищевая промышленность, 1979. 408 с.
6. Ricker W. Stock and recruitment // Journal Fisheries research board of Canada. 1954. Vol. 11. N 5. P. 559623.
моделей чудодейственного рецепта по увеличению уловов или восстановлению деградировавшей популяции. Прогноз на ближайшие 30 лет по восстановлению промысловых запасов севрюги и белуги оценивается автором как неблагоприятный.
Модель ценна как обобщение опыта анализа происходивших в экосистеме процессов, но вряд ли поможет повернуть их вспять. В учебнике Р. Риклефса «The Economy of Nature» описан эксперимент над гуппи, позволивший установить, какой максимальный процент изъятия сможет выдержать аквариумная популяция. Фактически аналогичный эксперимент производился с волжскими популяциями осетровых рыб. Определенный экспертами в 1977-79 гг. процент изъятия из запаса для этих долго живущих и поздно созревающих рыб (в те годы даже без разделения по видам) в размере 60 % превышает процент изъятия, при котором можно получать максимально устойчивый улов (MSY) от быстро созревающих гуппи.
Работа выполнена при поддержке Правительства Санкт-Петербурга.
7. Еремеева Е. Ф., Смирнов А. И. Теория этапности развития и ее значение в рыбоводстве // Теоретические основы рыбоводства. М.: Наука, 1965. С. 129138.
8. Пуанкаре А. О кривых, определяемых дифференциальными уравнениями. М.: Гос. изд-во техн.-теорет. лит., 1947. 392 с.
9. Сениченков Ю. Б. Численное моделирование гибридных систем. СПб.: Изд-во Политехн. ун-та, 2004. 206 с.
10. Колесов Ю. Б., Сениченков Ю. Б. Моделирование систем. Динамические и гибридные системы. СПб.: БХВ, 2006. 224 с.
11. Шмальгаузен О. И. Закономерности предличиноч-ных и личиночных периодов развития осетровых // Биологические основы осетроводства. М.: Наука, 1983. 197 с.
12. Васнецов В. В. Этапы развития костистых рыб // Очерки по общим вопросам ихтиологии. М., 1953. С. 218-226.
13. Фейгенбаум М. Универсальность в поведении нелинейных систем // Успехи физических наук. 1983. Т. 141. Вып. 2. С. 343-374.
