О взаимодействии имитационных и аналитических методов при моделировании эколого-биологических объектов Саранча Д.А. Вычислительный центр им. А.А.Дородницына РАН, Москва, Россия Тращеев Р.В. Институт фундаментальных проблем биологии РАН Г. Пущино

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

Первые опыты (аналитического) моделирования, заложившие основы количественной экологии, показали свою эффективность, и вместе с тем выявили важнейшую роль междисциплинарных взаимодействий, их определяющее значение при оценке адекватности модели. В аналитический период возможности исследователей были ограничены. «Мы максимально отвлекаемся от конкретных деталей моделируемого явления…Наиболее общие модели оказываются и наиболее простыми» … в силу этого допускают строго логическое математическое исследование». [стр. 13. Фомин С.В., Беркенблит М.Б. Математические проблемы в биологии. М. Наука с.].

Ситуация изменилась после создания Дж.Форрестером модели «Мировой динамики», с привнесенным им в экологические исследования культуры использования ЭВМ – метода создания имитационных моделей в диалоге с экспертами. Такой подход позволяет принимать к рассмотрению практически все предложения экспертов – в количественной или качественной форме, относительная простота модификации моделей такого типа позволяет проводить сравнительный анализ различных наборов исходных предположений, данных, гипотез. Моделирование пошло по пути «конкретизации». Использование вычислительного эксперимента позволило снять ограничения на степень детализации модели и на тип используемых уравнений, решать многопараметрические задачи, создало возможность изучения динамических режимов в многомерном пространстве параметров, позволило повысить требования к обоснованию уравнений.

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

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

Сбор, отбор, анализ и переработка исходной (биологической) информации. Обоснование и построение детальных имитационных моделей; анализ их свойств. Формирование имитационной системы – набора взаимосвязанных моделей разной степени детализации; включающего в себя упрощенные модели, допускающие проведение аналитического (портретного) исследования. Формулирования на их основе гипотез о ведущих механизмах исследуемого явления. Создание упрощенных (аналитических) моделей осуществляется посредством совместного анализа эколого-биологической информации и результатов вычислительных экспериментов на основании редукции базовых имитационных моделей.

Метод КОИС сформировался в процессе моделирования тундрового сообщества. динамики колебаний численности животных тундры. Был создан набор взаимосвязанных моделей – имитационная система, включающая разнообразные модели, как по экологическим объектам (тундровые популяция и сообщества), так и по математическому типу моделей.

Основа набора - модель тундрового сообщества «растительность- лемминги-песцы» (РЛП). Она была создана совместно с биологами, с привлечением экспертных оценок и с учетом сезонного изменения параметров. (Совместно с В.А. Орловым, Е.В. Байбиковым, Л.М. Шиляевой, И.В. Дмитриевой, Н.В.Белотеловым, В.Н. Глушковым и др.)

ЛЕММИНГ

Выбор объекта моделирования и структуры его математического описания осуществляется как компромисс между математическими и экологическими требованиями, на основе 1) принципа минимальности - использование минимально возможной математической структуры, необходимой для имитации изучаемого явления, с учетом два других – 2) системности - учет многообразия связей внутри и вне изучаемого объекта и 3) соответствия (экологичности) - привлекать предположения, не противоречащие доступным экологическим данным.

Модель РЛП была создана на основе экспертно оцененных зависимостей, учитывающих сезонные изменения параметров При этом были использованы междисциплинарные возможности компьютерных технологий, идея «экологического конструктора» при первоначальном формировании модели и ее последующих модификациях. Реализация этой идеи основана на объединении системной динамики Дж. Форрестера с гипотезой Вольтерра - Костицына – о возможности использования для описания экологических объектов систем обыкновенных дифференциальных уравнений. Наличие ЭК позволяет сравнительно просто (на формальном уровне) рассматривать альтернативные подходы. При этом реализуется тезис об устранении зависимости результатов моделирования от конкретной параметризации, восходящий к идее А.Н.Колмогорова.

Экологический конструктор На основании гипотезы Вольтерра - Костицина о возможности представления в модели взаимодействия популяций с помощью систем обыкновенных дифференциональных уравнений первого порядка, в которых одной переменной описывается биомасса (численность) популяции или трофического уровня модель выглядит следующим образом: dA/dt = R a - D a - M a (1) Здесь величины R a,D a, M a описывают скорости репродукции, отчуждения и естественного отмирания, а величина A является показателем одного из трофических уровней: W, L, F (растительность, лемминги, песцы). Каждая из компонент, в свою очередь, формируется произведением константы и на ряд сомножителей (соответствующих безразмерных функций), каждый из которых уже может обладать достаточно прозрачными свойствами и вид этих зависимостей может быть получен в результате диалога со специалистами соответствующей области.

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

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

При формировании базовой модели было создано достаточно много ее версий. Модификации моделей состояли в изменении констант, функций или в выборе нового набора функциональных зависимостей. Наиболее кардинальные изменения модели были связаны со сменой следующих ведущих предположений: взаимодействия по типу "хищник- жертва", пороговая зависимость скорости размножения леммингов от биомассы корма, определяющее влияние внутрипопуляционной динамики леммингов - принцип Олли: наличие плотности леммингов, оптимальной для размножения (Саранча, 1996).

Фазовый портрет модели РЛП Где L- биомасса леммингов, V- биомасса растительности. Жирная линия – одна из реализованных траекторий. Тонкие линии – фазовые траектории в различные сезоны.

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

Неудовлетворенность традиционным окончанием исследований имитационных моделей – прогнозом динамических характеристик модели при различных сценариях внешних воздействий, и стремление приблизиться к пониманию механизмов формирования динамики численности тундровых животных привело к формированию завершающих этапов метода КОИС. Была создана модель популяции леммингов, определяющих характер колебаний численности животных тундрового сообщества, что привело к обоснованию возможности использования в качестве упрощенной модели одномерного разностного уравнения, связывающего численности леммингов (ведущего блока в модели РЛП) в двух соседних годах.

Модель популяции леммингов Была построена и проанализирована модель популяции леммингов с учетом возрастной структуры c привязкой к конкретному региону: пос. Тарея (Западный Таймыр). В основу модели были положены данные В.А. Орлова (1985), дополненные по литературе, и экспертными оценками. Рассматривались два вида леммингов, обитающих в данном регионе – копытных и сибирских.

Проведенные на ЭВМ расчеты воспроизвели характерные для Западного Таймыра трехлетние колебания численности для сибирских и копытных леммингов (см. Орлов и др., 1986, Саранча, 1996). В результате аппроксимации результатов вычислительных экспериментов были получены две функции последования : 1, 1 – зависимость плотности в конце сезона размножения от плотности перезимовавших соответственно копытного и сибирского лемминга; 2 – зависимость плотности в начале сезона размножения от плотности в конце предыдущего сезона (выраженные в процентах от максимальной численности копытных леммингов). Пунктиром выделен 3-х летний цикл.

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

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

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

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

Анализ дискретных уравнений создал возможность выделения трех главных показателей, которые формируют динамические режимы колебаний численности леммингов. Такими показателями являются скорость прироста биомассы в благоприятный год; максимальная численность и выживаемость в наиболее неблагоприятных условиях (или же два безразмерных показателя - относительная скорость прироста численности популяции и доля гарантированно выживших зверьков).

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

Упрощенные модели, гипотезы.., Но удалось реализовать еще один аспект, который существенно повышает эффективность КОИС. Аналогов в доступной литературе обнаружить не удалось.

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

Это позволило нелинейная, неавтономная модель свести к набору систем из двух автономных линейных дифференциальных уравнений с постоянными коэффициентами:.

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

Описанный выше набор моделей приводит к созданию «имитационной системы». В рамках имитационной системы удается провести согласование моделей разного класса – моделей «растительность – лемминги – песцы» и соответствующего разностного уравнения. 1 уровень – детальная имитационная модель 2 уровень – упрощенная аналитическая модель 3 уровень – разностное уравнение

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

Что такое адекватность в период перехода описательных наук в точные? Затруднения в оценке адекватности математических моделей связаны с тем, что, как правило, теоретические построения опережают экспериментальные - порой даже не налажена процедура фиксирования соответствующих показателей. Полезность модели определяет, прежде всего, практика – доверие к ней биологов, возможность ее использования для планирования экспериментально-полевых исследований, а так же для теоретических построений. В этих условиях степень адекватности модели фактически заменяется степенью доверия биологов к качеству, как самой модели, так и процедуры ее создания. Методика КОИС дает пример такой процедуры

В результате применения КОИС может быть преодолено традиционное недоверие биологов к аналитическим моделям, с помощью обоснования их посредством процесса редукции базовых имитационных моделей, тех моделей через которые осуществляется непосредственное взаимодействие с экспертами-биологами, с исходным фактологическим и концептуальным материалом. Тем самым претензии к необоснованности исходных предположений модели переходят в плоскость доверия к экспертам и исходным биологическим данным, к возможности математических средств их представления. А возможности эти, с применением ЭВМ, весьма значительные.

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

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

Имитационные модели, нацеленные на поиск механизмов явления, поставляют и обосновывают новые классы задач. 1.Исследование нового типа дискретных уравнений. 2.Модели, учитывающие «сезонность». Потребность в расширении набора взаимосвязанных моделей приводит к использованию индивидуально- ориентированных моделей

Вид дискретных уравнений, связывающих численности леммингов в двух соседних годах, полученных в результате анализа моделей тундровых популяций и сообществ. Где P - скорость прироста биомассы леммингов в благоприятный год, X – равновесная численность, d – численность леммингов в оптимальном биотопе.

Зависимость траекторий расчётной модели от параметра d. Результаты вычислительных экспериментов.

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

Зависимость среднего расстояния между пиками численности от уровня оптимального биотопа

Тем самым удалось получить результат, предсказанный с помощью упрощенной модели

Модели с сезонностью. Совместно с И.Н. Юферовой, А.И. Лобановым, С.П. Поповым, Ю.И. Бибик и др.

Рассмотрена система: – в летний период – в зимний период где X, Y – биомассы жертв и хищников, a, d – коэффициенты, описывающие прирост биомассы жертв и смертности хищников, b, c – коэффициенты, описывающие взаимодействия видов, - коэффициент «самолимитирования» жертв, a 1, d 1 – коэффициенты вымирания. Параметры модели выбирались следующим образом: a = s, b = s, c = 0.01b, d = 0.01s, a 1 = 0.1, d 1 = 0.01 Параметр s характеризует скорость восстановления растительности.

«Кольца Сатурна» Колебания численности жертв и хищников происходят все время на одних кривых, около оптимальных траекторий. далее

Изменение биомассы жертв Изменение биомассы хищников Фазовый портрет и динамика модели назаддалее «кольца Сатурна»

Фазовый портрет модели Параметры модели: Величина параметра, описывающего зимнее вымирание жертв, сопоставима с величиной, описывающей их прирост в летний период

Динамика модели Изменение биомассы жертв Изменение биомассы хищников

Динамика развития популяций. U,V – плотность биомасс жертв и хищников. Шаг по времени – 5 с.

Динамика развития популяций. U,V – плотность биомасс жертв и хищников. Шаг по времени -1с.

Индивидуально- ориентированные модели. Исследование популяции леммингов. Задача распространения заболеваний. Создание универсальных подходов (разные виды, разные биотопы и т.д.). (Совместно с В.Д. Перминовым, Э.В. Недоступовым, А.В. Фроловой и др.)

Индивидуально-ориентированные модели.

Динамика численности леммингов

Распределение особей по возрасту и по потенциалу жизнестойкости

Заселённость территории в период депрессии и пика

Способы исследования свойств дискретных отображений. Порядок натурального ряда. Линия возврата. Отображение за положение равновесия. Совместно с Н.В. Белотеловым, И.В. Дмитриевой, В.Н. Глушковым, Э.В. Недоступовым, Ю.С. Юрезанской и др.

Расчет на ЭВМ модификации модельной задачи, если X t

Зависимость траекторий расчётной модели от параметра d. Результаты вычислительных экспериментов.

Для характеристики транзитной области [A, 1], в которой траектория не может находиться два такта подряд и для описания ее взаимодействия с областью длительного пребывания [0, A] введем две конструктивные конструкции. Будем использовать отображение области [А, 1] на себя - отображение за «положение равновесия» - отображение (ПР), а также метод построения «линий возврата» ЛВ.

Отображение ПР определяется как отображение отрезка [A, 1] на себя при котором каждому значению X из этого отрезка ставится в соответствии значение Y, равное значению X t+1 (X 0 =X) при первом возвращении траектории за положение равновесия, т.е. все X t A. ОПР позволяет изучать свойства исходного отображения с несколько иного ракурса. Метод особенно удобен для анализа разрывных функций

Новое отображение можно исследовать обычными методами (поиск стационарных точек, n-кратное отображение и т.д). Соответствующим линейным преобразованием отрезок [A, 1] можно привести к отрезку [0, 1].

Определение. ЛВ n-го порядка (ЛВn) для отображения F называется кривая в прямоугольнике A

Алгоритм построения ЛВn. Проведем горизонтальную линию для X t+1 из отрезка 0< X t+1

Линией ЛВ n-го возврата для отображения F называется кривая F c, определенная в области N t >A, обладающая следующим свойством: если начальное значение N o лежит на кривой F c, то траектория, из него исходящая, при n-ом возврате в область N t >A проходит через ту же точку N o. ЛВ нужна для определения периода цикла.

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

Исходное отображение (0DAD 2 ),его отображение за ПР (рис. справа) и их линии первого возврата (AD 1 A 1 D 2 A 2 )

Треугольное n-кратное отображение

Треугольное отображение и его линии возврата и отображение за ПР Линии возврата и отображение за ПР построены при помощи n-кратных отображений и их поворотов на 90 град.

Треугольное отображение, дополненное ступенькой справа Зависимость траекторий расчётной модели от высоты ступеньки d (0

Треугольное отображение Зависимость длины цикла от высоты ступеньки d (0

Треугольное отображение и его линии возврата

Симметрии Треугольное отображение и его повороты на 90, 180, 270 град. Линия первого возврата Любая точка пересечения лежит на периодической траектории

Симметрии Треугольное 1, 2, 3-кратное отображение и его повороты на 90, 180, 270 град Любая точка пересечения лежит на периодической траектории

Двухзонное дискретное отображение (ДДО) и его линия 1-го возврата

1, 2, 3, 4, 5-кратное отображение ДДО

Зависимость траекторий расчётной модели от высоты ступеньки d (0

Отображение ДДО Зависимость длины цикла от высоты ступеньки d (0

2 первых линии возврата для ДДО и порождаемые ими периодические точки Красные точки на биссектрисе – периодические. Точки на горизонтальных линиях имеют тот же период, что и точка на биссектриссе

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

Первые шесть линий возврата для ДДО Симметрии ЦИФРАМИ ОБОЗНАЧЕН ПОРЯДОК ЛИНИИ ВОЗВРАТА. ТОЧКА ПЕРЕСЕЧЕНИЯ ЛИНИЙ ДРУГ С ДРУГОМ ПОРОЖДАЕТ ПЕРИОДИЧЕСКИЕ ТОЧКИ.

Линии возврата с 7 по 14 ЦИФРАМИ ОБОЗНАЧЕН ПОРЯДОК ЛИНИИ ВОЗВРАТА. ТОЧКА ПЕРЕСЕЧЕНИЯ ЛИНИЙ ДРУГ С ДРУГОМ ПОРОЖДАЕТ ПЕРИОДИЧЕСКИЕ ТОЧКИ.

Координаты линий возврата (с 1 по 14) при трёх близких, но разных исходных значениях

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

б)

Логистическая функция, дополненная прямой близкой к вертикали после точки равновесия.

Результаты вычислительного эксперимента с логистическим отображением, дополненным ступенькой при изменении высоты ступеньки d (0

Зависимость периода цикла от высоты ступеньки d

Монографии Саранча Д.А. Биомоделирование. М.: ВЦ РАН, с. Саранча Д.А. Биомоделирование. Материалы по количественной экологии. Математическое моделирование и биофизические аспекты. М.: ВЦ РАН, с. Саранча Д.А. Количественные методы в экологии. Биофизические аспекты и математическое моделирование. М., МФТИ, с. Бибик Ю.В., Попов С.П., Саранча Д.А. Неавтономные математические модели экологических систем. М.: ВЦ РАН, с. Саранча Д..А., Сорокин П.А., Фролова А.А. Математическое моделирование динамики численности популяций животных. М.: ВЦ РАН, c. Глушков В.Н., Недоступов Э.В., Саранча Д.А, Юферова И.В. Компьютерные методы анализа математических моделей экологических систем. М.: ВЦ РАН, c.

Статьи Байбиков Е.В., Белотелов Н.В., Завьялова С.В., Обридко И.В., Орлов В.А., Саранча Д.А., Шелепова О.В., Шиляева Л.М. О моделировании тундровых популяций и сообществ. Математическое моделирование. Процессы в сложных экономических и экологических системах. М.: Наука, С Орлов В.А., Саранча Д.А., Шелепова О.А. Математическая модель динамики численности популяции леммингов (Lemmus, Dicrostonyx) и ее использование для описания популяций Восточного Таймыра/Экология. 1986, N 2. С Лобанов А.И., Саранча Д.А.; Старожилова Т.К. Учет сезонности в модели Лотки–Вольтерра//Биофизика. 2002, т.47, в. 2., с Перминов В.Д., Саранча Д.А. Фролова А.А. IBM based study of a disease propagation through the lemming population//Second International Conference on Matematical Ecology. Madrid, a, 09-b pp. Перминов В.Д., Саранча Д.А. Об одном подходе к решению задач популяционной экологии//Математичеcкое моделирование. 2003, т.15, N 11, с Бибик Ю.В., Попов С.П., Саранча Д.А. Численное решение кинетического уравнения Богоявленского и системы Лотки-Вольтерра с диффузией//Журнал вычислительной математики и математической физики, том 44, номер C Саранча Д.А. Дискретные отображения и их применение в одной задаче количественной экологии//Сб. научн. тр.: Математика. Компьютер. Образование. М. -Ижевск, 2006, вып. 13, т. 2, с Недоступов Э.В., Саранча Д.А., Чигерев Е.H., Юрезанская Ю.С. О некоторых свойствах одномерных унимодальных отображений//ДОКЛАДЫ АКАДЕМИИ НАУК, 2010, том 430, 1, с. 1–6