Качественные модели водных экологических систем

За последние годы количество аналитических качественных исследований в математической экологии ни в коей мере не уменьшилось, несмотря на всевозрастающие возможности современной вычислительной техники, которая позволяет разрабатывать и реализовывать сложнейшие имитационные модели. Это связано с рядом причин. В частности, при построении той или иной модели мы постоянно находимся в условиях острой нехватки информации об исследуемом объекте: экспериментальные данные либо измерены с некоторой неточностью (а ошибки при проведении гидробиологических исследований, как правило, велики), либо в экспериментах измерены не те параметры, что нужны в модели. Все это приводит к тому, что мы не можем построить имитационную модель лучше, чем те данные, на которых она основывается (Jorgensen et al., 1978; Jorgensen, 1979). Поэтому аналитические модели, допускающие некоторые упрощения внутренней структуры объекта, могут оказаться даже точнее и лучше тех данных, на которых они основаны, так как в основном используют лишь их средние статистические показатели.

Первые интересные аналитические результаты в математической экологии получены в классических работах А. Лотка (Lotka, 1925) и В. Вольтерра (Volterra, 1931; Вольтерра, 1976), независимо друг от друга предложившие модель экосистемы типа “хищник-жертва”, которая послужила демонстрацией существования периодических колебаний численности популяций в стационарной среде, независимо от внешних (в основном, климатических: солнечных, атмосферных и т. п.) факторов, которые первоначально представлялись экологам столь перспективными для объяснения колебаний в экосистемах. В результате сильных упрощений реальной картины автором построена модель сообщества, которая не дает точного описания какой-либо конкретной ситуации, но позволяет провести детальное исследование ее общих свойств. Труд В. Вольтерра заложил фундаментальные основы теории биологических сообществ, построенной именно как математическая теория.

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

1) пища имеется в неограниченном количестве, либо ее поступление с течением времени жестко регламентировано;

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

3) хищные виды поедают жертв, причем в единицу времени количество съеденных жертв всегда пропорционально произведению количества хищников на количество жертв;

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

5) если вид питается пищей, имеющейся в неограниченном количестве, то его размножение подчиняется условию, что прирост численности вида за единицу времени пропорционален ее наличному количеству;

6) если вид питается пищей, имеющейся в ограниченном количестве, то его размножение регулируется тем, что за единицу времени прирост пропорционален количеству съеденной пищи.

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

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

Глубокое аналитическое исследование системы типа “хищник-жертва” проведено А.Н. Колмогоровым (1936) (переработанный и дополненный вариант этой статьи: Колмогоров, 1972), который в своей работе ограничился лишь некоторыми качественными предположениями о характере межвидовых взаимодействий и получил условия, при которых в системе возникает устойчивое стационарное состояние или устойчивый предельный цикл. Однако, необходимого и достаточного условия наличия устойчивого предельного цикла автору получить не удалось.

Большой вклад в общую теорию математического моделирования и теорию моделирования водных экосистем был внесен В.А. Костициным (Kostitzin, 1937) и Г.Ф. Гаузе (Гаузе, Витт, 1934; Гаузе, 1935; Gause, 1935). Многие работы В. Вольтерра, В.А. Костицина имели гидробиологическую направленность, не говоря уже о том, что и первоначальный интерес В. Вольтерра к математической экологии возник после его знакомства с данными о колебании численности популяций рыб, а первые его модели были построены для объяснения того факта, что после длительного перерыва в рыболовстве, вызванного войной, в уловах рыбы из Средиземного моря стали преобладать хищники.

Качественные аналитические модели динамики водных экосистем с успехом развиваются и в настоящее время. Определенный интерес представляют исследования, посвященные вопросам управления биологическими системами (Свирежев, Елизаров, 1972; Зубер и др., 1972; Гончаров, 1984; Spencer, Collie, 1995), оптимизации структуры сообществ (Левич и др., 1996; Левич и др., 1997; Левич, Булгаков, 1997), а также работы, направленные на решение проблемы оптимального сбора урожая, построение алгоритмов поиска оптимальной стратегии кормления и выращивания рыб в рыбоводных прудах (Казанский, Михайлов, 1989; Решетников и др., 1990; Svirezhevet al., 1984).

На основе классических методов теории дифференциальных уравнений исследуют качественные свойства и устойчивость ряда моделей экосистем (Булгакова, 1968; Алексеев, 1976; Тонких, 1985; Базыкин, 1985; Свирежев, 1987; Баренблатт и др., 1993; Виноградов и др., 1993; Медвинский и др., 1997; Завалишин, Логофет, 1977; Петровский и др., 1997; Петровский и др., 1998; Suzuki, Fukuoka, 1979; Matsumura, Sakawa, 1980; Scheffer, 1990; Scheffer, 1991; Steele, Henderson, 1992; Spencer, Collie, 1995; Keitt, 1997; Steffen et al., 1997; Ghosh, Sarkar, 1998; Jonsson, Ebenman, 1998; Petrovskii, 1998; Varriale, Gomes, 1998; Edwards, Brindley, 1999; Jang, Allen, 1999; Murray, Parslow, 1999; Sun, Yang, 1999; Aota, Nakajima, 2000; Gurney, Veitch, 2000; Malchow, 2000; Medvinsky et al., 2000; Buffoni et al., 2001; Ruan, 2001; Venturino, Medvinsky, 2001; Chattopadhayay et al., 2002; Hart, 2002; Malchow et al., 2002; Medvinsky et al., 2002). Обычно в таких моделях рассматриваются отдельные участки трофической цепи, включающие биогенные элементы, детрит, фитопланктон, зоопланктон, реже - рыб. В пространстве параметров модели находят области устойчивости полученных стационарных точек, а затем исследуется динамика фазовых переменных в зависимости от значений этих параметров. Как правило, суммарное количество вещества в экосистеме считается постоянным. Вопрос об устойчивости режима поведения математических моделей природных сообществ имеет большое значение. Изучение этого вопроса важно хотя бы потому, что оно позволяет проанализировать различные механизмы функционирования экологических систем и отвергнуть те из них, которые порождают неустойчивые (негрубые) решения. Неустойчивые режимы просто не могут реализовываться в природе в течение длительного периода времени.

Обычно устойчивость рассматривается как атрибут самой системы, тогда как правильнее рассматривать устойчивость как атрибут того или иного режима функционирования экосистемы (Базыкин, 1985). Это придает актуальность систематическому исследованию таких моделей экосистем, в которых реализуется более одного локально устойчивого режима поведения. Такие модели в наибольшей степени соответствуют описанию качественных перестроек поведения экосистемы в разных ситуациях. Адекватным математическим аппаратом при анализе такого поведения служит качественная теория обыкновенных дифференциальных уравнений и теория бифуркаций (Арнольд, 1978; Хэссард и др., 1985; Эрроусмит, Плейс, 1997).

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

Методами качественного анализа показано также, как неточности, полученные в результате изменений начальных условий и параметров, могут приводить к катастрофическим ошибкам в модельных прогнозах (Воинов, 1981). Удалось решить и ряд принципиальных вопросов, связанных с оценкой погрешностей, которые возникают в результате агрегирования переменных моделей (Caleet al., 1983).

С помощью аналитических и анализирующих моделей исследователи, отталкиваясь от частных вопросов, ищут общие закономерности, присущие рассматриваемой системе, основываясь на той или иной научной гипотезе. Так, например, построив простейшую модель замкнутой экосистемы, описывающую динамику биомасс двух популяций фитопланктона, лимитируемых единственным биогенным элементом, и дополнив ее еще одной переменной, описывающей концентрацию биологически активного метаболита, Ю.А. Домбровский и Г.С. Маркман (Домбровский, Маркман, 1978) показали, как гипотеза о наличии у фитопланктона эктокринного регулирования позволяет объяснить широко известные качественные закономерности динамики фитопланктона (парадокс видового разнообразия фитопланктона, цикличность его динамики). Учет влияния метаболических выделений дает возможность получить неоднородные в пространстве решения (диссипативные структуры), которые имитируют существование устойчивых областей с повышенной плотностью (пятнистость распределения планктона) (Домбровский, Маркман, 1983).

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

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

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

, (1)

. (2)

Здесь P - фитопланктон, H - травоядный зоопланктон, потребляющий фитопланктонные организмы, a - скорость роста, d - скорость отмирания. Потребление P записано в виде простейшей функциональной зависимости , а рост H - как некоторая часть потребления корма. Уравнения (1), (2) неудовлетворительны потому, что они являются качественно негрубыми, приводят к появлению нейтрально устойчивых циклов в P и H и вызывают тем самым иллюзию цикличности динамики популяций (Рубин и др., 1977; Рубин и др., 1987. - 304 с.; Ризниченко, Рубин, 1993; Эрроусмит, Плейс, 1997). Кроме того, эта модель не принимает во внимание неоднородность распределения популяций и факторов, приводящих к ней. Введение в эту модель стохастических изменений, соответствующих изменчивости окружающей среды или популяций, обусловливает появление случайного шума, подавляющего P или H (Bartlett, 1957). По этой причине указанные выше уравнения необходимо модифицировать так, чтобы с их помощью можно было получить либо устойчивые стационарные решения, либо предельные циклы. Для этого можно использовать, например, следующие типы усложнений:

·можно сделать каждое из функциональных соотношений более реалистичным, введя, например, дополнительные члены в правые части уравнений (1), (2);

·можно добавить другие уравнения, в частности, для продуктов питания, от которых зависит фитопланктон;

·можно представить фитопланктон или зоопланктон в виде множества особей, обладающих различной реакцией;

·можно ввести в рассмотрение пространственные характеристики физической среды.

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

Функциональные соотношения. Соотношение между H и P, фигурирующее в (1) в виде , определенно, слишком простое. Согласно экспериментальным данным, при увеличении P скорость выедания на единицу H (например, копепод) стремится к некоторому пределу. Следовательно, можно заменить на . Однако новая система уравнений

, (3)

(4)

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

Что касается скорости выедания, то для нее мы располагаем, по крайней мере, многочисленными экспериментальными данными, которые могут быть использованы при обосновании модели. О справедливости представления эффекта вымирания хищника в виде можно только догадываться. В некоторых моделях смертность хищника задается в виде , или даже в виде , где , и n - эмпирические параметры. Влияние различных способов параметризации этого процесса на качественные особенности динамики двух- и трехкомпонентных систем рассмотрено в работах (Steele, Henderson, 1992; Edwards, Brindley, 1996; Edwards, Brindley, 1999; Murray, Parslow, 1999). На самом деле такая параметризация должна иметь, по крайней мере, вид , где F - выедание H хищниками, например, рыбами. В этом случае необходимо иметь уравнение еще и для F. Таким образом проблема здесь сводится к параметризации F в терминах H. В приведенных выше уравнениях F неявно принималось равным константе и тем самым предполагалось, что F пропорционально количеству хищников при всех концентрациях H. Выбор постоянного значения F можно оправдать для определенных временных масштабов, если допустить, что F растет значительно медленнее, чем H. Однако, как и в случае с P, естественнее было бы предположить, что существует некоторый верхний предел H, и тогда, как и раньше, скорость выедания H можно представить в виде

. (5)

Однако в этом случае мы вновь получаем неустойчивую систему. Более того, остается открытым вопрос о соотношении между P и H в стационарном состоянии . Действительно, уравнение (2) дает и, значит, получается, что состояние фитопланктона не зависит от его собственной скорости роста, характеризуемой в модели коэффициентом a. Используя (5), приходим к выводу об увеличении с уменьшением . Этот вывод может оказаться приемлемым при анализе возможных связей в пределах одной небольшой среды. Однако, имеются данные и о том, что , и изменяются синфазно (Blackburn, 1973). Таким образом, рассмотрение даже такой простой в математическом отношении модели показывает, что проблема параметризации F, которая может казаться второстепенной проблемой замыкания для верхнего уровня модели, в действительности оказывает решающее влияние на ее структуру и, стало быть, определяет временные и пространственные масштабы описываемого моделью явления.

Обобщение моделей. При построении любой экологически полноценной модели первое, на что необходимо обратить внимание, это на введение предела в питании. Для примера рассмотрим только один продукт питания N и очень простую, двухслойную, модель водоема. Предположим далее, что первичная продукция сосредоточена лишь в верхнем слое, а перемешивание между слоями характеризуется некоторой скоростью обмена k. N в нижнем слое остается постоянной величиной и равной , а . Это эквивалентно допущению о том, что верхний слой перемешивания имеет постоянную толщину, а на его нижней границе интенсивность света настолько мала, что роста фитопланктона за счет фотосинтеза в нижнем слое не происходит. Известно, что нехищный зоопланктон способен активно перемещаться из верхнего слоя в нижний, причем это перемещение имеет суточный цикл: в верхнем слое зоопланктон находится, главным образом, ночью, в нижнем - днем. Следовательно, можно предположить, что скорость выедания зоопланктоном фитопланктона в верхнем слое пропорциональна плотности его популяции H. Простейшей моделью, описывающей динамику N, P и H во времени будет модель следующего вида:

, (6)

, (7)

. (8)

Здесь N, P и H измеряются в единицах содержания продуктов питания (например, в мг/м³). Предполагается, что продукты питания, поедаемые H, но не используемые для роста, возвращаются в воду в верхнем слое.

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

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

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

Пространственные эффекты. Для пояснения характера реакций вернемся к простейшему их описанию:

, (9)

, (10)

где используется очень простое одномерное представление горизонтальной турбулентной диффузии с коэффициентом . Эффект диффузии должен способствовать устойчивости этих уравнений к любым осциллирующим возмущениям (Марри, 1983; Murray, 1975). Другими словами, средние значения P и H в интервале изменения x по-прежнему будут изображаться нейтрально устойчивыми циклами, но не такими, как в случае с P, когда подавлялись все возмущения с длиной волны, отличающейся от критической. Отличительной особенность уравнений (9), (10), как и всех достаточно развитых моделей экосистем, являются нелинейные члены, которые характеризуют взаимодействие между трофическими уровнями. Благодаря таким членам возмущения в любом небольшом диапазоне длин волн будут вызывать изменения на всех других длинах волн. Нелинейные эффекты могут быть связаны также и с возмущениями одного из коэффициентов системы в фиксированном интервале длин волн. К числу экологических проблем относится также вопрос о способах введения в модель возмущений. В частности, уравнения (9), (10) неадекватны в том смысле, что представление движения в верхнем слое в виде простой горизонтальной диффузии не соответствует особенностям биологических процессов.

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

Основой всех моделей динамики трофических сообществ является понятие трофической функции, описывающей индивидуальный рацион хищника. По определению, индивидуальный рацион есть количество жертв, потребляемых хищником за единицу времени. Вид функциональной зависимости между биомассой хищников и обилием их жертв имеет большое значения для изучения механизмов функционирования экосистем и до сих пор является темой продолжающихся споров в экологическом моделировании. В работе (Steele, Henderson, 1992) на примере пяти предельно простых N-P-Z-моделей, включающих в себя всего три переменные: питательные вещества, фито- и зоопланктон, показано, насколько важна правильная параметризация отдельных процессов в модели, в частности процессов трофических взаимодействий между компонентами экосистемы и процессов отмирания организмов. Выбор той или иной функциональной зависимости и соответствующий набор параметров оказывают очень большое влияние на функциональный ответ (тип и устойчивость особых точек, тип траектории на фазовой плоскости) всех без исключения звеньев такой системы и приводит к принципиально отличной временной динамике как питательный веществ, так и планктонных циклов. А это, в свою очередь, приведет к различным структурам энергетических потоков в водной экосистеме. Даже для относительно простых модельных ситуаций (например, в данной работе не рассматривались процессы оседания планктона, не учитывалось бактериальное звено, возрастная структура зоопланктона, разнообразие видов в пределах трофического уровня, субстратное лимитирование, вертикальное распределение фитопланктона и миграции зоопланктона, процессы турбулентного перемешивания) такого рода проблемы решать не так просто, поскольку заранее неизвестен конкретный механизм, управляющий динамикой трофических взаимодействий между видами и разными трофическими уровнями.

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

, (11)

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

, (12)

В - биомасса хищника. Для конкретизации вида корректирующего множителя предложено несколько функциональных зависимостей. Ниже будут рассмотрены наиболее часто используемые из них при построении математических моделей. В уравнении (11) функция pred записана в разностном виде. Однако ее аналогично можно записать и в дифференциальной форме.

Широко используются три типа функциональных зависимостей, которые были впервые введены Холлингом (Holling, 1959). Первый тип описывает широко известную динамику Лотка-Вольтерра. В наиболее простом случае предполагается, что функция зависит только от количества пищевого ресурса F, который доступен данному виду хищника (Prey-Dependent function). Классической формой трофической функции является зависимость индивидуального рациона хищника только от плотности популяции жертв F, то есть . Простейшим случаем PD трофической функции является классическая линейная зависимость Лотка - Вольтерра. Линейность является идеализированным представлением случайных столкновений хищников с жертвами в однородно распределенном сообществе (Романовский и др., 1984). Несмотря на корректность такого приближения при моделировании взаимодействий видов в однородных экосистемах, свойства точечной модели Лотка - Вольтерра противоречат наблюдениям за природными сообществами. Прежде всего трофическая функция не учитывает насыщения рациона хищника при увеличении плотности популяции жертв. Более адекватным является использование функции с фиксированным верхним пределом для того, чтобы выполнялось условие (12):

(13)

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

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

. (14)

Динамика роста в этом случае характеризуется константой полунасыщения К. Очевидно, что если F = K, то и скорость роста составляет половину от своего максимально возможного значения. Наиболее известная PD-модель, описывающая эффект насыщения, была предложена Холлингом (Holling, 1959): . Здесь a - коэффициент эффективности поиска жертв хищником, количественно характеризующий интенсивность его атак, h - величина, обратная максимальному индивидуальному рациону. Эффект насыщения индивидуального рациона эмпирически установлен для многих видов.

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

. (15)

Другие, реже используемые зависимости, имеют вид:

, (функция Гаусса) (16)

(функция Ивлева). (17)

Трофическая функция В.С. Ивлева была успешно идентифицирована как для некоторых видов рыб (Ивлев, 1955), так и для многих ракообразных с различными типами пищевого поведения (Сущеня, 1975). При малых значениях переменной F формула Ивлева эквивалента функции Холлинга, и обе зависимости насыщаются с ростом плотности жертв. Однако учет насыщения индивидуального рациона не устраняет всех несоответствий натурных данных прогнозам моделей, основанных на использовании PD трофических функций.

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

, (18)

так что

. (19)

Здесь предполагается, что при выедании жертв их биомасса без потерь трансформируется в биомассу хищников. При функция становится линейной функцией биомассы хищника В, если рассматриваются функции Холлинга первого, второго типов и функция Ивлева: . Это означает, что эффективность поиска хищником жертв не уменьшается при очень низких плотностях последних. Если рассматривается функция Гаусса, то , а значит эффективность поиска хищником редких жертв становится гораздо более высокой. Функция Холлинга III-го типа, наоборот, показывает, что , то есть эффективность хищничества снижается. Экологические причины для уменьшения эффективности хищничества могут состоять в следующем:

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

2) поиск становится гораздо менее эффективным, поскольку при очень низкой плотности жертв их трудно найти и поймать;

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

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

Рассмотрим в качестве примера точечную модель динамики системы хищник - жертва с логистическим воспроизводством жертв и трофической функцией Холлинга второго типа:

(20)

Здесь P - плотность популяции жертв; H - плотность популяции хищников; r - темп естественного прироста популяции жертв; K - максимальная плотность популяции, достигаемая жертвами при отсутствии хищника (емкость среды, характеризующая первичную продуктивность трофической системы); e - коэффициент эффективности хищничества, определяемый как отношение удельной скорости прироста популяции хищника к его индивидуальному рациону; m - коэффициент смертности хищника. Помимо тривиального равновесия , , соответствующего вымиранию обоих видов, и равновесия , , описывающего вымирание хищника, при выполнении условия

(21)

в системе (20) существует также и нетривиальное равновесие:

. (22)

Известно, что с ростом значения параметра емкости среды K, изначально устойчивый узел (22) становится устойчивым фокусом. Дальнейший рост значения параметра K приводит к бифуркации Пуанкаре - Андронова - Хопфа (Хэссард и др., 1985): равновесие (22) теряет устойчивость, что сопровождается рождением в его окрестности устойчивого предельного цикла. Амплитуда и период цикла быстро возрастают с увеличением параметра емкости среды. Колебания хищников и жертв, возникающие при увеличении продуктивности популяции нижнего трофического уровня сообщества, наблюдаются в лабораторных условиях. Другими словами, динамика искусственных трофических систем может быть достаточно адекватно описана при помощи точечных PD-моделей (Gause, 1935; Luckinbill, 1973; Bohannan, Lenski, 1997; Bohannan, Lenski, 1999; Kaunzinger, Morin, 1998). Однако крупномасштабные экосистемы не обладают этим свойством: увеличение первичной продукции повышает стабильность природных сообществ. Это несоответствие динамических свойств PD-моделей наблюдениям за природными трофическими сообществами получило название парадокса эвтрофирования (paradox of enrichment) (Rosenzveig, 1971). Аналогичная последовательность бифуркационных переходов и возникновение устойчивого предельного цикла с быстро растущей амплитудой наблюдается и при увеличении прожорливости хищника a или его эффективности e. С помощью других PD-моделей также нельзя получить одновременно низкое и устойчивое равновесие численности жертв. Это привело к так называемому парадоксу биоконтроля (biological control paradox), заключающемуся в невозможности смоделировать ситуацию успешного биоконтроля, при которой использование эффективного хищника или паразита позволяет контролировать популяцию вредителя на низком уровне численности (Luck, 1990; Arditi, Berryman, 1991).

Разрешение этих противоречий и необходимость согласования модельных свойств с натурными наблюдениями привели к развитию теории динамики трофических сообществ и появлению большого количества работ, описывающих различные поведенческие механизмы, включение которых в модели позволяет стабилизировать взаимодействия хищник - жертва. В качестве ключевых факторов стабилизации трофических систем называются пространственная неоднородность и агрегация хищников в местах скопления жертв, наличие убежищ для жертв, мутуализм и псевдо-интерференция хищников. Различные способы моделирования факторов, повышающих стабильность трофических сообществ, обсуждаются в работах Карейвы и Хастингса (Kareiva, 1990; Hastings, 1990).Одним из методов является модификация трофической функции: включение в нее плотности популяции хищника позволяет неявно учесть последствия пространственных взаимодействий популяций в рамках точечных моделей.

Наиболее простой вид трофической функции, зависящей от плотности хищника, был предложен Ардити и Гинзбургом (Arditi, Ginzburg, 1989). Они предположили, что пространственная неоднородность приводит к зависимости индивидуального рациона хищника не просто от численности жертв и хищников, а от их отношения (Arditi, Saїah, 1992; Akзakaya et al., 1995):

. (23)

Модели с такими свойствами будем называть RD-моделями (Ratio-Dependent). Наличие или отсутствие пространственной неоднородности обусловливает использование либо RD-, либо PD трофической функции для неявного описания агрегированной динамики сообществ.

RD-модели хищник-жертва разрешают как парадокс биоконтроля, так и парадокс эвтрофирования, позволяя без детального моделирования поведенческих механизмов описать сосуществование жертвы и хищника при низком уровне численности жертв (Arditi, Berryman, 1991). Модель хищник-жертва с трофической функцией Ардити - Гинзбурга - Холлинга (Arditi, Ginzburg, 1989) позволяет получить периодические колебания плотностей популяций и вымирание обоих видов (Jost et al., 1999; Berezovskaya et al., 2001).

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

Постепенно стало весьма очевидным, что необходимо использовать двумерные варианты трофической функции (A Two-Dimensional function). К сожалению, функцию практически невозможно измерить непосредственно в лабораторных экспериментах или с помощью полевых наблюдений. Кроме того, мы должны всегда понимать, что функциональная связь между хищниками и жертвами, на самом деле, отражает множество процессов их взаимодействия (Михайловский, 1988; Михайловский, 1992), которые мы пытаемся объединить в одном соотношении. Поэтому не должно быть неожиданностью то, что мы при различных условиях будем получать разные функциональные зависимости. Двумерная функция, основанная на третьем типе функции Холлинга, которая, как было показано выше, учитывает уменьшение эффективности хищничества при низком обилии ресурса, может быть записана в следующем виде:

, (24)

где В - биомасса хищника, К₂ - константа полунасыщения, К₃ - коэффициент, учитывающий степень зависимости от хищника,

, (25)

F - общее количество пищевого ресурса, К₁ - константа полунасыщения. F определяет только ту долю ресурса, который доступен хищнику. Функция Холлинга второго типа для F показывает, что при низком обилии ресурса (жертв) он менее доступен для потребления. При большая часть жертв защищена от хищничества. Из (24), (25) получим:

. (26)

Принимая, ради упрощения, что

2К₁ = 1/2K₂ = K, (27)

из (26) получим:

. (28)

Тогда функция будет иметь следующий вид:

. (29)

Введенные функции имеют следующие свойства.

·Для любых :

1) ;

2) .

·При :

1) ;

2) ;

3) ;

4) .

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

Наиболее важным и сложным шагом в математическом моделировании экосистем является выделение (или агрегирование) модельных переменных. Чисто биологический подход, основанный на таксономической классификации, в данной задаче оказывается неприемлемым. Когда мы имеем дело с большими экосистемными моделями, возникает потребность в формировании гораздо более крупных единиц. Одним из путей решения указанной проблемы является использование аллометрического масштабирования относительно веса организмов (Gin et al., 1998; Bussenschutt, Pahl-Wostl, 2000), поскольку, как известно, многие физиологические параметры можно описать степенным законом: (Barenblatt, Monin, 1983). Масштабный коэффициент зависит от рассматриваемого параметра, знак в показателе степени показывает, увеличивается или, наоборот, уменьшается параметр р с увеличением веса W организма. Физиологические скорости характеризуются отрицательными показателями степени. В рамках указанного подхода организмы разбиваются на дискретные весовые классы (Bussenschutt, Pahl-Wostl, 2000). Эти классы равномерно распределены на оси веса организмов, выраженном в логарифмическом масштабе. При переходе от одного класса к следующему вес удваивается. Предполагается, что все организмы, относящиеся к одному весовому классу, имеют одинаковые физиологические скорости процессов. Они различаются только с точки зрения их трофических взаимодействий и пищевых предпочтений. Трофические связи каждого компартмента определяются коэффициентами селективности, которые учитывают пищевые предпочтения видов. Важным преимуществом аллометрического способа агрегирования является небольшое количество свободных параметров. Так, если модели строятся на основе таксономического принципа агрегирования, то число параметров модели возрастает почти линейно, а иногда и быстрее, с ростом числа компартментов. В то же время аллометрическая модель может быть дополнена без добавления в нее новых параметров. Это может иметь принципиальное значение, когда моделируется, например, процесс инвазии новых видов в экосистему. Представляется разумным распространить данную аллометрическую концепцию и на рассмотрение временной динамики процессов. Организмы, которые отличаются друг от друга по размерам, не только растут с разными скоростями, но и их целостное взаимодействие со средой и другими организмами происходит на разных временных масштабах. Эти множественные временные масштабы включаются в модельное рассмотрение путем введения временных шагов, которые, как и физиологические параметры, масштабируются в зависимости от веса организмов по степенному закону при .

Если виды питаются на k весовых классах, то

, (30)

где - коэффициент селективности для i-го класса, , - биомасса в i-ом весовом классе, К₁ - константа полунасыщения, которая от i не зависит. Тогда

. (31)

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

. (32)

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

, (33)

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

, (34)

а частично выводятся из системы . Пул питательных веществ непрерывно обновляется за счет поступления из внешних источников и за счет поступлений из отдельных компартментов.

В реальных пищевых цепях взаимодействия компонентов экосистемы, а также сила этих взаимодействий не являются случайными, а имеют вполне определенную структуру, которая возникает (помимо всего прочего) вследствие распределения организмов по их размерам. Такого рода структурированность может влиять на устойчивость всей цепи. В работе (Jonsson, Ebenman, 1998) приводятся доводы в пользу использования отношений размеров тела и аллометрии для аппроксимации значений некоторых параметров в моделях Лотка - Вольтерра. Основная цель состояла в том, чтобы показать, как можно избежать динамической неустойчивости в реальных сообществах, принимая во внимание более реалистические распределения организмов по их размерам. В частности показано, как такого рода распределения влияют на локальную устойчивость в линейных моделях пищевых цепей. Если средние (относительные) отличия между размерами хищников и их жертв уменьшаются в зависимости от их позиции в пищевой цепи (с ростом трофического уровня), то это приводит к повышению локальной устойчивости такой системы (по сравнению с моделью с постоянным отношением размеров хищников и их жертв).

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

(35)

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

. (36)

Используя простые энергетические балансовые оценки, нетрудно получить, что

, (37)

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

, (38)

а определяется в зависимости от заданного трофического уровня в пищевой цепи:

. (39)

Скорости смертности видов задавались в соответствии с предположением о том, что .

Для простейших математических моделей, допускающих аналитическое исследование и, в то же время, дающих качественно верное описание основных процессов в экосистеме, в работе Н.Н. Моисеева и Ю.М. Свирежева (Моисеев, Свирежев, 1979) было предложено название “минимальные модели”, в котором подчеркивается стремление уловить в модели некую основную “минимальную” структуру, отвечающую за наиболее существенные динамические процессы в экосистеме. Минимальные модели являются одним из перспективных направлений развития математической экологии вообще и математического моделирования водных экосистем в частности. Такие модели позволяют подкрепить большую имитационную модель некоторым качественным анализом, рассмотреть топологию фазового пространства и основные ее динамические свойства, что в значительной мере облегчает дальнейшую работу с ней. Заметим, что математическая формализация минимальной структуры, очевидно, не единственна, и поэтому минимальных моделей может быть несколько для каждого конкретного объекта.

К числу детализированных минимальных моделей пресноводных экосистем можно отнести работу А.А. Воинова и Ю.М. Свирежева (Voinov, Svirezhev, 1984). В ней, на основе как общеэкологических соображений, так и асимптотического анализа процессов трансформации вещества, было показано, что динамику эвтрофицируемого фитопланктонного озера можно качественно описать системой четырех дифференциальных уравнений, в качестве фазовых переменных которой взяты концентрации фитопланктона, детрита, биогенных элементов и кислорода. Анализ и исследование динамики стационарных точек модели проведен при изменении суммарного количества вещества в экосистеме (так называемый “квазистационарный процесс”). Ввиду сильной агрегированности модели она не дает количественного описания процесса эвтрофикации, но позволяет качественно проследить изменение концентрации фитопланктона в озере в результате воздействия меняющейся антропогенной нагрузки.

Одним из подходов, реализующих концепцию минимальной модели в экологии стал компартментальный анализ (Завалишин, Логофет, 1997). Экосистема при этом разбивается на блоки, содержащие определенные запасы вещества и обменивающиеся им между собой и окружающей средой. На основе биологической информации задаются скорости обмена, а также скорости входных и выходных потоков. Получающаяся модель называется компартментальной, а блоки - компартментами. Достоинство такого подхода состоит в том, что, во-первых, нет нужды кропотливо собирать данные о взаимодействиях многих видов, обитающих во всякой реальной экосистеме, во-вторых, исследователь относительно свободен в выборе модельных переменных и предмета обмена (вместо биомассы организмов можно оценивать, например, концентрацию какого-либо биогенного элемента). На основе таких представлений в биологии разработаны методы экспериментального определения потоков и запасов вещества в естественных условиях (Сорокин, 1971; Виноградов, Шушкина, 1987; Сорокин, Джелли, 1996), в результате применения которых концептуальная схема круговорота вещества превращается в численную диаграмму: в блоках указываются запасы, а стрелки символизируют потоки.

В конце 70-х - начале 80-х годов прошлого века был предложен метод анализа этих статических диаграмм, получивший название энвирон-анализа (Patten, 1982; Patten, 1985). Входной или выходной энвирон данного компартмента - это совокупность блоков, с которыми он связан входящими или исходящими потоками. Используя матричные соотношения этого метода, можно вычислить, какая доля выходящего из i-го компартмента потока попадает в j-й, оценить среднее время, проводимое веществом в компартменте, сделать количественные выводы о роли того или иного процесса (например, дыхания) в функционировании экосистемы или круговороте химических элементов.

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

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

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

(40)

где P - биомасса фитопланктона, Z - биомасса зоопланктона, N концентрация биогенного элемента, с - коэффициент самолимитирования фитопланктона, e - коэффициент эффективности усвоения пищи зоопланктоном, F - максимальная скорость выедания рыбами зоопланктона, m - скорость смертности зоопланктона, - максимальная скорость роста фитопланктона, p - максимальная скорость выедания фитопланктона зоопланктоном, , , - константы полунасыщения. Уравнение, задающее в модели динамику фитопланктона, основано на логистическом уравнении роста с учетом биогенного лимитирования и соответствует кинетике Моно. Скорость выедания фитопланктона зоопланктоном зависит от концентрации микроводорослей. Таким образом, в данном случае используется PD-модель, описывающая эффект насыщения, и, соответственно, трофическая функция Холлинга второго типа. Этот процесс также соответствует кинетике Моно. Выедание зоопланктона рыбами описывается в модели трофической функцией Холлинга третьего типа. Зависимость скорость выедания от плотности зоопланктона имеет, следовательно, сигмоидальную форму. Использование именно такой зависимости позволяет учесть снижение эффективности выедания зоопланктона рыбами при малой численности зоопланктона. Функциональная реакция этого типа свойственна сравнительно высокоорганизованным хищникам, которые способны обучаться в ходе поиска жертвы, а именно, наращивать интенсивность охоты при увеличении плотности пищи. Такая реакция возникает и в тех случаях, когда хищник способен переключаться на альтернативную пищу.

В результате качественного анализа поведения модели и некоторых численных экспериментов с нею было показано следующее:

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

·постепенное увеличение плотности рыб в эвтрофных условиях (то есть при относительно высоких концентрациях биогенного элемента) может привести к скачкообразным изменениям плотностей фито- и зоопланктона;

·присутствие рыб ослабляет колебания в системе фитопланктон - зоопланктон, возникающие в условиях эвтрофирования;

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

Детальное исследование типов устойчивости данной системы выполнено в работе (Malchow, 1993). Показано, что для различных комбинаций параметров в такой системе возможны следующие типы стационарных состояний:

·при высокой скорости выедания зоопланктона рыбой стационарным состоянием системы является устойчивый узел и при этом формируется более высокая плотность фитопланктона;

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

·последующее увеличение плотности зоопланктона в такой бистабильной системе приводит к формированию устойчивого фокуса и возникновению бистабильности типа устойчивый узел - устойчивый фокус;

·дальнейшее увеличение плотности зоопланктона приводит к постепенному исчезновению указанной выше бистабильности в системе и формированию устойчивого узла с более высокой плотностью зоопланктона;

·при очень низкой (или нулевой) скорости выедания зоопланктона рыбами в такой системе может происходить постепенное формирование неустойчивого фокуса который, в результате бифуркации Пуанкаре - Андронова - Хопфа (Хэссард и др., 1985), сменяется предельным циклом.

Такая последовательность смены стационарных состояний в системе происходит в результате постепенного уменьшения выедания зоопланктона рыбами. Эта модель является прекрасной иллюстрацией влияния контроля “сверху” на состояние трофической цепи и особенностей ее временной динамики и устойчивости. Плотность рыб является важным фактором, который может как ослаблять колебания в системе фитопланктон - зоопланктон, так и приводить к скачкообразным изменениям динамики плотностей фито- и зоопланктона. Вместе с тем, не следует забывать о том, что типы динамического поведения системы будут зависеть не только от вида используемых функциональных зависимостей между переменными модели, но также и от структуры самой модели: количества модельных переменных, направленности и силы связей между ними (Sazykina et al., 2000; Hart, 2002).

Еще одним важным параметром, который может определять качественную перестройку временной динамики системы, является максимальная скорость роста фитопланктона . Ее величина сильно зависит от температурных и световых условий в водоеме. В работе (Steffen et al., 1997) исследовано влияние сезонной изменчивости этого параметра на поведение модельной системы (40). Предполагалось, что величина изменяется периодически в течение года в соответствии со следующим уравнением:

. (41)

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

Добавление в исходную модель (40) диффузионных слагаемых позволяет изучать особенности уже не только локальной, но и пространственно-временной динамики. В частности, в (Malchow, 1993; Malchow, 1994) были рассмотрены необходимые условия и показана возможность образования пространственно-временных структур типа бегущих волн планктона и структур Тьюринга. В работах (Медвинский и др., 1997; Медвинский и др., 2002; Medvinsky et al., 2000) рассматривается механизм самопроизвольного рождения и гибели вихря (в двумерных средах называемого также ревербератором). Была показана достаточность наличия подвижности рыбной стаи для формирования сложных пространственных структур (в том числе и вихревых) в популяциях планктона, а также продемонстрирована роль ревербератора в упорядочении локальной плотности планктона. В отношении экологических систем подобного рода исследования были выполнены, по всей видимости, впервые. В отличие от более ранних моделей, в которых предполагалось, что F является константой, в модели (Медвинский и др., 1997) от столь сильного ограничения авторы отказались, поскольку нельзя пренебрегать пространственной структурированностью рыбных популяций. Нельзя пренебрегать и таксисными реакциями рыб на изменения условий окружающей среды. Для учета этих факторов было принято, что: 1) максимальная скорость F потребления зоопланктона рыбой представляет собой не постоянную величину, а функцию пространственных координат; 2) поведение стаи рыб подчиняется так называемым правилам Эбенхоха (Ebenhцh, 1980), согласно которым стая после некоторого лаг-периода смещается в область с более высокой концентрацией пищи (зоопланктона), при условии, что концентрация пищи в пределах текущего месторасположения стаи оказалась ниже порогового значения. Рождение вихря происходит в тот момент, когда изначально однородные распределения фито- и зоопланктона под влиянием рыбной стаи, блуждающей по сложной траектории в соответствии с правилами Эбенхоха, становятся существенно неоднородными. Планктонный вихрь, как и вихри в любых активных средах, формируется в результате блокирования распространения автоволны в одном из направлений. Отметим, что если возникновение планктонного вихря связано с воздействием на планктон стаи рыбы, то в дальнейшем циркуляция планктонной волны не требует присутствия рыбы в области, охваченной вихрем. В этот период будут работать диффузионные слагаемые модели. В общем случае при этом образуются сложные, апериодические пространственные структуры, а локальные изменения во времени концентраций как фито-, так и зоопланктона оказываются неупорядоченными, хаотическими. Обратим внимание на то, что образование упорядоченных структур в данном случае никак не связано с воздействием внешних гидродинамических факторов, хотя они также могут играть важную роль в формировании пространственной неоднородности распределения планктона и рыб (Wroblewski, O’Brien, 1976; Steele, Henderson, 1979; Steele, Henderson, 1992; Malchow, Shigesada, 1994; Venturino, Medvinsky, 2001).

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

Предлагаемый способ описания поискового поведения хищника основан на использовании систем дифференциальных уравнений с частными производными типа реакция - адвекция - диффузия. В общепринятых моделях таксиса направленные перемещения хищника описываются адвективной компонентой скорости, которая в каждой точке пространства пропорциональна градиенту плотности популяции жертв (Березовская и др., 1999; Березовская, Карев, 1999). Подход, используемый в данной работе, отличается гипотезой относительно моделирования таксиса: предполагается, что поисковое поведение хищника определяется не скоростью как таковой, а ускорением, которое в каждой точке пространства пропорционально градиенту плотности распределения жертв, а в общем случае - градиенту некоторого миграционного стимула (Говорухин и др., 2000). Модель хищник-жертва типа реакция - адвекция - диффузия, основанная на гипотезе пропорциональности ускорения перемещения хищников градиенту плотности жертв, позволяет получать пространственно-неоднородные режимы даже для случая отсутствия воспроизводства и смертности хищников. Иными словами, такая модель способна описать пятнистость распределения сообщества, вызываемую пространственным поведением хищника (Тютюнов и др., 2001), а не демографическими процессами.

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

, (42)

где - плотность популяции жертв в точке в момент времени t; - скорость перемещения хищников; k - коэффициент поисковой активности хищника, характеризующий его чувствительность к неоднородности распределения плотности жертв; - коэффициент диффузии скорости. Диффузионный член в (42) интерпретируется как результат социального поведения особей: стайные эффекты выравнивают величины и направления скоростей находящихся рядом хищников (Flierl et al., 1999). Модель сообщества представляет собой систему уравнений в частных производных типа реакция - адвекция - диффузия. Локальная кинетика: воспроизводство видов, выедание жертв и смертность хищников моделируются как процессы, происходящие в небольших и, следовательно, однородных элементарных участках рассматриваемого ареала. Поэтому, следуя предположению Ардити и Гинзбурга (Arditi, Ginzburg, 1989) о соответствии PD трофических функций однородности распределения сообщества, локальное взаимодействие популяций будем описывать функцией Холлинга (Holling, 1959). Пространственная модель сообщества имеет следующий вид:

(43)

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

. (44)

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

Установление неоднородных режимов в модели зависит от баланса между направленными и случайными перемещениями хищника. В рассматриваемой модели рост поисковой активности хищника приводит к адвективной неустойчивости: рождению неоднородного периодического режима в окрестности однородного равновесия. рождению неоднородного периодического режима в окрестности однородного равновесия. Неоднородный режим усложняется с увеличением параметра k: цикл ® сложный цикл ® тор ® хаос. Численными экспериментами показано также, что при высоких значениях параметров интенсивности атак a и поисковой активности хищника k, неоднородный режим, реализующийся в системе, мало зависит от вида трофической функции. Определяющим фактором для возникновения такого режима является не нелинейность трофической функции Холлинга, а колебательная потеря устойчивости однородного нетривиального равновесия.

Предложенный подход к моделированию поискового поведения хищника позволяет описать пространственную динамику трофического сообщества как непрерывное движение пятен плотности взаимодействующих популяций. Согласно уравнению (42), хищник реагирует на неоднородность распределения популяций жертв, изменяя скорость своего движения в направлении градиента их плотности, что и описывает агрегированное поведение хищника. Достигнув точки максимальной плотности жертв, хищник замедляет свое движение, так как его ускорение меняет направление на противоположное. Агрегирование хищников приводит к локальному вымиранию жертв, в то время как в пятнах с низкой плотность хищников образуются локальные убежища, где происходит рост плотности популяций жертв и, согласно (42), хищник ускоряется в направлении вновь образовавшихся скоплений. Тем самым модель, основанная на гипотезе (42), явно описывает комплекс взаимосвязанных факторов, необходимых и достаточных для стабильности взаимодействий хищник-жертва. Благодаря этому модель (43) способна воспроизвести реальную ситуацию сосуществования хищника и жертвы на низком уровне плотности популяции жертвы. Высокая поисковая активность является механизмом, позволяющим популяции хищника адаптироваться к дефициту жертв. Активно мигрирующий вид способен повысить запас пищевого ресурса и, как следствие, увеличить свой рацион. Однако чрезмерно высокие величины параметра k приводят к установлению хаотической динамики с нерегулярными вспышками численности жертв. Фактически это означает, что хищник, обладающий слишком высокой поисковой активностью, не способен эффективно регулировать численность жертв. Представленная модель описывает динамику сообщества хищник-жертва на одномерном изолированном ареале. Двумерный вариант рассмотрен и исследован в работе (Arditi et al., 2001). В классе неявных точечных моделей парадоксы биоконтроля и эвтрофировния разрешаются, в частности, с помощью RD трофической функции (Arditi, Berryman, 1991), применение которой обусловливается влиянием пространственной неоднородности на динамику сообщества. В данной модели пространственная неоднородность не задана априори, а генерируется поисковым поведением хищника, что позволяет использовать даже самые простые PD-функции для описания трофических взаимодействий.

Неоднократно показано, что даже простые модели могут демонстрировать довольно сложное поведение, а при определенных условиях и хаотическую динамику. Тем не менее, подтверждения такой сложной динамики в природных популяциях редки. Возможно, что одной из причин этого является воздействие возмущений, как случайных, так и неслучайных. Ранее проведенные теоретические исследования показали, что хотя случайные возмущения и могут несколько видоизменять параметры периодических орбит некоторых простых популяционных моделей, таких, например, как логистическая модель или модель Морана-Рикера, однако совсем необязательно, что в таких системах действительно будут происходить существенные изменения их свойств (Crutchfield et al., 1982; Schaffer et al., 1986). В других работах, наоборот, обнаружено, что случайные возмущения могут приводить к изменениям некоторых свойств популяционных моделей (Ellner, Turchin, 1995; Kaitala et al., 1997). Хотя экологи и осознают, что случайные возмущения могут оказывать гораздо более сложное воздействие на динамику поведения популяционных моделей, но пока нет никаких математических методов для того, чтобы определить это, а если свойства системы и могут меняться, то не ясно, какого рода изменения могут и действительно будут происходить в системе.

В работе (Sun, Yang, 1999) рассматриваются двенадцать хорошо известных популяционных моделей, строятся бифуркационные диаграмы при воздействии случайных возмущений. Рассматривались случайные возмущения двух типов. В одном из них возмущения зависели от плотности популяции, в другом - нет. В обоих типах возмущение имело нормальное распределение с нулевой средней. Показано, что при воздействии возмущений в шести из рассмотренных моделей могут происходить существенные изменения свойств: 1) периодические циклы изменяться на хаотическое поведение; 2) хаос может измениться на периодические циклы; 3) может изменяться амплитуда колебаний; 4) сложная динамика может меняться на простую (правда это показано в исследованиях физической системы); 5) динамика системы может представлять собой совокупность нечетко выраженных циклов. Указанные изменения наиболее ярко выражены в многомерных моделях. В пяти моделях, в которых не выявлено существенных изменений свойств при воздействии случайных возмущений, тем не менее, имело место размазывание периодов в циклах. Все это говорит о том, что параметры моделей, которые мы определяем с помощью различных методов исходя из наблюдений реальных систем, могут не вполне характеризовать именно ту систему, которую мы на самом деле изучаем. Они в значительной степени зашумлены. Другой важный теоретический вывод, который можно сделать из проведенных исследований состоит в следующем. Если при воздействии случайных возмущений в системе произошли качественные изменения и появились элементы хаотического поведения, то прогноз такой системы исходя из наблюдений невозможен, по крайней мере в течение того периода времени, по истечении которого система не перейдет в новое состояние с предсказуемыми свойствами.

Система нелинейных уравнений типа конвекция - диффузия - реакция (или, в некоторых случаях, только диффузия - реакция) является важнейшей базовой минимальной математической моделью, применяемой для изучения водных экосистем (Марри, 1983; Романовский и др., 1984; Свирежев, 1987; Murray, 1975; Okubo, 1980). Вследствие нелинейности рассматриваемой модели даже малое локальное возмущение может стать причиной глобальных пространственно-временных изменений в системе. На основе таких моделей, например, в работах (Виноградов и др., 1993; Баренблатт и др., 1993; Petrovskii, 1998) изучается эволюция океанического планктонного сообщества после локализованного воздействия, нарушающего нормальное функционирование сообщества. В результате исследования математической модели, описывающей функционирование планктонной системы из девяти элементов (фитопланктон, бактерии, простейшие, зоопланктон, медузы, рыбы, РОВ, детрит, минеральный фосфор), показано, что при заданных биологических параметрах сообщества и определенных динамических свойствах среды (речь идет прежде всего о влиянии адвекции) в зависимости от размеров и интенсивности оказанного воздействия (которое называется повреждением или импактом) возможны два качественно различных режима. Если для заданной интенсивности импакта его размер меньше некоторого критического, то повреждение сохраняется только в ограниченной области и величина его убывает со временем - система восстанавливает свои свойства. Если же начальный размер импакта превышает критический, повреждение распространяется с постоянной скоростью, причем это распространение имеет характер бегущей волны. Моделировался случай внесения субстанта в период цветения фитопланктона. Было показано, что начальное возмущение сначала затухает, но затем, по мере развития фитопланктона, начинает распространяться. С некоторого момента происходит выход на волну, распространяющуюся с постоянной скоростью. Волна имеет в данном случае и задний фронт. После прохождения волны за задним фронтом довольно резко наступает “опустошение” в поле фитопланктона, распространяющееся по пространству в виде волны, затем постепенный его рост в центре, и с некоторого момента снова образуется импактная волна “опустошения” в поле фитопланктона. Скорость распространения импактной волны в численных экспериментах имела порядок 10 см c^-1. В реальной экосистеме (например такой, как экосистема Черного моря) это может означать, что внесенный локально в море субстант мог бы за счет подобного механизма взаимодействия с остальными компонентами экосистемы распространиться по всей акватории моря за несколько месяцев.

Механизм динамического удержания локализованного начального возмущения (ЛНВ) в системе “хищник-жертва” исследован в работе (Петровский и др., 1997). Динамика биологического сообщества описывается системой нелинейных уравнений в частных производных (рассматривается одномерный случай):

, (45)

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

Рассмотрим динамику двухкомпонентной активной среды, описываемой системой уравнений с нелинейностью типа “хищник - жертва”:

, (46)

, (47)

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

Динамика системы “хищник - жертва” в случае финитных начальных условий принципиально зависит от типа популяции жертвы, а именно, имеет ли популяция u критический порог плотности или нет. В рамках модели (46) - (47) отсутствие порога означает, что функция обладает следующими свойствами: , ; ; , . При отсутствии хищника в этом случае любая локальная “вспышка” численности жертвы приводит к возникновению распространяющейся неограниченно популяционной волны переключения (Свирежев, 1987). Поскольку при малых u, v нелинейность в уравнениях (46), (47) соответствует классической модели Вольтерра, то пресс со стороны хищника ни при каких начальных условиях не приводит к полному вымиранию жертвы. Наличие хищника не изменяет принципиально тип режима, но может менять характер распространения, приводя к появлению периодических популяционных волн. Существенно иная картина имеет место для случая популяции с критическим порогом плотности, когда для функции имеем: , ; , ; ; , . При отсутствии хищника ЛНВ численности популяции u приводит к ее неограниченному распространению (также сопровождающемуся движением границы заселенной области в виде волны) только при выполнении для начального распределения популяции определенных критических условий - достаточно больших протяженности и амплитуды ЛНВ. В противном случае внесенная популяция полностью вымирает. В частности, если для начального распределения выполняется , то популяция вымирает при любом значении протяженности ЛНВ. При наличии хищника, в том случае, когда пресс со стороны популяции v достаточно велик (например, при высокой пищевой активности хищника или при его большой концентрации, что означает большие значения функции Е), плотность популяции u может “сбрасываться” на уровень ниже порогового, и это означает полное вымирание жертвы. При этом, если в начальный момент времени плотность хищника невелика, его влияние на популяцию жертвы на ранней стадии распространения популяционной волны может практически не ощущаться. Таким образом, оказывается возможным эффект динамического удержания ЛНВ: возникающая (при превышении критических размеров) популяционная волна через определенное время оказывается полностью “съеденной” увеличившейся популяцией хищника. Еще раз подчеркнем, что затухание локальной вспышки численности вследствие влияния хищника возможно не во всяком сообществе “хищник - жертва”, а только когда популяция жертвы имеет критический порог плотности (то есть является популяцией типа Олли). Более того, и в этом случае эффект динамического удержания реализуется не всегда, а лишь при выполнении определенных условий: определенная ширина начального распределения популяции хищника, определенная для данной ширины амплитуда начального распределения популяция хищника. Отметим, что явление затухания начального возмущения вследствие динамического удержания принципиально отличается от хорошо известного “диффузионного” затухания. В случае динамического удержания сначала формируется волна, которая распространяется от области ЛНВ на значительные (по сравнению с размерами начального возмущения) расстояния с почти постоянной амплитудой, после чего происходит ее быстрое схлопывание. В случае же диффузионного затухания волна не образуется, и размер возмущенной области монотонно убывает со временем.

С помощью минимальных моделей исследуются и другие проблемы. Так, например, в работе (Jang, Allen, 1999) рассмотрена динамика простой пищевой цепи, состоящей из трех трофических уровней: (1) биогенного элемента или субстрата N, (2) потребителя питательных веществ или жертвы P, (3) хищника Z. Такая обобщенная форма пищевой цепи в данной работе конкретизируется в двух типах модельных систем: (а) планктонной системе (биогены - фитопланктон - зоопланктон) и (б) хемостатной системе (биогены - бактерии - простейшие). Известно, что ряд форм биогенных элементов, таких, например, как нитраты, нитриты и аммонийный азот могут лимитировать рост организмов при низких концентрациях и ингибировать рост при высоких своих концентрациях. Нитрат, например, лимитирует рост организмов при низких концентрациях, однако оказывается токсичным при высоких концентрациях. Такого рода ингибирующее воздействие может учитываться, например, введением дополнительного множителя в скорости роста. Нитритный и аммонийный азот при высоких концентрациях являются ингибиторами роста некоторых бактерий. Продукты промышленного загрязнения, такие как фенолы и тиоцианаты, часто подавляют процессы метаболизма у бактерий. Воздействие подобного рода ингибиторов роста на динамику простой пищевой цепи и является предметом исследования в данной работе. Особенностью данной модели является обобщенное рассмотрение скорости потребления биогенов и скорости выедания. Кроме того, скорость поглощения биогенов не имеет обычно используемой монотонной формы кинетики Михаэлиса - Ментен - Моно. Вместо этого рассматривается вариант немонотонной скорости поглощения. При исследовании динамики хемостатной системы (также как и для планктонной системы) предполагается, что субстрат лимитирует рост организмов при низких уровнях и ингибирует его при высоких концентрациях.

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

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

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

(48)

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

(49)

Эта система, как показано в данной работе, является диссипативной.

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

(50)

Здесь D - скорость разбавления, - доля субстрата, трансформированного (обращенного) в биомассу жертвы, - доля субстрата, трансформированного в биомассу хищника, - максимальные скорости роста жертвы и хищника соответственно. Смысл функций и такой же, как и в планктонной модели. Отличительной особенностью хемостатной модели является то, что в ней нет рециклинга субстрата. Таким образом, в общем случае условие не выполняется. Тем не менее, если скорость разбавления D постоянная, то для величины справедливо соотношение . В этом случае анализ трехмерной системы можно свести к анализу двумерной системы относительно P и Z при условии, что . Если обозначить , , то двумерную систему для хемостата можно записать в следующем виде:

(51)

Предполагается, что функции и удовлетворяют следующим условиям:

(i) . Существует такое , что для и для . Существует также такое , что для и для .

(ii) и .

Если , то функция f является монотонной. Однако, если , то функция f представляет собой выпуклую функцию с одним максимумом. Ингибирование имеет место при концентрациях . В данной работе концентрация определяется как токсический уровень, то есть это та концентрация биогенного элемента или субстрата, при котором имеет место ингибирование процесса роста и развития организмов. Таким образом, при исследовании динамики планктонной и хемостатной систем необходимо рассматривать три варианта условий: (I) , (II) , (III) . Величины порогов определяются следующим образом: , .

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

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

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

Предполагается, что рост фитопланктонных организмов, способных вырабатывать токсины, описывается логистической зависимостью со скоростью роста r. Математическая модель, учитывающая рост фито- и зоопланктона, их трофические взаимодействия, естественную смертность зоопланктона и негативное влияние на его развитие со стороны фитопланктона, имеет следующий общий вид (Chattopadhayay et al., 2002):

(52)

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

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

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

. (53)

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

, (54)

в то время как второе уравнение, после сокращения на Z, может быть записано в следующем виде:

. (55)

Из уравнения (55) можно найти Р. Следует заметить, что если особая точка является устойчивой особой точкой, то уравнение (55), а следовательно и уравнение (54), не имеет решений. Таким образом, для существования решений необходимо потребовать, чтобы эта точка была неустойчивой и выполнялось неравенство . Из (55) следует, что . Поскольку , то это означает, что и значит

. (56)

Рассматривая уравнение (54), можно заметить, что правая его часть зависит от выражения , которое не является монотонной функцией. Нетрудно видеть, что на интервале эта функция имеет максимум, который достигается при .

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

. (57)

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

. (58)

Следовательно, очевидно, что критерий устойчивости особой точки системы (52) будет целиком зависеть от множителя .

Поведение решений системы (52) существенным образом зависит от вида функций и . Возможны следующие девять случаев:

·Случай 1. Обе функции и являются линейными, то есть .

·Случай 2. Обе функции и являются функциями Холлинга II-го типа, то есть .

·Случай 3. Обе функции и являются функциями Холлинга III-го типа, то есть .

·Случай 4. является линейной функцией, а представляет собой функцию Холлинга II-го типа, то есть , .

·Случай 5. представляет собой функцию Холлинга II-го типа, а является линейной функцией, то есть , .

·Случай 6. является линейной функцией, а представляет собой функцию Холлинга III-го типа, то есть , .

·Случай 7. представляет собой функцию Холлинга III-го типа, а является линейной функцией, то есть , .

·Случай 8. представляет собой функцию Холлинга II-го типа, а является функцией Холлинга III-го типа, то есть, , .

·Случай 9. представляет собой функцию Холлинга III-го типа, а является функцией Холлинга II-го типа, то есть, , .

Во всех этих случаях параметр представляет собой константу полунасыщения.

Поскольку выделение токсина снижает скорость роста зоопланктона, вызывая существенную его смертность, то в этот период фитопланктон оказывается менее доступным для выедания. Кроме того, следует учитывать особенности трофических взаимоотношений между фито- и зоопланктоном при насыщении зоопланктона пищей и при недостаточном обилии фитопланктона. В эти периоды тактики выживания и добывания пищи существенно различаются. Для описания подобного рода явлений использование функций Холлинга II-го и III-го типов наиболее адекватно. Рассмотрение случаев 1) - 9) позволяет изучить различные механизмы и сценарии возможного возникновения цветения фитопланктона и его контроля.

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

В остальных случаях при определенных значениях параметров возможно появление периодических решений, которые приводят к возникновению цветения фитопланктона. Механизм прекращения (или торможения) развития фитопланктона пока не совсем ясен. Возможно, что его останавливают некоторые химические вещества, которые выделяются самим фитопланктоном. Имеются предположения о том, что в этом процессе могут участвовать некоторые вирусы (Beltrami, Carroll, 1994).

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

Одним из современных и бурно развивающихся направлений в изучении и качественном моделировании распределенных химических и экологических систем является использование клеточных автоматов, как детерминированных, так и случайных (Ванаг, 1999). Одно из основных отличий клеточных автоматов (КА) от системы обыкновенных дифференциальных уравнений (ОДУ) заключается в локальности правил, с помощью которых описывается динамика системы. В случае применения ОДУ мы пользуемся некоторыми правилами изменения усредненных по всей системе величин например, средних концентраций. При этом мы a priori полагаем, что такие правила существуют. В случае КА существование таких обобщенных правил (или макроправил) необязательно. Достаточно знать законы развития системы на микро- или мезоуровне в небольших пространственных областях (ячейках), из которых состоит макросистема. Важно лишь, что эти локальные правила одинаковы для всех ячеек. Другим важным отличием КА от ДУ является использование не только дискретных, но и, как правило, целочисленных переменных. Дискретность переменных позволяет рассматривать большой класс разрывных недифференцируемых функций. Надо отметить, что дискретные свойства КА заметно уменьшаются при работе с большими значениями переменных, но никогда не исчезают совсем. Всегда существует минимальный дискретный шаг изменения переменной. В случае же численного решения ОДУ или ДУ в частных производных можно уменьшать шаг дискретности до сколь угодно малых величин. Локальность законов поведения и дискретность переменных позволяют при использовании КА естественным образом учитывать флуктуации или внутренние шумы системы без каких-либо дополнительных предположений. Для учета флуктуаций в случае применения ДУ используются стохастические ДУ, для решения которых необходимо привлекать различные флуктуационно-диссипативные соотношения. Во многих случаях вывод таких соотношений является весьма трудоемкой задачей. Совокупность таких свойств, как локальность правил поведения и дискретность переменных порождает порой совершенно непредсказуемое поведение распределенной дискретной системы. Непредсказуемость поведения (не путать с хаотичностью!) является еще одной чертой КА и в данном контексте носит скорее качественный характер.

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

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

ВКА представляет собой регулярную решетку, состоящую из элементарных ячеек. Форма решетки может быть не только квадратной, но и прямоугольной с сильно различающимися длинами сторон. Каждая ячейка характеризуется некоторым набором целых чисел: числом молекул соответствующего сорта в данной ячейке и своими координатами (i и j). Ячейке приписывается также определенный объем и линейный размер . Объем используется при задании вероятностей протекания химических реакций в ячейках. Все ячейки считаются гомогенными. Это условие накладывает ограничения как на характерный размер, так и на характерное время: , где - характерное время химической реакции, а - коэффициент молекулярной диффузии. Размер ячеек не должен также превышать размера самых маленьких турбулентных вихрей. Следовательно, от не должен превышать размера Колмогорова , который представляет собой внутренний масштаб турбулентности, определяющий расстояния, на которых вязкость жидкости начинает играть заметную роль. Оценки показывают, что l не должно превышать 1 - 10 мкм.

Методом ВКА независимо моделируются три процесса: молекулярная диффузия, химические реакции, турбулентное перемешивание. При моделировании диффузии и турбулентного перемешивания используются, как правило, периодические граничные условия, то есть клетки одного края решетки считаются в контакте с клетками противоположного края. Рассмотрим блоки “Диффузия” и “Перемешивание” более подробно.

В блоке “Диффузия” предполагается, что массообмен между двумя соседними ячейками происходит следующим образом: на короткое время ячейки сливаются, за время слияния их содержимое полностью перемешивается, а затем пара распадается на две равновеликие ячейки прежнего размера. Для такого механизма вероятность того, что в момент времени в двух произвольных соседних клетках будет m и k молекул одного сорта, если в момент времени t в них находилось соответственно r и s молекул, где , а - некоторый постоянный момент времени, соответствующий одному временному шагу при моделировании диффузии, задается биномиальным распределением , где q - отношение объема ячейки к сумме объемов двух слившихся ячеек (в данном случае равновеликих ячеек, поэтому q = 1/2). Подчеркнем, что интервал времени является реальным временем, которое измеряется в секундах. Соседними в данной модели называются восемь клеток, окружающих центральную клетку. Нахождение чисел m и k производится следующим образом. Вначале случайным образом выбирается пара соседних ячеек и для нее определяется сумма r + s. Далее генерируется случайное число и ищется такое k, которое удовлетворяет неравенствам . Интенсивность массообмена регулируется величиной и параметром ND, определяющим число случайно выбранных пар соседних клеток за один временной шаг . Массообмен осуществляется независимо для каждого сорта молекул. Этому соответствуют различные числа ND. Как правило, оно не превышает 0,5. Отношение определяет константу (или частоту) перескока частицы в произвольную соседнюю клетку. Для гомогенных водных растворов величины , (или l) связаны через коэффициент диффузии : . В тех же случаях, когда нет прямой связи между , и , величины и могут быть независимыми переменными. Справедливость предложенного закона массообмена была проверена на примере различных задач математической физики.

Вблоке “Перемешивание” турбулентное движение рассматривается как результат наложения турбулентных пульсаций различных масштабов. При малых числах Рейнольдса Re в системе существуют только крупномасштабные пульсации, соизмеримые с размером всей системы. По мере возрастания Re появляются все более и более мелкие пульсации. Для моделирования крупномасштабных пульсаций на базовой решетке размера (например, 2⁸´2⁸, ) случайным образом выбираются квадраты со стороной , где n = 7, а затем четыре квадранта каждого выбранного квадрата случайным образом переставляются местами. Увеличение интенсивности перемешивания моделируется последовательным добавлением все более мелких масштабов турбулентного движения, то есть дополнительным выбором квадратов со стороной , где n = 6, 5, 4, 3, 2, 1 и соответствующей перестановкой составляющих их квадрантов. Числа выборок квадратов со стороной в единицу времени были подобраны так, чтобы скорость перемещения элемента жидкости для каждого масштаба движения удовлетворяла закону Колмогорова-Обухова. Было получено, что . Таким образом, в модели учитывается иерархия масштабов турбулентного движения. Этот блок запускается через равномерные интервалы времени , причем , как и измеряется в секундах, то есть процесс счета проводится в реальном времени. Использовался также несколько видоизмененных алгоритм блока “Перемешивание”. Сторона квадрата L , выбираемого в произвольном месте решетки, случайным образом может принимать одно из допустимых по программе значений. Интенсивность перемешивания регулируется временным шагом и масштабом наименьшего из выбираемых квадратов.

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

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

Рассмотрим систему уравнений в форме Лотка-Вольтерра:

, (59)

где - численность (плотность) i-го вида, - удельная скорость роста, - коэффициент, определяющий влияние j-го вида на i-ый. Мэй (May, 1972; May, 1974) исследовал статистическое распределение собственных чисел для случайной матрицы сообщества. Предполагалось, что виды подвержены самолимитированию, так что значения диагональных коэффициентов влияния равны -1. Некоторым недиагональным элементам (их доля равна С) приписываются случайные значения с нулевым средним и стандартным отклонением s. Знак коэффициента влияния определяет тип взаимодействия между видами. Отрицательный знак соответствует хищничеству или паразитизму, а положительный - мутуализму. Если виды не оказывают друг на друга никакого влияния, то соответствующим коэффициентам в матрице сообщества приписывалось нулевое значение. Таким образом, определенному множеству коэффициентов соответствует специфичная структура сообщества. Мэй высказал предположение о том, что сообщество будет неустойчивым, если выполняется неравенство:

, (60)

где m - число видов в системе. В соответствии с этим критерием устойчивости, устойчивость случайной системы Лотка-Вольтерра будет уменьшаться с увеличением числа видов m, количеством связей в системе C и возрастанием силы взаимодействия s. Позднее Коэн и Ньюмен (Cohen, Newman, 1985) показали, что критерий Мэя в целом не корректен, а справедлив только для достаточно узкого класса моделей пищевых цепей. Тем не менее, обратная зависимость между числом видов и связностью матрицы сообщества, подразумеваемая в соотношении (60), подмечена правильно (Cohen, Newman, 1988; Cohen et al., 1990).

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

. (61)

При описании взаимодействия среди отдельных организмов предполагалось, что в любой момент времени данную ячейку решетки способна занимать только одна особь. Взаимодействия между видами имели место только в том случае, если они в данный момент времени стремились занять одну и ту же ячейку. Исход этих встреч определялся совокупностью коэффициентов, аналогичных коэффициентам влияния в системе (59). Моделирование проводилось по двум вариантам: в режиме “локального взаимодействия” и в режиме “среднего поля”. В первом случае особи перемещались с равной вероятностью в одну из четырех ближайших соседних клеток. При этом задавались отражающие граничные условия. В режиме “среднего поля” особи перемещались с равной вероятностью в любую клетку решетки. Обновление ячеек решетки было асинхронным, то есть результат перемещения каждой особи или результат взаимодействия мгновенно изменял состояние ячейки прежде чем другая особь начинала свое движение. Все ячейки решетки посещались в случайном порядке в течение каждого этапа обновления. Когда особь i-го вида пыталась переместиться в ячейку, занятую особью j-го вида, результат взаимодействия зависел от набора коэффициентов взаимного влияния. Эти коэффициенты ранжировались от -1.0 до 1.0, причем полагалось, что в общем случае . С помощью этих коэффициентов рассчитывались две вероятности. Первая их них задавала вероятность того, что особь i-го вида будет двигаться к ячейке, занятой особью j-го вида. Эта вероятность перемещения определяется с помощью соотношения

. (62)

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

(63)

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

Моделирование проводилось как для однородной решетки, на которой все ячейки были доступны для заселения, так и для неоднородной решетки, на которой некоторые ячейки были организмам недоступны. В первом случае исследовалась возможность возникновения пространственных структур, обусловленная только межвидовыми взаимодействиями, сосуществованием видов и сложностью пищевой цепи. Было показано, что локальные взаимодействия между видами приводят к целому множеству пространственных структур, среди которых во многих случаях удается обнаружить структуры типа бегущих волн, состоящих из трех и более видов. Иногда формировались устойчивые участки, формируемые каким-либо одним-единственным видом. При моделировании, использующим режим “среднего поля”, показано, что между численностью видов, а также между силой взаимодействия и связностью существует значимая степенная зависимость. Тем не менее, регрессионные параметры значительно отличались от значений, ожидаемых согласно критерию устойчивости Мея. Наиболее явное отличие заключалось в слабой зависимости численности видов от силы взаимодействия, в то же время связность оказывала более значительное влияние на сосуществование видов. Различие же между критерием Мея и моделью в единицах устойчивого числа видов было относительно невелико. При моделировании, использующим режим “локального взаимодействия” без смертности, в качественном отношении выявлены аналогичные закономерности, хотя количество устойчиво сосуществующих видов в этом случае было наиболее высоким. Добавление в эту модель 5%-ого уровня смертности значительно ослабляло воздействие связности на видовое обилие, в то же время, как и ранее, сила взаимодействия оказывала относительно небольшое влияние на количество сосуществующих видов. Введение в модель пространственной неоднородности приводит к существенно нелинейной зависимости видового обилия от степени гетерогенности среды (решетки). Количество устойчиво сосуществующих видов сначала растет с ростом пространственной неоднородности и достигает максимального своего значения при 70%-ом уровне свободных ячеек, доступных для заселения. После чего число видов начинает уменьшаться с увеличением общего количества недоступных ячеек. Воздействие связности на число видов также меняется в зависимости от уровня доступных для заселения ячеек. Показано, что между связностью и степенью пространственной гетерогенности среды существует достоверная зависимость. Общий вывод заключается в том, что пространственно структурированные пищевые цепи оказываются более устойчивыми в отношении количества сосуществующих видов по сравнению с пищевыми цепями без пространственной структуры.

Большинство современных моделей водных экосистем сформулированы в терминах агрегированных, а не индивидуальных, переменных. В сущности, такие модели исследуют временную динамику, если так можно выразиться, “усредненных” по тому или иному признаку индивидуумов. Например, в простых качественных моделях, которые мы рассматривали выше, нас не интересовала размерная, половая и возрастная структура зоопланктонного сообщества, принадлежность организмов к разным видам. Вместе с тем, в последнее время получает развитие совершенно особый класс моделей водных экосистем - индивидуально-ориентированные модели (или IBM-модели) (Individual-Based Models) (DeAngelis et al., 1991; Crowder et al., 1993; Rose et al., 1993; Hinckley et al., 1996; Werner et al., 1996; Hermann et al., 2001). В таких моделях описание целостного поведения системы начинается на уровне отдельно взятого организма. При этом учитывается поведение и взаимодействие каждого индивидуума с другими организмами и с окружающей средой, текущее его состояние, а возможно, и ряд состояний в прошлом. Двигаясь “снизу - вверх”, мы в идеале должны понять, каким образом из частных, отдельных взаимодействий возникают целостные, системные свойства (Grimm, 1999). Эта задача, конечно, не может быть решена без движения “сверху - вниз”, то есть без понимания того, что такое системные свойства и какие причины приводят к их возникновению, то есть как на самом деле происходит самоорганизация системы.

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

Индивидуально-ориентированные модели можно использовать совместно с другими блоками моделей экосистем. Следует лишь заметить, что традиционное построение моделей (особенно термогидродинамических блоков и блоков, описывающих адвективный и диффузионный перенос компонентов экосистемы) и их численная реализация основываются на Эйлеровом описании, в то время как IBM-модели требуют Лагранжева подхода (Поддубный, Сухова, 2000; Поддубный, Сухова, 2002; Поддубный, Сухова, 2002; Hermann et al., 2001).

Дополнительная информация:

печатные издания и Интернет-сайты (1324)

В начало

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

Реестр моделей

Модели в экологии

Модели водных экосистем

Основные типы математических моделей водных экосистем

Качественные модели водных экологических систем