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

7.2. Численный эксперимент

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

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

Имеются два конкурирующих вида, численности которых обозначаются N₁(t) и N₂(t). Каким-то образом (это часть протокола эксперимента, которая будет обсуждаться в дальнейшем) задаются начальные численности N₁(0) и N₂(0) и шаг по времени h. Приращения

моделируются по формулам

где ξ_1i и ξ_2i — независимые между собой при всех i стандартные нормальные случайные величины, а B_i и D_i определяются таким способом, чтобы разность B_i – D_i представляла собой дискретный аналог уравнения конкуренции двух видов. Конкретно

В формулах (7.2) коэффициенты b_i — логарифмические скорости роста, K₁ и K₂ — предельные емкости среды для каждого из видов, α и β — константы конкуренции. При этом величины b_i, α и β — константы, специфичные для рассматриваемых видов и данной среды обитания (хотелось бы, чтобы эти величины были постоянными, насколько это возможно в биологии), а K₁ и K₂ жестко связаны с объемом среды, в котором производится опыт.

Уравнения (7.1) и (7.2) — математическая модель динамики численностей (модель Лотки-Вольтерра плюс случайные отклонения по Реншоу), но в нее не включена модель ошибок учета численностей видов. Учет производится путем отбора проб в какие-то моменты времени, которые удобно выбрать целыми кратными шагу по времени h. Объем пробы мы характеризуем примерным числом просчитываемых особей каждого вида m₁ и m₂; для простоты мы примем, что m₁ = m₂ = m и будем ориентироваться на среднеквадратическую относительную ошибку подсчета численностей ε = m^–1/2. Если число m велико, причем просчитанные особи уничтожаются, то нужно в уравнения модели внести соответствующую поправку (как это делал Гаузе), но мы пренебрегаем этим как для простоты, так и из-за того, что (как мы увидим ниже) можно ориентироваться на сравнительно небольшие количества просчитываемых особей. Таким образом, если моделированные численности равны N_i, то просчитанные в результате эксперимента численности принимаются равными

где ξ_3i — стандартные нормальные величины, не зависящие от ξ_1i и ξ_2i, которые участвуют в моделировании по формуле (7.1).

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

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

Перечислим сначала, какие варианты протокола априори представляются нежелательными. Во-первых, создание одновидовых популяций (как это делал Гаузе) нежелательно по той причине, что в этом случае мы экспериментируем на самой границе фазовой плоскости, принимая N₁ = 0 или N₂ = 0. Но на границе фазовой плоскости и где-то в ее средней части вполне могут действовать различные математические модели для динамики численностей. Лучше, если создаваемые в эксперименте численности видов будут величинами примерно одного порядка. Во-вторых, все, что нам в настоящее время известно об участке стабилизации численностей, указывает на нежелательность экспериментирования в этих условиях (с целью оценки параметров конкуренции). Гаузе наблюдал на этом участке значительные колебания численностей, которые не укладываются ни в его детерминированную модель, ни даже в модель со случайностями по Реншоу. Здесь происходит, по-видимому, нечто, не описываемое просто конкуренцией. Таким образом, тот протокол эксперимента, который применял Гаузе, ставится под сомнение уже по априорным соображениям. Итак, желательно экспериментировать на участке роста численностей с примерно равными по порядку величины численностями обоих видов.

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

Для численного моделирования по формулам (7.1) — (7.3) были приняты параметры модели, привязанные к данным Гаузе по инфузориям Paramecium aurelia и Paramecium bursaria, однако существенно поправленные в отношении скоростей роста. Приняты следующие значения параметров:

Шаг h по времени лучше сделать достаточно малым, и мы приняли его равным 0,1 суток. Периодичность учетов сохранили, как у Гаузе, 1 раз в сутки, т. е. через 10h. Емкости среды мы варьировали тремя способами:

I группа (большой микрокосм) — K₁ = 4100, K₂ = 3600 (это примерно соответствует экспериментам самого Гаузе);

II группа (средний микрокосм) — ровно в 10 раз меньше, т. е. K₁ = 410, K₂ = 360;

III группа (малый микрокосм) — еще в 10 раз меньше, т. е. K₁ = 41, K₂ = 36.

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

где k (номер суток) принимает значения 1, 2, 3, а i (номер вида) принимает значения 1 и 2. (Это означает, что в каждом отдельном варианте опыта мы получаем только шесть чисел — значений скоростей, но еще остается почти бесконечная возможность разнообразить различные варианты опыта за счет выбора различных начальных значений, числа повторностей и точности учета). Кроме скоростей, мы будем рассматривать также оценки средних за сутки численностей M_i(k) по формуле

Качественно явление конкуренции состоит в том, что большим значениям скорости роста одного вида отвечают меньшие значения средних численностей другого вида и наоборот. Напрашивается мысль о корреляции величин V₁(k) с M₂(k) и величин V₂(k) с M₁(k) при различных значениях k. Некоторые скажут, что вряд ли имеет смысл коррелировать между собой три числа, но не следует забывать, что в нашем распоряжении находится число повторностей опыта с фиксированными начальными условиями (для биолога каждая повторность означает новый счет особей). Взяв не слишком много повторностей, мы без труда получили устойчивые отрицательные значения коэффициентов корреляций порядка –0,9 или еще левее. Но вскоре мы сообразили, что эти корреляции бессмысленны. В течение опыта численности растут, а скорости роста падают, так что отрицательные корреляции будут наблюдаться и при отсутствии конкуренции между видами. Это было подтверждено, когда в численных экспериментах мы положили параметры конкуренции α и β равными нулю. Отрицательные корреляции остались (типичный случай «ложной корреляции», когда наличие отрицательной корреляции могло бы интерпретироваться как доказательство наличия конкуренции, а на самом деле объясняется иной причиной). В результате протокол эксперимента с корреляциями численностей и скоростей роста по ансамблю моментов времени был оставлен.

Вместо этих корреляций мы попытались рассмотреть некоторые другие. Создадим несколько повторностей эксперимента с одними и теми же начальными условиями. Наблюдаемые численности будут различными в различных повторностях (в математической модели — за счет выбора разных случайных чисел, используемых при имитации). Подобные повторности с одинаковыми начальными условиями создавал в своих опытах Гаузе — и спасибо ему за это, потому что только эти данные оказалось возможным использовать для тестирования модели Реншоу. Но хорош ли такой протокол эксперимента для определения параметров модели конкуренции? В частности, сможем ли мы хотя бы доказать наличие конкуренции, используя данные по подобным повторностям? Предлагается скоррелировать при каждом фиксированном k (номер суток) скорость V₁(k) с численностью другого вида M₂(k) и аналогично V₂(k) c M₁(k), ожидая отрицательных корреляций. Но оказалось, что даже при ошибке учета численности ε = 0 примерно половина получаемых коэффициентов корреляций положительна, а половина отрицательна. Это объясняется тем, что при фиксированных начальных условиях создается (в разных реализациях) слишком малый разброс численностей: например, численности порядка 1000 в разных повторностях колеблются лишь на несколько десятков. Поскольку конкуренция сказывается лишь при заметных относительных колебаниях численностей, то мы и не в состоянии ее заметить при таком протоколе эксперимента. Делается вывод, что экспериментировать надо активнее: менять начальные численности (не создавая, однако, одновидовых популяций, а соблюдая равенство численностей в смысле порядка величины).

Например, предлагается для каждой группы (т. е. определенного размера микрокосма) испытать по 4 эксперимента с разными начальными численностями (всего по одной повторности для каждого начального условия). Именно, для первой группы (большой микрокосм) положим N₁(0) = 250, а для N₂(0) возьмем 4 различных значения 100, 200, 300 и 400. Для второй группы сократим эти числа в 10 раз, а для третьей группы возьмем N₁(0) = 3, а N₂(0) сократим еще в 10 раз (т. е. возьмем значения 1, 2, 3, 4). Осуществим теперь с помощью стандартных программ линейную регрессию скоростей роста за данные сутки на значения средних численностей обоих видов за те же самые сутки. Более подробно это означает, что мы выписываем уравнения

Коэффициенты A_i, B_i, C_i в этих уравнениях определяются методом наименьших квадратов. (Для определения любой тройки, например, A₁, B₁, C₁, имеется 12 наблюдений, отвечающих трем различным значениям k и четырем вариантам начальных условий.) Такой вариант протокола эксперимента оказался на удивление удачным, но для того, чтобы это увидеть, лучше представить данные численных экспериментов в несколько пересчитанном виде — а именно, пересчитать полученные значения коэффициентов в уравнениях (7.6) непосредственно в параметры уравнений конкуренции.

Мы приняли самый простой приближенный способ пересчета. При достаточно малом h (шаг по времени) логарифмические скорости роста в модели конкуренции приближенно представимы в виде

Допустим, что эти выражения достаточно точны и при h = 1, если в их правых частях вместо N₁ и N₂ подставить полусуммы M₁ и M₂. Тогда получится система линейных уравнений, связывающих логарифмические скорости роста со средними численностями, аналогичная системе уравнений линейной регрессии (7.6). Поэтому связь между параметрами регрессии и параметрами уравнений конкуренции должна быть следующей:

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

Таблица 1

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

 

Из таблицы видно, что предложенный довольно грубый способ оценки параметров модели конкуренции дает результаты с точностью, вполне приемлемой для биологических параметров. Ни малый объем микрокосма (группа III), ни большие ошибки при подсчете проб (вплоть до ε = 0,35, что отвечает примерно десятку особей, фактически просчитываемых в пробе) не делают бессмысленными оценки параметров модели. Общее количество подсчетов численностей при данном протоколе эксперимента составляет 4 × 4 × 2 = 32 (4 повторности эксперимента; 4 подсчета для каждой повторности, т. е. начальные численности плюс три подсчета по суткам; 2 вида). Если в каждом подсчете участвуют даже не 10, а 100 особей (это соответствует примерно 10-процентной точности подсчетов), то общее количество просчитанных особей 3200 еще вполне приемлемо. Но стремиться к такой точности отдельных подсчетов нет необходимости, так как при рассматриваемом варианте протокола эксперимента (с применением стандартного регрессионного анализа) независимые случайные ошибки отдельных подсчетов компенсируются статистической обработкой. А для успеха регрессионного анализа следует создать значительно отличающиеся (в разных вариантах эксперимента) начальные численности видов.