EcoSimulation — 03. Усложнение модели: имитация логистического роста

 

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

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

 

3. Усложнение модели: имитация логистического роста

3.1. Сравнение экспоненциального и логистического роста

На прошлом этапе мы сделали модель, описывающую экспоненциальный рост популяции. Как вы помните из курса экологии, усложнением экспоненциальной модели является логистическая (рис. 3.1). Бельгийский математик Пьер Ферхюльст (1804–1849) ввел в уравнение Мальтуса компонент, который обеспечивает торможение популяционного роста по мере приближения к величине K — емкости среды. 

Рис. 3.1. Логистический рост (Шабанов, Кравченко, 2009)

Когда мало, выражение в скобках практически равно единице, и рост мало отличается от экспоненциального; когда N достигает величины K, прирост останавливается, и численность популяции становится стабильной, соответствующей K.

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

Таблица 3.1. Добавление к модели экспоненциального роста ячеек, имитирующих логистической рост

  A B C D
1        
2        
3 N0= 1000    
4 r= 0,5    
5 К= 10000    
6 t N (exp) N (log)  
7 0 1000 {=B3} 1000 {=B3}  
8 1 1500 {=B7+B7*B$4} 1450 {=C7+C7*B$4*(B$5-C7)/B$5}  
9 2 2250 {=B8+B8*B$4}    
10 3 3375 {=B9+B9*B$4}    

 

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

Вводя формулу в ячейку C6 (именно ее потом надо будет «растянуть» на весь столбец), учтите, что ссылка на предыдущую ячейку (C5) должна быть относительной, а на ячейки для ввода параметров (B2 и B3) — абсолютной, по крайней мере, в отношении строки (B$2 и B$3 — этого достаточно, если формула будет «растягиваться» внутри столбца).

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

Рис. 3.2. Два способа сравнения экспоненциального и логистического роста рядом; обратите внимание: ячейка C8 выделена, курсор на строке формул, ячейки, на которые ссылается формула, выделены цветовыми рамочками (при необходимости модель можно скачать с сайта)

 

3.2. Учет смертности в определенном возрасте

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

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

В экспоненциальной модели учет смертности будет не столь обоснован, как в логистической. Экспоненциальная модель построена на утверждении dN/dt=r×N. Мы можем отдельно рассмотреть рождаемость и смертность и учесть, что рождаемость увеличивает численность популяции, а смертность — снижает. Мы можем рассматривать r как разницу между рождаемостью и смертностью.

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

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

Таблица 3.2. Сравнение логистического роста без учета смертности особей и с ее учетом

  A B C D
1        
2        
3 N0= 1000    
4 r= 0,5    
5 К= 10000    
6 t N (логистический) N (логистический со смертью на 4-м шаге)  
7 0 1000 {=B3} 1000 {=B3}  
8 1 1450 {=B7+B7*B$4*(B$5-B7)/B$5} 1450 {=C7+C7*B$4*(B$5-C7)/B$5}  
9 2 2069,9 {=B8+B8*B$4*(B$5-B8)/B$5} 2069,9 {=C8+C8*B$4*(B$5-C8)/B$5}  
10 3 2890,6 {=B9+B9*B$4*(B$5-B9)/B$5} 1890,6 {=C9+C9*B$4*(B$5-C9)/B$5-C7}  
11 4 3918,1 {=B10+B10*B$4*(B$5-B10)/B$5} 2207,2 {=C10+C10*B$4*(B$5-C10)/B$5-C7*B$4*(B$5-C7)/B$5}  

 

В таблице 3.2 расчеты, соответствующие логистическому росту, перенесены в столбец B (в Calc, как вы помните, достаточно выделить столбец и перетащить его). В столбце C построена новую зависимость. 

В ячейках, соответствующих 7-й строке, всюду стоят одни и те же значения. Если бы адресация в ячейке B7 была абсолютной, формулу в ней можно было бы просто «растянуть» «» на ячейку C7. А поскольку адресация в этой ячейке относительная, формулу в ней можно скопировать. Однако, обратите внимание, что если вы просто скопируете ячейку B7 и вставите ее в C7, Calc перепишет содержащуюся в ней формулу (она будет указывать на ячейку C5). Действовать надо иначе: выделить ячейку B7, войти в строку формул, выделить формулу, скопировать ее в буфер (Cntr+C), выйти из режима редактирования формулы с помощью клавиш Enter или Esc, войти в ячейку C7, перейти в строку формул, вставить нужную формулу (Cntr+V) и выйти с помощью клавиши Enter.

Теперь построим нужную зависимость. Формулы в ячейках C8 и C9 аналогичны таковым в B8 и B9. Их нужно не копировать, а растягивать (ссылки на входные данные в этих ячейках носят абсолютный характер, а обращения к предыдущим шагам времени должны адресоваться к ячейкам в столбце C, а не B). В ячейке C10 нужно учесть гибель особей, появившихся в ячейке C7. Начиная с ячейки C11, из формул надо вычитать прирост, который добавлялся три шага назад. Так, из численности особей в ячейке C11 надо вычесть C7*B$4*(B$5-C7)/B$5, то есть прирост в ячейке C8. Затем формулу, которая получилась в ячейке C11, можно «растянуть» на весь столбец.

Исследуйте получившуюся зависимость (рис. 3.3). Вы можете убедиться, что если r достаточно велико (например, r=0,5) численность стабилизируется на уровне, более низком, чем K. Если r невелико (например, r=0,2), численность популяции вначале возрастает, а потом, с появлением смертности, снижается.

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

Рис. 3.3. Мы учли смертность особей по достижении ими определенного возраста. Показана формула в ячейке С11; при необходимости, модель доступна для скачивания

3.3. Численность популяции не может быть отрицательной

Используем только что полученную модель еще для одного эксперимента (не вполне корректного с точки зрения биологии). Что будет, если r примет отрицательное значение (например, r=–0,1)? Биологического смысла тут немного, так как отрицательное значение r может принимать в том случае, если в нем уже учитывается смертность, которую в рассматриваемой модели мы и так уже предусмотрели. Тем не менее рассмотрим и такой случай. Результат показан на рисунке 3.4. 

Рис. 3.4. Может ли численность популяции быть отрицательной?

Модель, не предусматривающая смерть особей через 4 года, ведет себя адекватно: ее численность шаг за шагом снижается (и рано или поздно дойдет до нуля). А в модели с гибелью четырехлетних особей численность популяции на определенном этапе (например, на 3-м шаге в показанном на рис. 3.4 варианте) численность популяции становится отрицательной (а потом опять на какое-то время становится положительной)! Имеет это какой-либо биологический смысл? Вероятно, нет (в каком-то случае мы можем вложить в отрицательную численность какой-то смысл, но, вероятно, это не тот случай). Вероятно, нужно сделать так, чтобы численность популяции в нашей модели не могла принимать отрицательные значения.

Для решения этой проблемы используем одну из функций Calc. До сих пор в формулах мы использовали действия. Кроме того, мы можем вводить в формулы функции, выбирая их из достаточно большого набора. Функции имеют названия, которые набираются прописными буквами, а после названий в них всегда находятся круглые скобки, в которых могут стоять необходимые данные. Calc и, например, Excel имеют совместимые наборы функций. Названия функций могут вводиться и показываться на разных языках. Например, сейчас мы обсуждаем функцию, которая на английском языке называется IF. Русский вариант этого названия — ЕСЛИ. При переходе с одного языка на другой электронные таблицы сами будут переводить названия функций. В этом пособии мы будем использовать английские названия (например, потому, что при написании формул удобно не переключать языки при необходимости использовать символы, присутствующие только в английской раскладке, как, например, $). Для того, чтобы вводить английские названия функций, надо включить их в параметрах Calc. Пройдите по пути «Засоби/Параметри/LibreOffice Calc/Формули», и вы попадете на диалог, показанный на рис. 3.5. Поставьте «флажок» на опции «Використовувати англійські назви функцій». Кстати, на вкладке «Параметри мови» можно выбрать украинский язык (для которого даются приведенные в этом пособии примеры), если в системе по умолчанию установлен какой-то другой.

Рис. 3.5. В этом пособии используются английские названия функций

Формат обсуждаемой нами функции таков: IF(условие;результат_если_условие_истинно;результат_если_условие_ложно). Условием может быть любое равенство или неравенство, использующее знаки =, >, <, >= (больше или равно; не меньше; ⩾; ≥), <= (меньше или равно; не больше; ⩽; ≤). После первой точки с запятой следует указать, что функция вставит в ячейку в случае если равенство или неравенство справедливо, после второй — что должно стоять в ячейке, если равенство или неравенство неверно.

Предположим, в ячейке C8 стоит {=C7+C7*B$4*(B$5-C7)/B$5}. Напишем там {=IF((C7+C7*B$4*(B$5-C7)/B$5)<0;0;C7+C7*B$4*(B$5-C7)/B$5)}. Если результат вычисления для численности популяции окажется отрицательным, эта функция подставит вместо него 0; в противном случае будет вставлен результат вычислений в соответствии с исходной формулой.

Снабжаем подобным усложнением все строки в расчетах, — и результат вычислений никогда не окажется отрицательным. Чтобы сделать это, надо вставить указанную формулу в ячейку C8 и «растянуть» ее на ячейку C9, а затем изменить выражения в ячейках C10 и C11, и, наконец, «растянуть» ячейку C11 до конца столбца.

К сожалению, на этом сложности не заканчиваются. Количество особей, которое умирает на шаге, имитируемом в ячейке C11, высчитывается как {C7*B$4*(B$5-C7)/B$5} — это прирост, который добавился на шаге C8. Если величина r отрицательная, прирост в ячейке C8 тоже отрицательный. Вычитая отрицательною величину, мы получаем положительную. После того, как численность популяции снизилась до 0, она продолжает расти! Естественно, это не соответствует биологической сути имитируемых явлений. Как быть? Вставить еще одну функцию IF в уже имеющуюся! Например, это можно сделать так. Вставим в ячейку C11 выражение {=IF((C10+C10*B$4*(B$5-C10)/B$5-C7*B$4*(B$5-C7)/B$5)<0;0;IF(C7*B$4*(B$5-C7)/B$5>0;C10+C10*B$4*(B$5-C10)/B$5-C7*B$4*(B$5-C7)/B$5;C10+C10*B$4*(B$5-C10)/B$5))}. Его смысл таков: если четыре шага назад численность модельной популяции увеличивалась, имеющаяся численность сокращается на эту величину; если нет  — остается без изменений. Конструкция, которую мы использовали, имела такой вид: IF(условие;результат_если_условие_истинно;IF(условие;результат_если_условие_истинно;результат_если_условие_ложно)).

3.4. Округление до целых

Обратите внимание, что во многих случаях результаты вычислений представлены дробными числами. Эта проблема сложнее предыдущей; дело в том, что в ряде случаев дробные значения могут иметь смысл. Впрочем, для начала определим, как избавиться от нецелых значений.

Во-первых, десятичные дроби можно просто не показывать. Это можно сделать, например, с помощью кнопки «Вилучити десяткове місце», которая вместе с кнопкой «Додати десяткове місце» находится на панели инструментов по умолчанию (рис. 3.6). Другой способ — можно щелкнуть правой клавишей мыши на ячейке, в контекстном меню выбрать «Формат комірок...», а там — первую вкладку, «Числа». В окошке «Десяткові знаки» можно указать количество знаков после запятой. В то же самое диалоговое окно можно попасть в результате выполнения команд «Формат/Формат комірок...». В описанных случаях речь идет не об округлении чисел, а просто о том, сколько знаков будет демонстироваться. Если Calc при этом не и будет показывать десятичные разряды, он все равно будет учитывать их при вычислениях: они никуда не денутся, они просто станут невидимыми.

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

Во-вторых, нецелые значения можно округлить. Для этого можно использовать команды ROUND(аргумент;число разрядов), ROUNDUP(аргумент;число разрядов) и ROUNDDAWN(аргумент;число разрядов). ROUND(аргумент;число разрядов) округляет в соответствии с принятыми правилами округления, ROUNDUP(аргумент;число разрядов) округляет до нижнего значения, а ROUNDDAWN(аргумент;число разрядов) — до верхнего. Работу этих команд можно показать на примере, показанном в таблице 3.3. Числа, стоящие в заголовках строк, в каждом столбце здесь округлены командой, стоящей в заголовке столбца.

Таблица 3.3. Сравнение разных функций Calc, обеспечивающих округление данных 

Число ROUND(число;2) ROUND(число;0) ROUND(число;) ROUNDUP(число;) ROUNDDOWN(число;)
1,345 1,35 1 1 2 1
1,5 1,50 2 2 2 1
1,666 1,67 2 2 2 1

 

Как вы можете увидеть, если не указывать количество разрядов, аргумент будет округлен до целых значений.

Итак, указав для всех строк нашей модели округление, мы получим только целые значения. Лучше всего использовать именно округление, при котором значения, меньшие, чем 0,5 округляются в нижнюю сторону, а большие, чем 0,5 — в верхнюю. Такое округление окажет наименьшее влияние на конечный результат (округление вниз будет его занижать, а вверх — завышать). Приведенная в прошлом подпункте формула для ячейки C11 будет в таком случае выглядеть так: 

C11 {=ROUND(IF((C10+C10*B$4*(B$5-C10)/B$5-C7*B$4*(B$5-C7)/B$5)<0;0;IF(C7*B$4*(B$5-C7)/B$5>0;C10+C10*B$4*(B$5-C10)/B$5-C7*B$4*(B$5-C7)/B$5;C10+C10*B$4*(B$5-C10)/B$5));)}.

3.5. Вероятностное округление до целых

Давайте теперь рассмотрим ситуацию, когда дробные значения имеют биологический смысл. Например, известно, что простая воспроизводящая рождаемость (количество детей на одну мать, при котором численность поколения детей равна численности поколения родителей) составляет для высокоразвитых стран 2,03, а для слаборазвитых — 2,2 (предположим, что для Украины эта величина будет близка к 2,1). Если бы у каждой матери (пары родителей) было ровно 2 потомка (суммарный коэффициент рождаемости=2), неизбежная смертность привела бы к тому, что поколение детей оказалось бы меньше, чем поколение родителей.

А что означает суммарный коэффициент рождаемости 2,1? Матерей, которые рожают одну десятую потомка, к счастью, нет. Но такая рождаемость будет характерна для случая, когда в 9 семьях из 10 будет по 2 ребенка, а в одной — 3 (21 ребенок на 10 матерей). Округляя такую величину с помощью функций ROUND и ROUNDDAWN, мы всегда получим 2, а ROUNDUP — всегда получим 3. А как сделать, чтобы величина 2,1 округлялась до 2-х 9 раз из 10, а до 3-х — 1 раз из 10?

Для этого мы можем воспользоваться функцией RAND(). Эта функция генерирует случайное число, находящееся между 0 и 1. Функция RANDBETWEEN(нижняя_граница;верхняя_граница) задает диапазон, внутри которого должно находиться целое случайное число.

Для нас важно то, что формула {=ROUNDDOWN(2,1+RAND();)} задает именно такую зависимость, которая в 9-ти случаях из 10-ти будет равна 2, и в 1-м — 3. При каждом пересчете случайное число будет принимать новое значение. В 9 случаях из 10 случайная величина будет меньше 0,9, и значение суммы округлится до 2; в 1 случае из 10 значение случайной величины будет находиться между 0,9 и 1, и значение суммы округлится до 3. Запустить такие пересчеты можно, нажимая клавишу F9.

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

C11 {=ROUNDDOWN(IF((C10+C10*B$4*(B$5-C10)/B$5-C7*B$4*(B$5-C7)/B$5)<0;0;IF(C7*B$4*(B$5-C7)/B$5>0;C10+C10*B$4*(B$5-C10)/B$5-C7*B$4*(B$5-C7)/B$5;C10+C10*B$4*(B$5-C10)/B$5))+RAND();)}.

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

Два варианта логистической зависимости, построенные со случайным округлением, приведены в этом файле. Скачайте (а лучше — сами постройте!) его и исследуйте, как меняется итоговый результат при каждом нажатии клавиши F9  или выполнении команды «Дані/Обчислити/Обчислити» (рис. 3.7).

Рис. 3.7. Вероятностное округление логистического роста с учетом и без учета смертности. Обратите внимание, что с этими условиями логистический рост иногда «тормозит», как 12-м шаге модели; рост со смертностью «застревает» на численности в 1 особь. При нажатии на клавишу F9 результаты изменятся

3.6. Пример полученного при моделировании артефакта

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

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

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

Что делать в таком случае? Один вариант — использовать имеющуюся модель только в той части пространства ее возможных состояний, где она адекватна. Второй вариант — усложнять модель, приближая ее к действительности.