CC BY
166
Панъгина Нина Николаевна
МОДЕЛИРОВАНИЕ ДИНАМИКИ ЧИСЛЕННОСТИ ПОПУЛЯЦИИ
Любая экосистема состоит из популяций - сообществ организмов одного вида в выделенном пространстве и взаимодействующих между собой и окружающей средой.
Описание популяции обычно представляется через ее численность (биомассу). Это величина целочисленная, но мы будем описывать ее непрерывной функцией от времени - Х((). В качестве функции Х^) можно выбрать плотность популяции (количество особей на определенной площади выделенного ареала).
В своей работе я всегда пытался объединить истину с красотой.
Герман Вейль
Запишем закон изменения численности популяции за время Аt:
АХ = Прирост - Потери Прирост = Я X
Я - коэффициент прироста. Считая, что численность популяции, заполняющей данную экологическую нишу, определяется ограниченным пищевым ресурсом среды Q, полагают, что зависящий от размера популяции коэффициент прироста Я пропорционален величине, характеризующей потребление ресурса Q - X. Полагаем также, что
Рисунок 1.
жасиайемл ссойашй шшрлщий...
-_. т:> -и-.
■ :■■...........
иооотн*
'1 КЧВДГттНг/.
■ ■ I ■!■ ■ - 1 м л
■
1 ■ . _ -и^ : >г ■ | ришщ
потери составляет только естественная убыль особей, которая пропорциональна численности популяции. Следовательно, имеем (уравнение Ферхюльста):
AX = к ■ (Q - X) ■ X - d X (1)
Смоделируем, как изменяется плотность популяции в зависимости от единицы времени time с помощью электронной таблицы Excel. На рисунке 1 представлен экран рабочего листа Excel. Исходные данные: константы Q, к, d (ячейки A5, B5, C5), малая начальная плотность популяции X(0) (ячейка E5) хотя и выбраны условно, но они отражают характер происходящих процессов. Необходимые формулы, коэффициенты, график изменения плотности популяции также представлены на экранном листе. При построении модели используются следующие элементы Excel:
- форматирование ячеек (заливка, объединение, прорисовка границ, выбор шрифтов, выравнивание);
- графическое представление данных (выбор данных и типа графика, оформление графика);
- различный вид формул для задания данных;
- элементы управления (счетчики);
- оформление рабочего листа (примечания, размеры ячеек и др.).
Рассмотрим, как изменяется рост популяции в зависимости от входных параметров. Для варьирования значений коэффициентов Q, к, d используются счетчики, связанные текущими значениями с ячейками B11, B15, B19. Для простоты будем менять лишь один параметр к (влияние остальных - для самостоятельных исследований и выводов).
При увеличении рождаемости (с шагом 0.001 и начальным значением 0.04) популяция, очевидно, быстрее (во времени) достигает асимптотического (равновесного) значения (модель Bio_1.xls).
Упражнение 1.
а) Найдите асимптотическое (во времени) значение численности популяции.
б) При каких значениях к плотность популяции остается постоянной, равной из-
Рисунок 2.
начальной плотности? Определить вид графика при меньших значениях к.
Указание. Используйте соотношение (1) при АХ = 0.
Вид графика на рисунке 1 означает, что адаптивность популяции велика и поступление ресурса достаточно. При дальнейшем увеличении параметра (к = 0.305) возникают затухающие колебания численности популяции (рисунок 2). Это характеризует плохо приспособленную популяцию в неблагоприятной среде (бедной ресурсами).
При последующем увеличении параметра к наблюдаются интересные особенности развития процесса. Возникают (при к = 0.396) устойчивые периодические осцилляции около асимптотического значения (рисунок 3), а также (например, при к = 0.420) «непрогнозируемые» значения численности (рисунок 4).
На рисунке 5 представлена так называемая бифуркационная диаграмма развития популяции в зависимости от параметра к. Для каждого значения параметра (в соответствии с размером окна просмотра) наносятся на плоскость рисунка 120 итераций значения X. Первые 500 точек итера-
Рисунок 3.
Паньгина H.H.
График изменении .4
-1 It U № 4L м
Рисунок 4.
ций переходного процесса пропускаются. Вначале имеем одну точку асимптотического значения численности популяции, затем 2, 4, 8 и т. д., вплоть до области хаоса, когда точки заполняют целые полосы или, наоборот, образуют пустые окна. Немыслимо, но совершенно непредсказуемые значения возникают из полностью детерминированного процесса! (Можно привести и обратный пример, когда из хаотичных аплодисментов возникает бурная овация). Изучением подобных явлений занимается наука синергетика при исследовании различных процессов в самоорганизующихся системах [1].
Данная диаграмма имеет фрактальный характер [2] (части рисунка имеют подобную структуру). На рисунке 6 показана ук-
рупненным планом выделенная рамкой часть основного рисунка. Подобную картину можно наблюдать на большой глубине вложенности частей рисунка.
Упражнение 2. Составьте программу для графического изображения бифуркационной диаграммы (рисунок 5).
Указание.
a) Используйте преобразование координат, чтобы вписать точку X в графическое окно
Yw = Yw . + (Yw - Yw .)(X - X .)/( X - X.),
mm 4 max mm/4 mm7 4 max mm"
где Yw - координата графического окна, соответствующая значению плотности популяции X.
b) Преобразование видоизменяется, если графическое окно использует экранные координаты
Упражнение 3. Дополните программу (из упражнения 2) функцией указателя координат (в исходной системе k-X) любой точки окна и функцией увеличения любой прямоугольной подобласти.
Указание.
a) Для выбранной подобласти соответственно изменяется шаг приращения параметра k.
...w&eftMeftfta ftmfcefctcofifeMUe ¿ЯлгеНсф ßtyjHutcAioeH «f палЛос&ш фе&ержиЛироёаЯЛага арвфссл!
Рисунок 5.
В приложении на CD приведен пример создания программы в среде Visual Basic (Bifurcation.doc), выполняемый модуль, созданный средствами Delphi (bio.exe).
Процессы, описываемые уравнением (1), универсальны в том смысле, что находят свое приложение для более широкого круга явлений, например, в исследованиях по лазерной физике, автокаталитических реакций в химии. Математики открыли универсальное число (подобное числу л), оно называется числом Фейгенбаума (5) и определяется как предел последовательности
5 = (кг - кг-1)Ккг+1 - к)
где ki — значение параметра, при котором асимптотика из 2' точек переходит в 2'+1 точки.
Упражнение 4.
Используя программу (приложения), найдите значения ki удвоения периода для начальных ' и определите приближение числа Фейгенбаума.
Указание. 5 = 4.669...
КОНКУРЕНЦИЯ И СОСУЩЕСТВОВАНИЕ ВИДОВ
Рассмотрим, как «работает» уравнение (1) для случая появления в данном ареале другого вида, потребляющего тот же источник пищи Q. Уравнение (1) записывает-
Рисунок 6.
ся для каждого вида со своими (свойственными виду) коэффициентами:
АХг = кх ■ Хх ■ ^ - Хх- Х2) - йх ■ Хх,
= k9 • • (Q2 — Xi— X2) — d • X2,
(2)
"-2 ^2 ^2 ^2> 2 2' причем Q1 = Q2 = Q, а выражение ^ - Х) уравнения (1) заменяется на Ш-Х1 - Х2) в уравнениях (2), где Х1 и Х2 - численности особей соответствующего вида популяции. (В общем случае, потребление ресурса отличается для различных видов, и надо использовать выражение ^ - а1Х1 - а2Х2), но мы не станем усложнять исследование дополнительными параметрами). В этих условиях выживает только один вид (рисунок 7, модель Bio_2.xls), а другой вымирает, так как вид с большей константой размноже-
W айих условиях вифивгеЛ ¿мошну ofuft вир...
Рисунок 7.
Паньгина Н.Н.
Таблица 1.
Вид K Q d
1 0.22 6 0.2
2 0.2 4 0.1
ния (например к2) поедает пищу гораздо быстрее, и, как следствие этого, даже при малой доле начальной численности вытесняет другой вид (как случилось, например, при заселении Австралии кроликами). Сосуществование видов становится возможным, если количества подводимой пищи различны для разных видов Ф Q2). Проведем адаптацию популяций путем подбора индивидуальных констант к, Q, й. На рисунке 8 приведены графики численности конкурирующих популяций для параметров модели, представленных в таблице 1 (модель Bio_2C.xls).
Проведенное моделирование показывает, почему для выживания видов важную роль играют так называемые экологические ниши и почему выживающие виды столь специализированы. Хорошо известный пример сосуществования и конкуренции - распределение растительности по высоте в горах.
Упражнение 5. (Для самостоятельного исследования).
Смоделируйте симбиоз двух популяций, когда кооперация различных видов облегчает их существование (например, деревья - пчелы). Найдите условие устойчивости видов (постоянство количества особей).
Указание.
а) Предполагается, что скорость размножения одного вида линейно зависит от количества особей другого.
~PJ4*KH пыепенлп И
О Я « Н
Рисунок 8.
б) Если в уравнении (1) не учитываются члены, описывающие внутривидовое подавление - kX2, то образующие уравнения имеют вид:
DX1 = (k1 + k\ ■ X2) • X1 - d1 ■ X1 DX2 = (k2 + k2 ■ X1) ■ X2 - d2 ■ X2
МОДЕЛЬ СИСТЕМЫ ХИЩНИК-ЖЕРТВА
Из старинных охотничьих книг, которые велись со скрупулезной педантичностью, была отмечена временная осцилляция количества зайцев и рысей в заданном ареале (в пересчете к охотничьим трофеям). Смоделируем отмеченную динамику.
Пусть вид-жертва имеет индекс 1. При достаточном количестве пищи данный вид размножается по закону:
Приросту = k1 ■ X1 Скорость убывания, очевидно, пропорциональна числу жертв и хищников: Потери1 = d1 ■ X1 ■ X2 Рассмотрим уравнения для хищников (индекс 2). Так как хищники живут за счет жертв, то скорость размножения хищников пропорциональна их собственному числу и количеству пищи, которое определяется численностью вида 1, а естественная убыль хищников пропорциональна их количеству: Прирост2 = k2 ■ X1 ■ X2 Потери2 = d2 ■ X2 Таким образом, уравнения имеют вид (модель Лотки-Вольтерра):
DX1 = k1 ■ X1 - d1 ■ X1 ■ X2, (3) DX2 = k2 ■ X1 ■ X2 - d2 ■ X2.
Разработка модели в Excel представлена на рисунок 9.
Исходная модель чрезмерно общая и потому очень неустойчива по отношению к входным параметрам. Биологические связи между популяциями более сложны. Например, биологи указывают на возможность выживания жертв, когда определенное минимальное число их выживает (они могут перебраться в другие районы, недоступные хищникам). Данный факт может быть уч-
Рисунок 9.
...биологи ука^-иватй На вурмлфКос&ь вифиваНия
тен в модели добавлением постоянного члена V1 в правой части уравнения для определения прироста вида-жертвы. На рисунке 9 представлен график изменения во времени численностей популяций при указанных на том же рисунке параметрах и начальных численностях X2(0) и X1(0) хищников и их жертв. Отметим запаздывающий характер динамики численности вида-хищника по отношению к численности вида-жертвы (модель Bio_2X.xls).
Приведенные в статье модели в табличном представлении MS Excel приложены на CD.
ЗАДАНИЯ ДЛЯ САМОСТОЯТЕЛЬНОЙ РАБОТЫ
Задание 1
Нарисуйте примерный график временной зависимости суммарной протяженности железных дорог на территории отдельного развитого государства (например, США) в начале эпохи развития этого вида транспорта. Приведите обоснование вида графика.
Указание. Начальный этап развития железнодорожного транспорта имеет вид, подобный логистической кривой (рисунок 10), так как чем больше дорог, тем быстрее
они способствуют строительству новых дорог. Спад темпов строительства вызывается уменьшением территории, где железные дороги еще не построены.
Задание 2
Приведите простую модель развития эпидемии вируса гриппа в большом городе без профилактического и врачебного вмешательства. Смоделируйте в Excel процесс во времени. Выделите временные участки начального распространения заболевания, эпидемического криза и спада эпидемии. Объясните объективное наличие этих фаз даже при проведении (каких?) профилактических мероприятий. Почему вводится жесткий карантин службами МЧС даже при малых зафиксированных случаях заражения опасными вирусами?
Рисунок 10. Логистическая кривая.
Паньгина H.H.
Рисунок 11.
Указание. Если N - количество населения в городе, 1(г) - число инфицированных («вирусоносителей») в момент времени г, то каждый потенциально может передать вирус части здорового населения ^ -I). Рассмотрим относительные величины, то есть, например, г = I / N - доля инфицированного населения, я = 1 - г - доля, подверженная возможному инфицированию (уязвимая к вирусу), тогда доля вновь заболевших за единицу времени (Аг) станет
Аг = к ■ г ■ я, где к - контактный показатель, учитывающий среднее количество контактов одного человека за единицу времени и риск передачи вируса в каждом из этих контактов.
Условные начальные данные для задания:
п = 1, г'(0) = 2.0е - 4 (к примеру, инфицировано 10 человек в городе из 50000 человек),
к = 0.15,
Аг = 1 «условная» неделя болезни.
В данном случае развитие эпидемии также прослеживается «логистической» кривой вида
Литература
i(t) = 1/(1 + exp(-kt)). Задание 3
Классическая модель (SIR) развития эпидемии предусматривает исцеление в условиях задания 2 при наличии антивируса (обзор публикаций в [3]). Предполагая, что увеличение доли (r) исцеленных (и невосприимчивых в дальнейшем к вирусу) индивидуумов происходит пропорционально (коэффициент g) доли инфицированного населения (антивирус активизируется именно в этой среде), получим следующие образующие уравнения
As = -k•i • s Ai = k • i • s - g • i Первое уравнение характеризует уменьшение (знак минус в правой части) доли населения, которая может быть подвержена инфекции в единичный интервал времени. Второе уравнение показывает увеличение на такую же величину доли инфицированных людей за вычетом части инфицированных людей, которые выздоровели через некоторое время и приобрели стойкий иммунитет к данному вирусу (при этом r = 1 - s - i).
Построить данную классическую модель эпидемии с помощью таблиц Excel.
Указание. При условных начальных данных, указанных ниже, получим графики развития эпидемии, представленные на рисунке 11.
i(0) = 0.0002 s(0) = 1 - i(0) k = 0.15 g = 0.3
1. Г. Хакен. Синергетика. М.: «Мир», 1980.
2. Х.-О. Пайтген, П.Х. Рихтер. Красота фракталов. М.: «Мир», 1993.
3. Захарченко А. Черводинамика: причины и следствия http://az 13.mail333.com/4/ch6.html
Паньгина Нина Николаевна, директор Центра информационных технологий г. Сосновыш Бор Ленинградской области.
© Наши а вторы. 2005. Our authors, 2005.
