Математическое моделирование в экологии:
Историко-методологический анализ

6.4. Нелинейная демографическая динамика

В заключение дадим пример работы, которая с методической стороны представляется нам достаточно совершенной. В ней есть всё — обширный экологический эксперимент, математическая модель с реалистической оценкой ее возможностей, статистическая обработка результатов с доказательствами статистической однородности; всё, разумеется, кроме народнохозяйственной пользы. Речь идет о работе [85], носящей несколько пышное название: «Нелинейная демографическая динамика: математические модели, статистические методы и биологические эксперименты». Всё перечисленное в заголовке в работе действительно есть, хотя и не в применении к широкому кругу экологических проблем, а лишь к динамике популяции мучного жучка Tribolium castaneum. Жучок Tribolium является классическим экспериментальным объектом, которым биологи занимаются более 60 лет. В списке литературы к статье приводится двухтомная монография по биологии Tribolium (издана в 1974 году в Оксфорде). Судя по этому же списку, двое из авторов статьи (R. F. Costantino и R. A. Desharnais) — биологи, работавшие с Tribolium не менее 20 лет, а двое других авторов — специалисты по статистической обработке. Что же может сделать хороший коллектив биологов и математиков с хорошими же экспериментальными данными о классическом лабораторном объекте?

Рассматриваемая работа посвящена исключительно математической обработке. Сами экспериментальные данные получены еще в 1980 году (см. [86]). Были проведены 13 экспериментов с культурами Tribolium castaneum Herbst с одинаковым во всех 13 случаях начальным количеством личинок, куколок и взрослых особей. Каждая культура содержалась в полупинтовой (237 мл) молочной бутылке с питательной смесью. Через каждые 2 недели содержимое бутылки пересчитывалось (кроме яиц, которые пересчитать невозможно) и целиком, включая яйца, помещалось в новую питательную среду. Это делалось 19 раз (т. е. в течение 38 недель). С четырьмя из 13 культур не делалось больше ничего (группа контроля), а в 9 оставшихся (разбитых на 3 группы по 3 повторности) вносились в конце 10-ой недели демографические изменения: в первой групп добавлено по 100 взрослых особей; во второй группе выброшены все взрослые особи; в третьей группе оставлены лишь взрослые особи.

Работа, которую мы излагаем, представляет собой исключение из сложившейся вредной традиции, о которой мы говорили еще в связи с работами Гаузе: результаты эксперимента (т. е. таблицы учетов численностей) опубликованы, хоть и не сразу (см. [87]).

В жизненном цикле жука Tribolium различают стадии яйца (продолжительность 2-4 дня), кормящейся личинки (14 дней), некормящейся личинки, куколки и молодого взрослого (вместе три стадии примерно 14 дней), и наконец, взрослого жука. Яйца в работе не учитывались, кормящиеся личинки называются просто личинками (larvae), число их в момент t обозначается L_t; третья, четвертая пятая стадии объединены под названием куколки (pupae), число их P_t; число взрослых (adult) особей обозначается A_t (t = 1, 2, ..., 19 — порядковый номер просчета культуры).

Важнейшим фактором, определяющим динамику популяции, является каннибализм: личинки едят яйца и куколок (но не личинок и не взрослых), взрослые едят яйца, личинок и куколок. Каннибализм сходен с хищничеством в смысле влияния на динамику популяции, но в то время как в лабораторных сообществах с хищниками и жертвами трудно добиться незатухающих циклических колебаний численностей, в культурах Tribolium подобные колебания возникают легко. Количественно каннибализм предлагается описывать экспоненциальной моделью: если c_el — вероятность того, что за единицу времени (т. е. 2 недели) данное яйцо съест данная личинка, то при наличии L_t личинок в момент t через единицу времени останется доля exp(–c_elL_t) несъеденных яиц. Из этого (довольно грубого из-за перемены L_t в течение одной единицы времени) рассуждения получаются уравнения модели. На этот раз это не дифференциальные уравнения, а рекуррентные (с дискретным временем).

Например, пусть в момент t имеется L_t личинок, P_t куколок и A_t взрослых. Единица времени выбрана так, что в момент t + 1 все L_t личинок превратятся в куколок (кроме той части, которая погибнет), все же куколки (кроме погибших) превратятся во взрослых. Каким же будет число личинок L_t+1? Они могут возникнуть лишь из яиц, отложенных A_t взрослыми (практически все эти яйца — кроме погибших — за 2 недели превратятся в личинок). Если бы яйца не ели, то получилось бы bA_t личинок (где b — коэффициент плодовитости с учетом гибели яиц). С учетом же выедания яиц получаем

(Выедание самих личинок, а не яиц, взрослыми авторы считают возможным не учитывать.)

Для числа куколок предлагается уравнение

где μ_l — коэффициент, учитывающий смертность личинок (их выедание в модели не учитывается).

Наконец, для числа взрослых особей предлагается уравнение

(учитывается выедание куколок взрослыми с коэффициентом c_pa и смертность взрослых с коэффициентом μ_a).

Модель (6.4.1) — (6.4.3) получается нелинейной за счет присутствия экспонент. Она имеет 6 параметров, которые нужно определять по данным наблюдений: b, μ_a, μ_l, c_ea, c_el, c_pa. Авторы работы сначала довольно долго излагают все сложности динамики этой детерминированной математической модели, которые могут возникнуть в зависимости от значений параметров. Очевидно, для того, чтобы взяться за новую обработку наблюдательных данных, пролежавших с 1980 года, нужно было сначала осознать красоту детерминированной нелинейной динамики. Но детерминированная модель, как вскоре обнаруживается, вообще не может описать наблюдательные данные: в модель приходится ввести случайные возмущения (как увидим ниже, весьма значительные). Детерминированная динамика, очевидно, необходима лишь для вдохновения.

Случайные возмущения, о которых идет речь, представляют собой просто множители exp(E_1t), exp(E_2t), exp(E_3t), на которые нужно домножить правые части, соответственно, уравнений (6.4.1) — (6.4.3), чтобы получить реально рассматриваемую модель со случайными воздействиями. Считается, что вектор E_t = (E_1t, E_2t, E_3t) при каждом t имеет многомерное нормальное распределение со средним 0 и некоторой матрицей ковариаций Σ = (σ_ij), причем при различных значениях t векторы E_t независимы. Таким образом, добавляются еще 6 параметров — элементы (симметричной) матрицы ковариаций.

Во введении авторы пишут (несколько пышным стилем, как и в заглавии работы), что они собираются «строго связать нелинейную демографическую модель с биологическими данными с помощью заново развитых статистических методов для нелинейных случайных процессов». Но вообще-то нелинейных случайных процессов в используемой модели не оказывается, если перейти к логарифмам, что авторы вскоре и делают. Некоторым усложнением является предполагаемая зависимость (при каждом t) случайных вмешательств E_1t, E_2t, E_3t друг от друга. Однако даваемые в работе оценки ковариаций σ_ij приводят к весьма скромным величинам для соответствующих коэффициентов корреляции: r₁₂ = 0,081, r₁₃ = 0,176, r₂₃ = 0,118. Практики знают, что коэффициент корреляции, меньший (по абсолютной величине), чем 0,3, плохо оценивается по наблюдательным данным и не имеет практических последствий. Иными словами, мало бы что изменилось, если бы случайные вмешательства считать независимыми друг от друга.

Если же считать E_1t, E_2t, E_3t независимыми, то задача оценки параметров модели делается очень простой (обычный регрессионный анализ). Действительно, например первое уравнение модели (домноженное на случайный фактор) после логарифмирования приобретает вид

Таким образом, беря из наблюдений значения A_t, L_t, L_t+1, получаем, что задача оценки параметров lnb, c_ea, c_el является стандартной задачей линейного регрессионного анализа.

Авторы статьи и сами в качестве первого приближения пользуются регрессионным анализом, но основным для оценки параметров считают более сложный метод максимума правдоподобия. В данном случае его, видимо, нельзя реализовать просто путем обращения к какому-нибудь стандартному статистическому пакету программ, и действительно потребовались новые программные разработки. Однако сопоставим окончательно принятые авторами значения параметров с теми значениями, которые получились бы, если бы просто усреднить 4 варианта стандартного регрессионного анализа, которые также приводятся авторами (см. [85], стр. 266, таблица 1).

 

Достаточно ясно, что гора в виде «заново разработанных методов» в данном случае родила мышь.

Теперь самое главное — каковы же реальные достижения?

Выпишем, скажем, первое уравнение системы с учетом случайностей, а именно

Дисперсия случайной добавки E_1t оценена в работе как 0,2771, следовательно, стандартное уклонение составляет 0,526. Допустим, что параметры модели, входящие в уравнение (6.4.4), оценены без ошибки и поставим задачу прогноза числа куколок L_t+1 по известным значениям величин A_t и L_t, т. е. на один шаг вперед. Единственное, что можно сделать для такого прогноза — это положить равной нулю принципиально непредсказуемую величину E_1t. При этом порядок ошибки в прогнозе логарифма lnL_t+1 определяется стандартным уклонением 0,526, что соответствует ошибке прогноза в самой величине L_t+1, равной exp(0,526), т. е. в 1,69 раза. Итак, ошибка прогноза на один шаг вперед может составлять несколько десятков процентов.

Для второго уравнения вместо дисперсии ошибки 0,2771 нужно взять 0,4284 (т. е. вместо 0,526 взять 0,654, а вместо 1,69 — 1,92). Иными словами, возможные ошибки прогноза еще несколько больше.

Для третьего уравнения исходная цифра иная: 0,01112, а возможная ошибка прогноза составляет примерно 10%.

Для оценки качества прогноза возможные ошибки важны не сами по себе, а в сопоставлении с тем, насколько может (за один временной шаг) измениться сама прогнозируемая величина. Рассматривая графики изменения величин L_t, P_t, A_t, можно увидеть, что величины L_t и P_t за один шаг по времени могут измениться в несколько раз, так что прогноз для них с ошибкой в два раза все-таки есть прогноз. Величина же A_t за один шаг по времени меняется на десятки процентов, так что прогноз для нее с ошибкой порядка 10% есть прогноз примерно того же качества, конечно, не особенно хорошего.

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

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

Правда, ориентация авторов на детерминированную нелинейную динамику приводит к излишне тщательному изучению вопроса, в какую область попадают параметры системы. Дискретная детерминированная система вообще имеет, как теперь говорят, «неоднозначное» отношение к реальной динамике, прежде всего потому, что дискретное приближение является очень грубым, потому что выбранный шаг по времени довольно велик. Например, пусть численность личинок L_t = 50, численность L_t+1 = 200 (такой скачок численностей за один временной шаг вполне возможен). Если интенсивность выедания яиц личинками равна 0,01, то какую же долю яиц оставили несъеденными эти личинки : то ли exp(–0,5), то ли exp(–2)? Вероятно, желательно было бы иметь модель, более близкую к непрерывной (с меньшим шагом по времени).

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