EcoSimulation — 04. Учет демографической структуры популяций

 

Д.А. Шабанов
Конспект курса
"
Имитационное моделирование сложных биосистем
(с использованием LibreOffice Calc)
"

Усложнение модели: имитация логистического роста Учет демографической структуры популяций Имитация ограничения доступного количества ресурсов
Имитационное моделирование биосистем-03 Имитационное моделирование биосистем-04 Имитационное моделирование биосистем-05

 

4. Учет демографической структуры популяций

4.1. Новые возможности

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

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

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

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

Чтобы присвоить ячейке (или диапазону ячеек) уникальное имя, надо выделить эту ячейку и пройти по последовательности команд "Вставка/Іменовані вирази/Задати...". Ячейка получит уникальное имя. Некоторые сочетания символов использовать в качестве имени нельзя (например, имя не может совпадать с обозначением какой-либо иной ячейки, т.е. не может состоять из сочетания буквы и номера). В именах символов можно использовать подчеркивания (_), а также, в середине имени, точки и пробелы. Чтобы увидеть список присвоенных в данном листе Calc имен и иметь возможность его отредактировать, надо пройти по пути "Вставка/Іменовані вирази/Керування".

Во-вторых, в качестве средств управления используем полосы прокрутки. Это одна из разновидностей элементов управления, предусмотренных электронными таблицами. Самый простой способ вставить в файл полосу прокрутки — скопировать ее из того файла, где она уже есть. Более сложный (но более правильный) способ таков. Пройдите по пути "Перегляд/Панелі інструментів" и поставьте "галочку" напротив панели "Засоби". На листе Calc появится эта панель. Третий слева элемент в ней — "Засоби керування". Разверните его. Перейдите в режим разработки с помощью кнопки, на которой написано "ОК": "Перемкнути режим розробки". Станут доступны элементы управления, среди которых можно выбрать полосу прокрутки "Смуга прокрутки" (рис. 4.1.). Выбрав ее, очертите какое-то пространство на листе; Calc вставит туда этот элемент управления. Вначале полоса прокрутки является вертикальной. Ее надо настроить. Оставаясь в режиме разработки, выделите полосу прокрутки. Станет активной кнопка "Властивості елемента керування", на которой изображено зубчатое колесо. Нажав ее, вы попадете в диалог "Властивості: Смуга прокрутки". На первой вкладке, "Загальні", можно настроить свойства полосы, в число которых входят размеры, положение, минимальное и максимальное значения, задаваемые с помощью этой полосы, шаг изменения, а также вертикальная или горизонтальная ориентация. Как видно по рис. 4.1, в модели, которую мы строим, используются горизонтальные полосы. На второй вкладке, "Дані", следует указать ту ячейку, значение в которой будет изменяться с помощью этой полосы прокрутки. До того, как вы выйдете из режима разработки, вы не увидите результатов сделанных вами изменений.

Есть одна неприятная особенность полос прокрутки, которая, возможно, проявляется не во всех случаях. Управляя полосой прокрутки, можно задать максимальное значение величины, значение которой изменяет эта полоса "Макс. значення прокрутки". На компьютере автора этого текста реальное максимальное значение, которое можно выставить с помощью полосы прокрутки, ниже того, которое задается в диалоге. Если в диалоге выставить 1000, выйти из режима разработки и проверить, что получилось, можно увидеть, что крайнее левое положение "бегунка" соответствует значению 995 в целевой ячейке. Вероятно, это пример т.н. программного "бага". Чтобы реальное максимальное значение составляло 1000, в окошке "Властивості: Смуга прокрутки/Загальні/Макс. значення прокрутки", приходится выставлять значение 1005.

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

Рис. 4.1. Курсор показывает на кнопку вставки полосы прокрутки. Выделена ячейка B9, для которой задано имя, соответствуюшее обозначению переменной; это имя можно увидеть в поле слева вверку (первая в том ряду, в котором находится строка формул)

Начиная с этой модели, мы будем придерживаться стиля обозначений, который можно назвать "рекомендованным" (в данном курсе), по крайней мере, для популяционно-экологических моделей. Список обозначений, которые используются в таких моделях, приведен здесь. Рассматриваемая совокупность особей в таких моделях состоит из отдельных групп. Для обозначения совокупности в целом используются прописные (большие) буквы, для обозначения больших групп, которые делятся на части, — тоже большие буквы, но с указанием в скобках, к кому они относятся, а для отдельных групп — строчные (маленькие). Для обозначения численности групп на разных этапах цикла могут использоваться разные обозначения. Например, в начале цикла в следующей модели численность группы обозначается буквой αn, а после сокращения — ωn ("движение" от начального состояния до конечного символитзируют первая и последняя буквы греческого алфавита). Начальная численность всей рассматриваемой нами популяции обозначается αN, начальная численность самок — αN(F), численность самок первого возраста — αn(f1). Чтобы обозначить момент времени, используются подстрочные индексы; например, окончательная численность самок второго возраста на 14-м цикле работы модели обозначается ωn(f2)14. В таком случае численность самок после сокращения может обозначаться как ωN(F). Численность пар (pair) с самками второго возраста обозначим как ωN(PF2). Можно было бы использовать и строчную букву p, но тогда ее можно будет спутать с общепринятым обозначением вероятности. Поскольку численность пар вычисляется с учетом количества самцов, будем считать, что мы "имеем право" использовать прописную букву. 

Другие характеристики групп также обозначаются подобным образом. Например, в следующей модели к числу начальных параметров относится плодовитость самок. Для ее обозначения использована буква b (от breed); плодовитость самок второго возраста — b(f2).

4.2. Описание простейшей модели

Итак, рассмотрим структуру простейшей модели, имитирующей динамику популяции из трех поколений (рис. 4.2).

Рис. 4.2. Модель с тремя поколениями и сокращением численности особей до K на каждом цикле

Рассмотрим эту модель подробнее (кроме прочего, ее можно скачать, хотя, конечно, полезнее сделать ее самостоятельно).

В информационном блоке даны короткие пояснения относительно работы модели.

В блоке для входных параметров — 6 ячеек, пояснения к которым видны на рисунке, а присвоенные им имена (сверху вниз) таковы: n_f1_0, n_m1_0, K, b_f2, b_f3 и F_M

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

В столбце A представлены шаги модели (циклы, соответствующие в действительности, например, годам). Первая строка расчетов — нулевая, A17{0}. Вторая ячейка этого столбца и все последующие ("растянутые" из ячейки A18) содержат простейший счетчик: A18{=A17+1}.

В столбце B указывается численность самок первого года в начале цикла — αn(f1). Формулы в ячейках таковы: B17{=n_f1_0}, B18{=ROUNDDOWN(S17/2+RAND();)}; все последующие ячейки этого столбца получены из B18 "растягиванием". Смысл этой формулы прост. В ячейке S17 указано количество потомства, появившегося на предыдущем шаге модели. Половина этого количества — самки. Округление (важное, в том случае, если появляется нечетное количество потомков) носит вероятностный характер (как объяснялось ранее).

В столбце C указывается численность самок второго года в начале цикла — αn(f2). Формулы в ячейках таковы: С17{0} (по условиям все особи в начале работы модели относятся к первому возрасту), С18{=I17}; все последующие ячейки столбца получены из С18 "растягиванием". В ячейке I17 указывается численность самок первого возраста, которые жили в системе год назад (после сокращения численности). Теперь они достигли второго возраста.

Аналогичные формулы введены в столбец D, где указывается численность самок третьего года — αn(f3). Формулы в ячейках таковы: D17{0}, D18{0} (третьегодки не успевают появиться и через год), D19{=J18}; все последующие ячейки столбца получены из D19 "растягиванием". Объяснение этой формы аналогично предыдущему случаю — это численность самок предыдущего возраста годом раньше.

Формулы в столбце E, где указывается численность самцов первого года, αn(m1), определяют число самцов как разность между общим количеством потомков и числом самок в их составе; E18{=S17-B18}. Столбец F аналогичен столбцу C, а G — аналогичен D.

Движемся далее. В столбце H (αN) указывается общая численность всех особей; H17{=SUM(B17:G17)}. Функцию SUM мы еще не рассматривали, но ее применение интуитивно понятно и не требует многословных пояснений.

В столбцах I - N численность всех возрастных классов пересчитывается после сравнения αN и K. Их обозначения аналогичны таковым для столбцов B - G, за тем исключением, что мы рассматриваем уже численность групп после сокращения (ωn). Для пересчета используются формулы, аналогичные следующей: I17{=ROUNDDOWN(IF($H17>K;B17*K/$H17;B17)+RAND();)}. В функции IF определяется, выполняется ли условие H17>K, то есть превышает ли суммарная численность особей емкость среды. Если не превышает, сокращение не производится. Если превышает, численность всех групп сокращается пропорционально, домножаясь на величину K/H17, то есть на разность от деления емкости среды на имеющуюся численность. К примеру, если численность особей превосходит емкость среды в два раза, то K/H17=1/2 и численность всех групп особей уменьшается вдвое. Вычисленная с помощью функции IF величина подвергается вероятностному округлению с помощью функций ROUNDDOWN и RAND, как это описывалось в пункте 3.6.

В столбцах O и P происходит вычисление общего количества половозрелых самок и самцов, ωN(F) и ωN(M) соответственно. Формулы очень просты: O17{=SUM(J17:K17)} и P17{=SUM(M17:N17)}, то есть складывается численность половозрелых самок и самцов всех возрастов.

В столбцах Q и R вычисляется количество пар с самками второго и третьего возрастов, ωN(PF2) и ωN(PF3) соответственно. Пример формулы для таких расчетов виден на рис. 5.1 в строке формул. В данном случае предусмотрено, что самцов может не хватить на всех самок, и вычисления проводятся так: Q17{=IF(P17=0;0;ROUNDDOWN(IF(O17/P17>F_M;P17*F_M*J17/O17;J17);))} и R17{=IF(P17=0;0;ROUNDDOWN(IF(O17/P17>F_M;P17*F_M*K17/O17;K17);))}. Первая из функций IF предохраняет от деления на 0, поскольку далее необходимо делить численность самок на численность самцов. Если самцов в популяции нет, то и пар нет. Если самцы есть, то надо определить, достаточна ли их численность для того, чтобы они покрыли всех самок. Если отношение числа самок к числу самцов меньше величины F_M, то численность пар с самками второго возраста равна численности самок второго возраста, а численность пар самок третьего возраста равна численности самок третьего возраста. Если самцов недостает (надо сказать, что при предусмотренных в модели условиях это может произойти только на самом первом этапе, если задать избыток самок и недостаток самцов), то численность пар самок обоих возрастов сокращается пропорционально наличествующему количеству самцов: количество_самцов * граничное_отношение_самок_к_самцам * количество_самок_соответствующего_возраста / количество_самок, или для самок второго возраста, P17*F_M*J17/O17.

Наконец, в столбце S вычисляется количество потомства (обозначим его как B). Для этого к произведению количества пар с самками второго возраста на плодовитость самок второго возраста следует добавить произведение количества пар с самками третьего возраста на плодовитость самок третьего возраста S17{=Q17*b_f2+R17*b_f3}. После этого все вычисления данного цикла закончены, можно переходить к следующему циклу и выполнять ту же последовательность действий опять.

4.3. Анализ полученной модели

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

Можно ли вычислить соотношение, изменение которого вызывает биения численности? Можно ли построить график такой величины? Можно ли наложить график такой величины на график динамики численности?

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

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

Программа дальнейших действий начинает вырисовываться...

4.4. Выяснение причины колебаний численности

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

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

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

При тех входных параметрах, которые отражены на рисунке 4.2, в первые "годы" существования модельной популяции ее рост заторможен недостатком самцов. "Биения" начинаются примерно с 5-го "года" (т.е. цикла модели). Введем в одну из ячеек сбоку от зоны вычислений формулу, которую видно на рис. 4.3. Над ней сделаем подпись, поясняющую, что именно мы вычисляем. Выберем для этой ячейки процентный формат данных. Затем "растянем" ячейку с формулой и ячейку с заголовком на две ячейки вправо. Calc сам поймет, что в новых ячейках с подписями надо заменить единицу на двойку и тройку. Затем "растянем" ячейки с подписями вниз (до конца модели или на ее начальную часть, на протяжении которой сильнее всего проявляются биения. Отразим полученные данные на графике. Результат показан на рис. 4.3.

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

Интересный способ представления тех же результатов показан на рис. 4.4. Данные столбца, описывающего долю молодых особей, отложены по оси ординат, а данные о доли трехлеток — на оси ординат. Ломаная кривая описывает траекторию системы в показанном фазовом пространстве. Самая первая точка — в нижнем правом углу приведенного графика (100% молодых). Следующая — в начале координат (все особи — двухлетки). На третьем году соотношение поколений 50%/50% (особи первого поколения оставили потомство, и на каждую пару родителей пришлось, в соответствии с показанными на рис. 5.3 параметрами, по двое потомков). 

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

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

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

4.5. Модель с различиями групп особей по выживаемости

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

Обратите внимание: задавать смертность (или обратную ей величину — выживаемость) можно по-разному. Мы будем использовать следующий вариант. Смертность 0,2 означает, что на каждом шаге гибнет 20% особей. Такой смертности аналогична выживаемость (s, от survival) 0,8. В наборе параметров, который мы будем использовать, задается именно выживаемость, а не смертность.

Перестройте предыдущую модель следующим образом. В число входных данных моделей введите показатели выживаемости для каждой из шести рассматриваемых в модели групп. Так, s(f1) — это доля самок первого возраста (потомства на предыдущем шаге модели), которые сохранятся после сокращения численности. Аналогично следует определить s(f2), s(f3), s(m1), s(m2) и s(m3). Поскольку полоса прокрутки не позволяет задавать дробные значения, используйте просто ячейки, значения в которых будут редактироваться вручную. Конечно, можно поступить и иначе: с помощью полосы прокрутки изменять целое число (ячейку с этим числом можно спрятать под саму полосу), а считывать значение выживаемости из ячейки, в которой число, заданное с помощью полосы прокрутки, будет уменьшаться, допустим, в 100 раз.

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

В этой модели, как и в предыдущей, придется, кроме прочего, в числе входных параметров задать b(f2) и b(f3) — плодовитость второго и третьего возрастов, а также K — емкость среды, т.е. такое количество особей, которое может прокормить данное местообитание.

Численность представителей всех групп после сокращения (желтые ячейки в модели) следует вычислять так: ωn(f1)=αn(f1)×s(f1) (используются те же обозначения, что и в столбцах предыдущей модели). Чтобы эта численность не оказалась нецелой, можно использовать вероятностное округление. Используя комбинированную форму записи — с символьным обозначением величин в сочетании с функциями Calc, можно предложить такой вариант: ωn(f1)=ROUNDDOWN(αn(f1)×s(f1)+RAND();). Кстати, в этом пособии мы в дальнейшем будем использовать именно такую форму записи.

В предыдущей версии модели общая численность особей (αN) вычислялась до сокращения численности (так как эта величина использовалась при самом сокращении). Теперь надо вычислять ее после сокращения, то есть ωN=ωn(f1)+ωn(f2)+ωn(f3)+ωn(m1)+ωn(m2)+ωn(m3).

Способы вычисления ωN(F) и ωN(M), т.е. общего количества половозрелых самцов и самок после сокращения численности, а также ωN(PF2) и ωN(PF3), т.е. численность приплодов от самок второго и третьего возрастов можно не изменять. Осталось понять, как следует вычислять количество потомков.

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

n(f0)+n(mo)=ROUNDDOWN(IF((ωN+ωN(P2)×b(f2)+ωN(PF3)×b(f3))<K;ωN(PF2)×b(f2)+ωN(PF3)×b(f3);K-ωN)+RAND();).

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