Моделирование динамики численности планктоновых сообществ с помощью модели «Хищник-жертва» с учетом эффекта ограниченного перемещения хищника Текст научной статьи по специальности «Математика»

электронное научно-техническое издание

НАУКА и ОБРАЗОВАНИЕ

Эл № ФС 77 - 30569. Государственная регистрация №0421100025. ISSN 1994-040S

Моделирование динамики численности планктоновых сообществ с помощью модели «хищник-жертва» с учетом эффекта ограниченного перемещения хищника.

77-30569/367901 # 04, апрель 2012

Кренке А. Н., Чернышев В. Л., Коломенский Д. С.

УДК. 574.583+574.34+51-76

Институт географии РАН К.гепке-10И.А8@уапёех. ги [email protected] [email protected]

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

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

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

Зоопланктон - это планктонные животные. Зоопланктон представлен как травоядными, так и хищниками; при этом травоядные питаются фитопланктоном и являются пищей для зоопланктонных хищников. Как основа трофической цепи система «фито- и зоопланктон» и ее динамика существенно влияют на все трофические связи и сети в океане.

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

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

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

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

Возникает вопрос: могут ли биологические факторы, такие, как рост и взаимодействия в системе «хищник-жертва», быть причиной формирования планктонных структур без каких бы то ни было гидродинамических воздействий?

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

Методы

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

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

В работе [6] предложен способ описания динамики системы двух взаимодействующих популяций планктона с эффектом взаимного нелокального ингибирования по времени и пространству посредством свертки функций ингибирования и аллелопатии. Такой же метод можно применить и для построения модели типа «хищник-жертва» [3], реализующей взаимодействие фито и зоопланктона на локальном уровне, при предположении условной однородности среды, то есть моделируется ограниченное «рысканье» массы зоопланктона.

В данной работе рассмотрена следующая постановка задачи. Пусть ЩхД) -пространственно-временное распределение ресурса (фитопланктона), а ^"(хД) -пространственно-временное распределение потребителя (зоопланктона). Взаимодействие видов происходит пропорционально количеству N согласно некой трофической функции У^). Кроме того, введем в рассмотрение функцию «активности» или «дальнодействия» потребителя Р(|х-^|), которая показывает, как убывает по пространству потребление данного ресурса потребителем, расположенным в точке х (рис. 1). Необходимым условием для функции Р(|х-^|) является то, что она стремится к нулю по ходу удаления от точки х. Также необходимо определить естественный прирост ресурса 0^), смертность потребителя как (-М^ и коэффициент переработки ресурса в потребителя (К).

Рис. 1. Распределение ресурса (синий), потребителя (красный) по ареалу (X) и функция

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

В соответствии с полевыми наблюдениями за динамикой планктонных систем [18], локальный рост популяции «жертвы» или ресурса является логистическим, а трофическая функция хищника принадлежит второму Холлинговскому типу, то есть имеется некий порог насыщения хищника [1]. Таким образом, используя безразмерные параметры, получим

Q(R) = R(1-R), V(R) = R/(R+h), где h - параметр насыщения хищника.

Теперь рассмотрим уравнение динамики планктоновой популяции с учетом пространственного взаимодействия. Для этого осуществим свертку функции P(|x-^|) и N(x,t) в уравнении, определяющем динамику потребителя, и свертку P(|x-^|) и V(R) уравнении, определяющем динамику ресурса.

Или:

Q(R) — V(R(x, t) J N(£, t)P(|x —

—да

+» (1)

= —MN + к ■ N(x, t) J V(£, t)P(| x — £ |)d£

dR

Yt

8N_ dt

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

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

дх

¿>(£ 0*41*1 - = + ^ ПГ|*1 - к

дн г+оо ,

>41',- + (вд О + V С ^ «I* - (+

(3)

Предположим, что функция P(|x-ф симметрична относительно х и

=0,

тогда можно приравнять нулю первый интеграл в (2-3). Если обозначить через о2 второй момент (дисперсию) функции Р(|х-^|) , то интеграл / \х — %|2Р(|х — £= а2, что следует из определения дисперсии для непрерывных значений. Отбрасывая все члены выше третьего порядка, получим следующую систему уравнений в частных производных:

дЯ

дЫ

= 0 ( Я ) - V ( Я ( х, ? ))( N ( X, t ) + д

2 дх _2

= -^ + К ■ N(х,t){(V(х,t) +^УГ(^)2 + V г

,д 2 К

(4)

]}

дt 2 и ""ч дхх г дх2

Подобные системы называют «шредингеровскими», теория этих систем относительно слабо развита, и они могут обладать весьма интересными свойствами. В частности, здесь могут существовать периодические по пространству решения и решения с «мозаичной» пространственной структурой, где из начального всюду плотного распределения возникает устойчивая структура чередований наличия и отсутствия потребителя и ресурса [4]. Дополним задачу граничными условиями. В данной задаче естественно выбрать граничные условия периодического типа, так как рассматриваемый здесь ареал в 5-10 км является лишь частью значительно большего ареала распространения планктона.

Для численного решения полученной системы уравнений применяется метод конечных разностей. Выбрана частично неявная схема. Такой выбор обусловлен поиском компромисса между скоростью вычислений и предполагаемой устойчивостью схемы [2, 9]. Временная координата обозначена ^ а пространственная х; участок пространства-времени, на котором выполняется численное решение, подразумевается конечным.

дс

£Т2 Я&^-МГР+ЧН^

Ах2

(5)

-м(^)-к-м*- +£ (^т^)2+к

гк'

(6)

Где 1 =1 ,.,пх- шаг схемы по пространственной координате, пх - количество точек разбиения по х, а к=0,..., пь — 1 - шаг схемы по времени, где п1: количество точек разбиения по t.

Решение системы уравнений (9)-(10) в некой точке 1 относительно будет выглядеть как:

(7)

где матрица коэффициентов Л имеет вид:

Ас

а вектор правой части выглядит как:

Решая уравнение (7), для всех точек разбиения по х необходимо решить систему из 2пх уравнений. Удобно записать матрицу этой системы через блочные матрицы, в которых верхние пх строк заполнены коэффициентами ¿-того уравнения (5), а нижние пх строк - коэффициентами для ¿-того уравнения (6). Блочная матрица будет содержать четыре блока: Ак - матрица коэффициентов при в уравнении 5, Вк - матрица

коэффициентов при в уравнении 5, Ск - матрица коэффициентов при в

уравнении (7), Ок - коэффициентов при в уравнении (7). В соответствии с видом

матриц в системе (7), матрицы Ак и Вк будут диагональными, а Вки Ск -трехдиагональными. Учитывая периодические граничные условия, последний элемент первой строки и первый элемент последний строки матриц Вки Ск равны, соответственно, коэффициенту при и и Д^1 и

Решение данной системы линейных уравнений относительно векторов

N

к+1

сводится к решению следующего матричного уравнения вида:

(ак Вк\(Як+1Л = (Рк\

(8)

где

■-Д = ■;

;

= л^ + Агк^

- ДtM^V/г;

/ 1,..., пх, к=0,..., щ - 1.

Решение системы уравнений (8) получается элементарным образом:

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

Сходимость разностной схемы оценивается через невязку относительно наиболее точного решения с количеством узлов равным 212 по пространственной координате и 214 по временной. Невязка оценивается по норме Ь2 и 1/'. Оценка невязки осуществляется отдельно относительно шага по времени и отдельно относительно шага по пространству, таким образом, получается две оценки. Шаг уменьшения точности решения в обоих случаях соответствует 2-1, минимальная точность для пространственного и временного разбиения составляет 23 и 25 соответственно. Из результатов проверок сходимости можно сделать вывод, что невязка относительно наиболее точного решения монотонно уменьшается при увеличении точности проверочных решений. Таким образом, можно предположить применимость предложенной вычислительной схемы для решения системы

(4).

Результаты

Исходя из постановки задачи, рассматривается ареал популяции планктона, условно гомогенный по свойствам среды. Следует отметить, что большие амплитуды изменчивости плотности планктона в мезомасштабе (1-100 км) и микромасштабе (100 м -10 км) часто наблюдаются в отсутствие значительных физических воздействий на него. В данной работе рассмотрен масштаб 10 км. Что касается временного интервала, то разумно выбрать его похожим на время теоретически возможного существования структур, порожденных внутренними процессами системы фито-зоопланктон. В океане значительные изменения среды за счет гидрологических процессов происходят раз в несколько дней даже на спокойных участках. Таким образом, временной интервал в 30 дней более чем достаточен.

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

Коэффициенты h=H/b, т=М/а и к=Ки/а должны быть подобраны таким образом, чтобы, с одной стороны, происходило сосуществование фито и зоопланктона, а с другой, можно было бы сделать вывод об эффекте аллелопатии на динамику системы. Если

рассмотреть локальную динамику системы, то есть без эффекта аллелопатии, то легко показать, что система

^ = - Ю--— N = 0

а R + h

йЫ R

-= -тЫ + к-N = 0

а R + h

имеет три стационарных состояния (0,0) - полное вымирание, (1,0) - вымирание хищника и стабилизация за счет самоингибирования жертвы, и сосуществование хищника и жертвы, при котором

—

гк 1 —г

= (1 - ЯД

где г =m/k.

Отметим, что состояние (Д*, /V*), при котором виды сосуществуют и в биологически

1—г

значимой области Я>0 и N>0, соответствует значению параметра К < —. Это легко проверить подстановкой.

На основании вышеизложенного, в качестве параметров можно взять 11=0.4, г=0.3, к=2, учитывая соотношение г =ш/к, ш=0.6 [7,1]. Дисперсию определим как а2 = 0.22, стандартное отклонение, равное 0.2. Пространственно временная динамика системы хищник жертва, очевидно, в значительной степени определяется выбором начальных условий. Простейшей формой распределения можно считать пространственно однородные начальные условия. Однако в этом случае распределение видов будет однородным в любой момент времени. Пространственные структуры в ходе развития событий появиться не могут. В данной работе нас интересует, как разные формы неоднородности начального распределения ресурса и жертвы, вызванные экзогенными факторами, могут повлиять на развитие динамики системы. Рассмотрим случай, в котором на некоторой области ареала присутствует начальное возмущение, например след от прошедшей стаи китов, выевшей планктон в данной области. Выберем в качестве начального распределения следующие функции

Д, = 0.4 -

0.1

1+10 О;-ЗУ

= 0.3-

0.1

1+10(л:;-5)2

(9)

На рисунке 2 приведены начальные условия, заданные этими функциями.

Рис. 2. Распределение, заданное уравнениями (9).

Рассмотрим поведение системы при данных начальных условиях.

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

Необходимо отметить, что, по всей видимости, одним из эффектов аллелопатии является взаимное перемещение пиков плотностей ресурса и потребителя, т.е. взимодействие между ними носит не просто пропорциональный или обратно пропорциональный характер, а происходит блуждание потребителя, который премещается вдоль наиболее «изобильных» участков ареала. В соответсвии с этим происходит распространиение волн в пространстве относительно центра возмущения. В целом можно сказать, что волнообразная структура распределения плотности планктона с эффектом затухания колебаний по ходу удаления от эпицентра довольно часто встречается в природе. Такие структуры, вероятно, можно объяснить внезапным изменением численности планктона на неком небольшем участке.

Рис. 3. Распределение плотностей ресурса и потребителя в различные моменты времени для начальных условий, заданных функциями (9).

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

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

В данном случае выберем в качестве начального условия функции следующего

вида:

Д, = 0.4 -

ЛГ( = 0.3 -

0.15

0.1

0.1

1+10(зе;-3)2 1+50;-б)2 1+10(>;-9)2' 0.15 0.1 0.1

1 + 10(л:;-3)2 1+5(>;-б)2 1+10 О;-8)2

На рисунке 4 приведен график соответствующего начального распределения.

(10)

Рис. 4: распределение, заданное уравнениями (10)

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

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

Рис. 5. Распределение плотностей ресурса и потребителя в различные моменты времени для начальных условий, заданных функциями (10).

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

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

Ниже приведена формула вычисления кросс-корреляции двух рядов:

Здесь п - смещение относительно элемента ^ а М - мат. ожидание, вычисляемое как среднее ряда.

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

Для случая с несколькими источниками возмущения кросс-корреляции Я и К, вычисленные с лагом вплоть до 1/2 исследуемого ареала в момент времени t=25, не превосходят |0.54|, наблюдающегося при смещении одного ряда относительно другого на ~1/3. В случае с одним источником возмущения максимальная кросс-корреляция для t=25 составляет |0.68|.

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

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

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

Заключение

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

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

1,1 [(я(0 - Мй) * (N(1-71)- МАО]

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

Работа выполнена при частичной финансовой поддержке грантов МК-4255.2012.1, РФФИ 10-07-00617-а и 12-07-00267-а, а также программы «Научные и научно-педагогические кадры инновационной России» (госконтракт 14.740.11.0794).

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

1. Медвинский А.Б., Петровский С.В., Тихонова И.А., Тихонов Д.А., Ли в Б. Л., Вентурино Э., Мальхё Х., Иваницкий Г.Р. Формирование пространственно-временных структур, фракталы и хаос в концептуальных экологических моделях на примере динамики взаимодействующих популяций планктона и рыбы // Успехи физических наук, 2002. Т. 172, №1. c.31-66.

2. Самарский А.А. Введение в теорию разностных схем.- М.: Наука, 1971.- 552 с.

3. Свирежев Ю.М. Нелинейные волны, диссипативные структуры и катастрофы в экологии.- М.: Наука, 1987. - 366 с.

4. Хакен Г. Синергетика: пер с англ.- М.: Мир, 1980. - 406 с.

5. Levin S.A., Segel L.A. Hypothesis for origin of planktonic patchiness// Nature, 1976. № 259, c. 659-659.

6. Mukhopadhyay A., Tapaswi P.K. A space-time state-space model of phytoplankton allelopathy// J.Nonlinear Analysis:Real World Applications, 2003. № 4, p. 437-456.

7. Mimura M, Murray J.D. On a diffusive prey-predator model which exhibits patchiness// J. Theor. Biol., 1978. №3. p. 249-262.

8. Segel L.A., Jackson J.L. Instability of a layer of chemotactic cells, attractant and degrading enzyme// J. Theor. Biol., 1972. №3 p. 399-415.

9. Steele J.H., Henderson E.W. The role of predation in plankton models// J Plankton Res, 1992. №1. p. 157-172.

10. Scheffer M .J. Should we expect strange attractors behind plankton dynamics - and if so, should we bother? // J. Plankton Res., 1991. № 6. p. 1291-1305.

electronic scientific and technical periodical

SCIENCE and EDUCATION

_EL № KS 77 -3()56'J..VaU421100025. ISSN 1994-jMOg_

Simulation of the plankton's population dynamics using a predator-prey model with weak mobility of the predator.

77-30569/367901 # 04, April 2012

Krenke A.N., Chernyshev V.L., Kolomenskii D.S.

Institute of Geography Krenke-IGRAS@yandex. ru [email protected] [email protected]

The paper considers the possibility of organizing the spatial patterns of plankton's population densities defined solely by biological factors, under the condition of homogeneous environment. To simulate self-organization of two plankton's species populations a mathematical model of the "predator-prey" type with regard to the effect of limited seeking by the predator was studied. The simulation results showed the fundamental possibility of spatial structures organization through limited self-movement of plankton.

Publications with keywords: population dynamics, trophic model, plankton, Schrodinger-type systems

Publications with words: population dynamics, trophic model, plankton, Schrodinger-type systems

References

1. Medvinskii A.B., Petrovskii S.V., Tikhonova I. A., Tikhonov D.F., Li v B. L., Venturino E., Mal'khe Kh., Ivanitskii G. R. Formirovanie prostranstvenno-vremennykh struktur, fraktaly i khaos v kontseptual'nykh ekologicheskikh modeliakh na primere dinamiki vzaimodeistvuiushchikh populiatsii planktona i ryby [The formation of spatial-temporal patterns, fractals and chaos in conceptual ecological models on the example of the dynamics of interacting populations of plankton and fish]. Uspekhi fizicheskikh nauk, 2002, vol. 172, no. 1, pp. 31-66.

2. Samarskii A.A. Vvedenie v teoriiu raznostnykh skhem [Introduction to the theory of difference schemes]. Moscow, Nauka, 1971. 552 p.

3. Svirezhev Iu.M. Nelineinye volny, dissipativnye struktury i katastrofy v ekologii [Nonlinear waves, dissipative structures and catastrophes in ecology]. Moscow, Nauka. 1987. 366 p.

4. Haken H. Sinergetics: An Introduction. Nonequilibrium Phase-Transitions and Self-Organization in Physics, Chemistry and Biolog. lrd ed., Berlin, Springer, 1977. (Russ. ed.: Khaken G. Sinergetika, 1980, Moscow, Mir. 406 p.).

5. Levin S.A., Segel L.A. Hypothesis for origin of planktonic patchiness. Nature, 1976, no. 259, pp. 659-659.

6. Mukhopadhyay A., Tapaswi P.K. A space-time state-space model of phytoplankton allelopathy. J. Nonlinear Analysis: Real World Applications, 2003, no. 4, pp. 437-456.

7. Mimura M, Murray J.D. On a diffusive prey-predator model which exhibits patchiness. J. Theor. Biol, 1978, no. 3, pp. 249-262.

8. Segel L.A., Jackson J.L. Instability of a layer of chemotactic cells, attractant and degrading enzyme. J. Theor. Biol, 1972, no. 3, pp. 399-415.

9. Steele J.H., Henderson E.W. The role of predation in plankton models. J. Plankton Res, 1992, no. 1, pp.157-172.

10. Scheffer M .J. Should we expect strange attractors behind plankton dynamics -and if so, should we bother?. J. Plankton Res, 1991, no. 6, pp.1291-1305.