Алексеев В.В. Физическое и математическое моделирование экосистем

Подождите немного. Документ загружается.

На

плоскости параметров

(С, U)

изоклины

О = 0

имеют

вид

прямых 1—4 (рис. 7.7), изоклина С

= 0

имеет вид кривой

ABMN.

В положениях

1 и 4

стационарные точки системы асимптотиче-

ски

устойчивы.

положении

3 на

участке возрастания изоклины

Рис.

7.6.

Зависимость скорости

про-

изводства цАМФ

от

концентрации

цАМФ

С и

концентрации ингибитора

U(Ui<Un<U

•

точка

колебательно неустойчива. При должном виде функции

/(С,

U) в

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

поддерживающихся колебаний концентраций

С и U.

Режим соот-

С-0

Рис.

7.7.

Нульизоклины системы

(7.23)

на

фазовой плоскости

(U,

С).

а: 1 —

положения изоклины

—

при четырех значениях управляю-

щего параметра S\

> S

; б

— фазовая траектория системы

(7.23)

ходе

релаксации

устойчивой

т. В

после внешнего возму-

щения.

ветствует

автономному ритмическому испусканию цАМФ

дости-

гается при определенном истощении запаса предшественника инги-

битора

S. В

таком режиме работают клетки, автономно испускаю-

щие

цАМФ

образующие центры волнового процесса.

На

падающем участке изоклины С

= 0

система,

как

нетрудно

видеть, способна

возбуждению внешним стимулом цАМФ конеч-

16 Заказ №

241

ной

амплитуды. Действительно, в окрестности точки Q (см.

рис.

7.7 6) сравнительно небольшого начального возмущения АС

достаточно, чтобы выйти за пределы притяжения стационарной

точки В.

В этом

случае

при релаксации к равновесному состоянию си-

стема совершает цикл, в

ходе

которого сначала нарабатывается,

а затем распадается цАМФ. Эта ситуация соответствует реакции

возбуждения амеб в ответ на внешнее возмущение.

Таким

образом, оба наблюдаемых режима в модели естественно

следуют

за монотонным изменением параметра S, уменьшение

которого зависит от истощения запасов пищи. В местах клеточных

сгустков значение S первым становится ниже порога возбуждения

(точка М на рис. 7.7),

тогда

как в остальной части пространства

клетки

находятся в возбудимом состоянии.

Определим условия распространения автоколебаний в прост-

ранстве. Рассмотрим отдельный сгусток (кластер) клеток, пола-

гая,

что приближенно в нем 5 = const. Поскольку перемещение

клеток

мало по сравнению с диффузией цАМФ в среде, движе-

нием

клеток в модели пренебрегаем.

Добавляя к локальной кинетике

(7.24)

член, отвечающий диф-

фузии цАМФ, получим систему

(7.24)

vSC - at/.

Граничным условием для кругового кластера

будем

считать от-

сутствие диффузионного потока цАМФ:

= _D

C = 0.

(7.25)

Стационарное

состояние системы — точка (С

, U

Рассмотрим «мягкий» режим зарождения автоколебаний, огра-

ничиваясь

наинизшими нелинейными вкладами в функцию /(С, U).

В приближении малых амплитуд функция f(C, U) может быть

представлена в виде

/(С,

U) =

/(Со,

С/о)

fcc-\fu\u-

ус\

(7.26)

где с — С — С

; и = U — U

При

заданной f(C, U) система

(7.24)

приобретает вид

•W=(fe ~

Я)

с - | fa \и - ус

(7.27)

vcS — мы,

242

или

где

Ctt

— pc

3-yCfC

—

£>

cov

= 0,

(7.28)

= SpurF,

fc-Ч

Re*-

Прежде

чем

работать

нелинейной моделью (7.27), определим

характер поведения малых отклонений

от

стационарной точки

линейном приближении.

Это

дает

представление

ведущих

сте-

пенях свободы-—«нормальных

модах»

си-

стемы. Пространственная зависимость

нормальных

мод

определяется собствен-

ными

функциями оператора Лапласа

1 д д . 1 д

= з— г

-=—

у- , „ в

области

Г ОГ ОГ Г

O(f

с круговой границей; временная зави-

симость экспоненциальна. Таким обра-

зом,

решение линеаризованной систе-

мы

(7.27)

для

малых отклонений

от

Рис.

7.8.

Графический

вид

дисперсионной зависи-

мости

ReA(£)

для

линеаризованной системы (7.30).

стационарной

точки

(С

, £/о)

может быть задано

виде суперпози-

ции

нормальных

мод

Упт

= е

Чкпт)

(к

пт

г)е

(7.29)

где I

—

функция Бесселя

n-го

порядка;

—

разрешенные

значения

волнового числа

при

заданных граничных условиях

(т —

номер локального экстремума функции

/„).

Дисперсионное уравнение

для X(k)

находится после подстановки

(7.29)

линеаризованную около

(С

, U

)

систему (7.27). Графи-

ческий

вид

зависимости

X{k)

приведен

на рис. 7.8;

кривые

/, //,

///

на

рисунке соответствуют устойчивости однородного состоя-

ния,

критической точке

длинноволновой колебательной неустой-

чивости.

Существенные нормальные моды возмущений однородного

со-

стояния—

это те из

набора (7.29),

для

которых

при

заданной

сте-

пени

закритичности

< £

кр

. При

любом размере кластера

дан-

ному условию удовлетворяет мода, соответствующая однородному

возмущению

k = 0.

Пока

мало,

других

мод с

положительным

инкрементом

нет. С

увеличением

значения

убывают

как 1/R.

Захождение

интервал

(0, &

кр

)

нетривиальных значений волнового

16*

243

числа происходит

следующем порядке:

ku,

&02, &12, соответственно

моды

Vn, Y02, V12

последовательно приобретают положительный

инкремент.

Иной

способ возбуждения

той же

последовательности

состоит

уменьшении бифуркационного параметра

р, в

этом

слу-

чае порог

кр

сдвигается

область больших

По

мере увеличения числа нерелаксирующих

мод

анализ урав-

нения

(7.28)

сильно осложняется [164]. Если оставить лишь

две

моды

Too и Тц,

возбуждающиеся первыми,

то уже в

этом

при-

ближении возможна теоретическая трактовка экспериментально

наблюдаемых режимов.

Полное

представление решения уравнения (7.28), включающее

Рис.

7.9.

Диаграмма ветвления

автоколебательных режимов

при

увеличении параметра

р в

(7.30).

—

устойчивые режимы;

—

неустой-

чивые;

А —

синхронные автоколебания;

—

стоячая волна,

С —

круговая

бегу-

щая

волна.

, имеет

все возможные действительные комбинации

мод

вид

(г,

ф, t) = S(t)slnQt + C(t)cosQt-\-XM(t)I

(t) /,

cos

sin

(t) /,

sin

-f

X^(t)[

r)sinQtcos((

(7.30)

В силу близости

порогу колебательной неустойчивости зависи-

мость амплитуд

от

времени должна быть медленной

по

сравнению

с осциллирующей частью.

После

подстановки решения

(7.29)

(7.28), использования соот-

ношений

ортогональности

2л

jdqnM>

I, fe = 0; 1; 2, гфк,

где

Yo = 1,

~V\

—1\ sin

ф,

= h

cos

ф, отбрасывания высших

гар-

моник,

получаем систему обыкновенных дифференциальных урав-

нений

для

амплитуд

S, С, X^(i= 1; 4).

Укажем существенные стационарные решения [161].

Первыми

при

прохождении бифуркационного параметра через

нуль

из

исходного стационарного состояния возникают простран-

ственно-однородные колебания

амплитудой

С =

(2p/3y)'b(S,

Х<*>

0):

(г,

ф, 0 =

Ccosfi(

= a +

p(0,

(j>0.

(7.31)

Данное решение

существует

при

закритической бифуркации

р > 0,

если

v > 0

(рис.

7.9).

244

Два

других

решения ответвляются от уже неустойчивого одно-

родного состояния при р — pi

kw.

Это круговая бегущая

волна

(7.32)

{г, Ф, 0 = A",/,

(fe,,r)

cos (Qe + ф),

где

стоячая волна вида

(г, ф, О =

cos

Of cos ф,

(7.33)

Рис.

7.10. Спиральный фронт

концентрацион-

ной

волны цАМФ, возбуждаемой круговой вол-

ной

внутри кластера.

Стрелками показаны направления перемещения

амеб

на

переднем фронте волны цАМФ.

где

Остальные стационарные решения получаются из уже перечислен-

ных путем сдвига точки отсчета фазы ф и времени /.

Исследование на асимптотическую устойчивость

дает

следую-

щие

результаты: режим синхронных колебаний

(7.31)

устойчив

области своего существования, стоячая волна

(7.33)—всюду

не-

устойчива, бегущая круговая волна

(7.32)

становится устойчивой

вследствие вторичной бифуркации на определенном удалении от

порога зарождения.

В возбудимой среде, окружающей кластер, круговое пятно,,

колеблющееся как целое, порождает разбегающиеся концентри-

ческими

кругами волны, а бегущая по круговой границе кластера

волна возбуждения — спиральные вращающиеся волны цАМФ.

Градиент концентрации цАМФ в каждой точке агрегационного

поля

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

фронта

концентрационной волны цАМФ перемещения амеб про-

исходят по спиральным траекториям, сходящимся к круговой гра-

нице

кластера (рис.

7.10).

Клетки, первоначально находившиеся

внутри кластера, двигаясь по градиенту цАМФ, перемещаются

периферии круга, так что на месте прежнего клеточного кластера

образуется пустое пространство. Описанные результаты теоретиче-

245.

«кого анализа хорошо соответствуют экспериментально наблюдав-

шимся

картинам агрегационного поля амеб D.

dlscoideum

[161].

Таким

образом, модель реакционно-диффузионного типа каче-

ственно хорошо описывает спонтанное возникновение коллектив-

ных процессов, связанных с химической сигнализацией в популя-

ции

подвижных одноклеточных амеб

Dictyostelium

discoideum.

Возникновение

таких явлений нередко

«портит»

результаты экоси-

стемных экспериментов, так, например, колебания «хищник—

жертва»

в системе

Dictyostelium

discoideum

— Е.

coli,

описанные

п. 4.4, прекратились, по-видимому, именно из-за агрегационного

поведения амеб D.

discoideum.

Агрегационные явления в популяциях исследуются в настоящее

время не только на модельной системе D.

discoideum,

но и в сооб-

ществах подвижных бактерий

Escherichia

coli,

а также на неко-

торых

других

биологических объектах.

7.4.

«Пятнистые структуры»

в биологических

сообществах

Агрегирование особей в пространстве может наблюдаться не

только в отдельные моменты времени при неблагоприятных усло-

виях среды, но существовать достаточно длительное время, обра-

зуя стационарные пространственно неоднородные структуры. Воз-

никновение

в популяциях неоднородных структур, обычно имеющих

форму «пятен», отражает взаимодействие

двух

противоположно

направленных факторов:

необходимости концентрирования особей одного вида в группы

для повышения выживаемости за счет кооперативных эффектов

(совместное выделение метаболитов для кондиционирования среды

химической защиты; увеличение надежности сигнального опове-

щения

об опасности, улучшение условий размножения и т. д.);

внутри- и межвидовой борьбы особей за индивидуальное выжи-

вание в условиях ограниченных ресурсов пищи и пространства.

Ниже

будет

рассмотрено явление образования пространствен-

ной

структуры в водной экосистеме, отражающее особенности про-

странственной группировки особей в популяциях.

«Пятна»

фитопланктона.

Внимание гидробиологов давно при-

влекало необъяснимое явление, наблюдающееся во

всех

доста-

точно крупных водоемах — несмотря на подвижность и хорошую

перемешиваемость поверхностных слоев воды,

«цветение»

микро-

водорослей планктона происходит не равномерно по всей площади

водоема, а отдельными «пятнами» размером от нескольких метров

до нескольких километров, разделенными пространствами, обед-

ненными

планктоном. Наблюдения показали, что «пятнистость»

развития водорослей обусловлена не только неоднородностью рас-

пределения биогенных элементов и температуры и

существует

также в местах с исходно равномерным распределением этих по-

казателей. Предпринимались многочисленные попытки объяснить

246

это

явление выеданием водорослей зоопланктоном, которому при-

писывалась тенденция мигрировать в направлении градиента плот-

ности

пищи. Математические модели типа «хищник—жертва»

с диффузией и миграцией хищников указывают на возможность

существования неоднородных

структур

[11, 129, 130], тем не ме-

нее с биологической точки зрения такое объяснение не может счи-

таться удовлетворительным. Действительно, экспериментально

установлено, что зоопланктон активно избегает больших скоплений

фитопланктона,

поскольку в области пятна цветения создаются

высокие концентрации экзометаболитов микроводорослей (фе-

нолы,

антибиотики и

другие

аллелопатические вещества), что

ряде случаев приводит к угнетению и гибели не только зоопланк-

тона, но и рыб [134]. Выеданию зоопланктоном обычно подверга-

ются разреженные или отмирающие популяции микроводорослей.

Таким

образом, «пятнистость» фитопланктона представляет

собой,

по-видимому, чисто популяционный эффект. Словесная мо-

дель появления пятна цветения может быть сформулирована сле-

дующим образом. Как было показано в п. 2.8, для перехода по-

пуляции

микроводорослей в стадию активного роста необходимо

наличие определенной минимальной концентрации собственных

экзометаболитов в среде. Экзометаболиты выделяются в водную

среду

самими клетками, при этом время достижения необходимой

концентрации

определяется, с одной стороны, начальной локальной

концентрацией

клеток, а с

другой

стороны — скоростью рассеяния

экзометаболитов в окружающей среде. При возникновении локаль-

ных благоприятных условий (слабая перемешиваемость вод, мик-

роскопление

клеток, повышенная температура) успевает нараба-

тываться необходимая концентрация экзометаболитов, и начина-

ется экспоненциальное размножение клеток, приводящее к обра-

зованию пятна. Косвенным подтверждением описанного механизма

является тот факт, что интенсивное цветение в водоеме начинается

после нескольких дней теплой и безветренной погоды; с

другой

стороны,

для предотвращения цветения вод в водохранилищах

используется метод искусственного интенсивного перемешива-

ния

вод.

Дальнейшее увеличение пятна цветения определяется запасом

биогенных элементов в воде и скоростью их выедания водорос-

лями.

Первоначально в центре пятна микроводоросли размножа-

ются наиболее быстро, по краям — медленнее, затем, после исто-

щения

биогенов в центре пятна, водоросли начинают отмирать,

пятно

принимает вид расходящегося кольца, которое постепенно

становится неустойчивым и разрушается.

Полная

модель образования и разрушения пятна цветения

пока

не создана, наиболее близкой по механизму описания про-

цесса является модель Ю. А. Домбровского и Г. С. Маркмана

[63,64].

Двумерная модель распределения фитопланктона в водоеме

учетом

диффузии, взаимодействия с метаболитами и потребле-

ния

лимитирующего биогенного элемента питания может быть опи-

24 Г

сана

следующей

системой уравнений:

(7.34)

-=0

где X(r, t)—плотность фитопланктона; S(r, t)—концентрация

лимитирующего биогенного элемента; М(г, /)—концентрация

метаболитов; D

и D

— коэффициенты диффузии для фитопланк-

тона и биогенов; £>

— коэффициент диффузии метаболитов; С —

скорость поступления биогенов; |3 — коэффициент экскреции ме-

таболитов; б — коэффициент их распада. Зависимость скорости

роста а от концентрации метаболитов М предполагается линейной

функцией

a = a

+ ?М-

Если

предположить, что метаболиты не влияют на рост водо-

рослей (Y =

0),

TO первые два уравнения системы

(7.34)

отделя-

ются и поведение решений оказывается простым: величины X и S

независимо от начальных условий и параметров системы стремятся

распределению, однородному по пространству и стационарному

но

времени. Аналогичный вид имеют решения модели при у < О

(ингибирующее действие метаболитов).

Наиболее интересные динамические режимы появляются в мо-

дели при стимулирующем действии метаболитов: у > 0.

При

математическом исследовании модели целесообразно пе-

рейти к безразмерным переменным

у/

(далее

штрихи у безразмерных переменных

будут

опускаться).

В безразмерных переменных система

(7.34)

приобретает вид

(7.35)

248

где

а = а„С/б

, р = ус

/Ь\

Уравнения

(7.35)

имеют единственное однородное по пространству

решение .

=\,

М, = 4", S

= (a +

-^-)"'.

(7.36)

Изучение модели начнем с рассмотрения точечного варианта

системы (7.35), отражающего поведение траекторий во времени.

Соответствующая

(7.35)

система обыкновенных дифференциаль-

ных уравнений имеет вид

X = X (а + рМ) S - X,

S = —X(a +

pM)S+\,

(7.37).

М = ЬХ — М.

Устойчивость пространственно-однородного решения

(7.36)

опре-

деляется в линейном приближении характеристическим уравнением

А»

+ (ц +

А.»

[ц(1

+ /) - / +

-2L]

А.

/ц

О,

(7.38)

При

\х, близких к нулю, все Re X < 0, особая точка устойчива.

увеличением р устойчивость нарушается, критические значения

(х можно определить из равенства нулю условия Раусса—Гурвица

/)ц

+ /(/-1)ц

+ /(а-/)ц + а/

= 0. (7.39).

Это уравнение в зависимости от значений параметров (а, /) либо

не

имеет ни одного положительного

корня,

либо имеет два положи-

тельных

корня:

0 < ц** < ц*.

Из

(7.38),

(7.39)

следует,

что при ц < щ* особая точка

(7.36)

есть устойчивый фокус, при ц, > ц** — устойчивый узел. При

(j,

= ц* или ц = щ*,

(7.36)

— есть особая точка типа «центр».

При

(А** < (J. < (л* особая точка (7.36)—седло-фокус второго рода

(входящий одномерный

«ус»

и двумерный выходящий); как пока-

зали численные расчеты, особая точка в этом

случае

окружена

предельным циклом, на который наматываются все траектории.

Форма

предельного цикла в проекции на плоскость (X, S) и дина-

мика

решений точечной модели

(7.37)

показаны на рис. 7.11.

Рассмотрим, далее, пространственно-временную динамику мо-

дели (7.35). Устойчивость однородного по пространству решения

(7.36)

будем

определять при различных значениях коэффициента

экскреции

Ь, который

будет

бифуркационным параметром.

249

Вводя возмущения

g(p, t),

(р, /) (X = Х

+ |, S =

+ £) в

(7.35),

линейном приближении получим

дг\

^ ... . ,

<Э

т1

(7.40)

X,S,M

Рис.

7.11. Траектории модели (7.37) на фазовой плоскости (а) и во

времени

(б).

X — биомасса; М — концентрация лимитирующего биогена; S — концентрация

метаболита; а — 0,1; р = 3; Ь = 0,1; 1 = 1.

Считая для простоты, возмущения периодическими по простран-

ственной координате р, решение (7.40) ищем в виде

(Р.

0. Л(Р> 0.

(P.

(I, rj, г)созшр,

(7.41)

где |, г], ^ — неизвестные числовые коэффициенты, волновое число

w изменяется непрерывно от нуля до бесконечности. Протяжен-

ность водоема предполагается достаточно большой.

Подставляя

(7.41)

в (7.40) и приравнивая нулю определитель

получившейся однородной системы, получим характеристическое

уравнение для определения X:

+ к' [\i + I + и (1 + &

+ d

)] + I [и

(d, + d

) + d,d

] +

и [(d,n + l(d

+ d

) + a) + ц + la + la/a - I] +

(a + ud,a + u44*) (t + u) - I (1 -

fl/,i)

= 0, (7.42)

где и = до

При

6 = 0 (|i, = a), все ReX < 0 и однородное решение (7.36)

устойчиво. С ростом Ь свободный член в (7.42) меняет знак, про-

ходя при некотором Ъ* через нуль, что соответствует переходу од-

:250