К ТЕОРИИ НЕОДНОРОДНЫХ ПОПУЛЯЦИЙ [1]

Карев Г.П.

(Москва)

1. Постановка задачи.

Рассмотрим модель динамики численности популяции N(t) в предположении, что коэффициент размножения F(N,a) зависит от некоторого параметра a:

dN/dt = N F(N,a). (1)

Как правило, оказывается недостаточным изучить решения уравнения (1) при каком-либо фиксированном значении параметра. Вариация параметра (в области допустимых значений) может привести к перестройкам качественного поведения решений уравнения (1), которые изучает теория бифуркаций. Известны различные обобщения модели (1) на случаи, когда значение параметра является не фиксированным, а распределенным: а) во времени, б) в пространстве, в) по особям популяции. Чаще всего предполагается, что параметр возмущается некоторой случайной функцией времени типа белого шума. Эти возмущения воздействуют на популяцию в целом, то есть являются одинаковыми для всех особей популяции. Задачи а) и б) обычно сводятся к уравнениям в частных производных и в принципе достаточно изучены.

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

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

2. Модель динамики неоднородной популяции.

Совокупность всех особей популяции в момент t, обладающих заданным значением параметра a, будем называть a-группой; пусть l(t,a) – численность a-группы в момент t. В соответствии с уравнением (1) предположим, что коэффициент размножения a-группы равен F(N,a), то есть зависит от «группового» значения параметра и общей численности популяции, но не зависит от численности a-группы. Динамика такой популяции описывается системой

dl(t,a)/dt = l(t,a) F(N,a),

N(t)= l(t,a)da , (2.1)
l(0, a) = l₀(a),

где A – множество возможных значений параметра. Начальное распределение l(0, a) = l₀(a) параметра a в популяции считается заданным. Очевидно, что при фиксированном значении параметра (то есть при l₀ (a) = d(a-a₀)) модель (2.1) эквивалентна (1.1).

Модель (2.1) в общем виде является достаточно сложным математическим объектом. Далее предполагается, что коэффициент размножения линейным образом зависит от параметра, то есть имеет вид

F(N,a) = f(N)+ ag(N). (2.2)

Этот формально частный случай соответствует большому количеству реальных моделей. Более общий случай рассмотрен в [3].

3. Основные свойства моделей с распределенным линейным параметром.

Обозначим P_t(a) плотность распределения вероятностей параметра a в момент t:

P_t(a) = l(t,a)/N(t).

Математическое ожидание по распределению вероятностей с плотностью P_t(a) будем обозначать E_t, то есть

E_t g = g(a)P_t(a)da .

Обозначим производящую функцию моментов (ПФМ) распределения P_t(a)

M_t(l) =exp(la) P_t (a) da ,l>0 .

Производящая функция моментов (если она существует) взаимно однозначно соответствует своему распределению вероятностей. ПФМ начального распределения M₀ (l) играет ключевую роль при исследовании неоднородной модели (2.1), (2.2).

Введем вспомогательные функции p(t), q(t) как решения уравнений

dp/dt= g(N), p(0)=0, (3.1)
dq/dt=q f(N), q(0) =1.

Справедлива следующая основная теорема.

Теорема 1. Для модели (2.1), (2.2) при всех t

1) текущая численность N(t) удовлетворяет уравнению

dN/dt=NF(N, E_ta) =N[f(N) + E_ta g(N)] (3.2)
и определяется формулой

N(t) = N(0)q(t)M₀(p(t)); (3.3)
2) текущее среднее значение параметра E_t a определяется формулой

E_ta = d[lnM₀(l)/dl?_l=p(t)] (3.4)
и удовлетворяет уравнению

dE_ta/dt = g(N) s²(t); (3.5)
3) текущая дисперсия s (t) определяется формулой

s²(t) = d²M₀(l)/dl²?_l=p(t)]/M₀(p(t)) – ([d M₀ (l)/dl?_l=p(t)]/M₀(p(t))²_. (3.6)

Следствие 1. Плотность распределения вероятностей параметра a в момент t равна

P_t(a) = exp((p(t) a) l₀(a)/N(0)M₀(p(t)). (3.7)

Теорема 1 содержит полное описание модели (2.1)-(2.2) как в виде явных формул для общей численности, среднего и дисперсии параметра и т.д. в любой момент времени, так и виде системы дифференциальных уравнений. Основной шаг состоит в замене фиксированного значения параметра в исходной модели на текущее среднее значение параметра E_ta, которое зависит как от начального распределения, так и от коэффициента размножения модели f(N,a).Теорема 1 и следствие 1 позволяют в принципе решить эту задачу для любых конкретных распределений с известной ПФМ.

4. Важнейшие случаи начальных распределений.

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

Предложение 1. При всех t

M_t(l) = M_o(l+p(t))/M_o(p(t)). (4.1)

Теорема 2. Пусть начальное распределение P₀ (a) является

1) нормальным со средним a₀ и дисперсией s₀. Тогда при всех t текущее распределение параметра P_t(a) будет также нормальным со средним

E_ta = a₀ + s₀²p(t) (4.2)
и той же дисперсией s₀². При этом

N(t) = N₀ q(t) exp(a₀ p(t) + p²(t) s₀²/2); (4.3)
2) пуассоновским со средним a₀. Тогда при всех t текущее распределение параметра P_t(a) будет также пуассоновским со средним E_ta = a₀ exp(p(t)), при этом

N(t) = N₀ q(t) exp(a₀(e^p(t) -1); (4.4)
3) G -распределением с коэффициентами s, k, b, то есть

P₀(a)=l₀(a)/N₀ = s ^k (a – b)^k-1exp[-(a – b) s ]/G( k), (4.5)
где G(k) – G-функция, a ³ b, -ђ < b < ђ, s, k>0. Тогда в любой момент времени t<T₁ = inf{t>0: p(t) = s} параметр a будет также иметь G- распределение с коэффициентами (s- p(t)), k, b, так что

E_ta = b + k/( s – p(t)), s² (t) = k /(s- p(t))². (4.6)
Текущая численность N(t) удовлетворяет уравнению

dN/dt = N[f(N)+ g(N)( b + k /(s- p(t)) (4.7)
и определяется формулой

N(t) = N₀ q(t) exp(p(t)b) /(1- p(t)/s) ^k (4.8)

Из (4.6)-(4.8) следует, что модель (2.1), (2.2) с начальным G-распределением параметра обладает неожиданным свойством: численность (а также среднее и дисперсия) могут «уйти на бесконечность» в момент T₁: p(T₁)=s, и модель теряет смысл при t³T₁.

Список интересных начальных распределений параметра, которые встречаются в реальных ситуациях, может быть продолжен. Так, модель с лог-нормальным начальным распределением «уходит на бесконечность» в момент T₂ = inf{t>0: p(t)>0}. Можно показать, что, например, равномерное начальное распределение в заданном интервале в любой момент времени t>0 становится экспоненциальным в этом же интервале. Это означает, что равномерное распределение параметра не может возникнуть и сохраниться в результате естественного развития популяции (описываемой моделью (2.1), (2.2)), но может быть лишь искусственно созданной «стартовой позицией», сохранение которой требует постоянного внешнего вмешательства в естественную динамику системы. (Это обстоятельство вызывает очевидные социально-экономические аналогии).

5. Простейшие классы неоднородных моделей.

5.1. Неоднородная модель Мальтуса. Качественно новые свойства появляются уже в простейшей мальтузианской модели dN/dt=aN, если каждая особь имеет собственный коэффициент размножения a. Переменные p(t), q(t) для этой модели равны соответственно p(t)=t, q(t)=1. В соответствии с теоремой 1,

N(t) = N₀exp(at)P₀(a)da; dN/dt= E_ta N. (5.1)

Из (5.1) видно, что N(–t) совпадает с преобразованием Лапласа распределения вероятности P₀(a). Это позволяет полностью описать множество всех функций N(t), которые могут быть решением неоднородной модели Мальтуса. Напомним

Определение. Заданная на [0,¥) функция j(l) называется вполне монотонной, если она имеет производные всех порядков jⁿ(l) и при всех n

(-1)ⁿjⁿ(l)³0, l>0.

Заданная на [0,¥) функция j(l) называется абсолютно монотонной, если она имеет производные всех порядков jⁿ(l) и при всех n jⁿ(l)³0, l>0.

Используя классическую теорему Бернштейна, получаем следующее

Предложение 2. 1) Любая абсолютно монотонная функция N(t) является решением некоторой неоднородной мальтузианской модели dN/dt=aN, и обратно. 2) Любая вполне монотонная функция N(t) является решением некоторой неоднородной мальтузианской модели вымирания dN/dt= – a N, и обратно.

Это предложение демонстрирует широту возможных применений неоднородных мальтузианских моделей.

5.1.1. Пусть, например, начальное распределение параметра a нормально со средним a₀ и дисперсией s₀. Тогда, согласно теореме 2, распределение параметра в любой момент времени будет также нормальным со средним E_ta = a₀ + s₀²t и той же дисперсией s₀. Текущая численность определяется формулой

N(t) = N₀ exp(a₀t+t²s₀²/2) (5.3)

и удовлетворяет уравнению

dN/dt = N[a₀ + s₀²t]. (5.4)

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

5.1.2. В приложениях часто встречается экспоненциальное распределение, которое является частным случаем G-распределения при k=1 то есть

P₀(a)= sexp[-(a – b)s], a ³ b, s>0.

Если начальное распределение мальтузианского параметра является экспоненциальным с коэффициентами b, s, то согласно теореме 2, в любой момент t<s распределение параметра будет также экспоненциальным с коэффициентами b, (s‑t). При этом

N(t) = N₀ s exp(tb) /(s – t), (5.5)
E_ta = b+1/(s- t), (5.6)
s(t) =1/(s- t) . (5.7)

Численность популяции определяется также уравнением

dN/dt = N[b+1/(s-t)]. (5.8)

Отсюда следует важный вывод: даже сколь угодно малая, но ненулевая дисперсия в экспоненциальном распределении коэффициента размножения приводит к тому, что модель «взрывается», то есть численность популяции становится бесконечной в конечный момент времени T=s. На первый взгляд это означает, что мальтузианская модель с экспоненциально распределенным коэффициентом роста имеет мало смысла. Однако именно эта модель имеет прямое отношение к моделям глобальной демографии (см.п.6). Существенно, что явление «взрыва» исчезает в модели с экспоненциальным распределением в ограниченном интервале.

Предложение 3. Пусть параметр a принимает значения из интервала А = [b,с] и распределен в начальный момент времени по экспоненциальному закону

P₀(a) = Cexp(-sa), (5.9)
где s =const, C=s/(exp(sc)-exp(sb)) – нормирующая постоянная.

Тогда в любой момент t>0 распределение мальтузианского параметра будет также экспоненциальным в интервале [b,c] с коэффициентом s-t, то есть

P_t(a)=C(t)exp(-(s-t)a), (5.10)
где C(t) =(s-t)/[exp((s-t)с) – exp((s-t)b)]. При этом

E_ta = b+1/(s- t) +(c-b)/ [1-exp((c-b)(s – t)]. (5.11)
Текущая численность N(t) удовлетворяет уравнению

dN/dt =N E_ta (5.12)
и определяется формулой

N(t) =N₀ s exp(bt)[exp(s(c-b)) – exp(c-b) t)]/[(s- t)(exp(s(c-b)) -1)]. (5.13)

Подчеркнем, что модель с экспоненциальным распределением в ограниченном интервале (5.9) качественно отличается от модели с экспоненциальным распределением (4.5). В последней модели происходит «взрыв» в момент T=s (то есть при t®s), после которого модель теряет смысл. Этого явления не происходит в модели с распределением (5.9). Пользуясь формулами (5.11), (5.13), можно показать, что

N(t) ® N₀ s (c-b)exp(sc)/[exp(s(c- b)-1], E_ta® (c+b)/2 при t ® s. (5.14)

Модель сохраняет смысл и после момента t=T, в частности, при t® ¥ E_ta® с.

5.2. Неоднородная логистическая модель задается уравнением

dN/dt = bN(a – N), (5.15)

Оба параметра, от которых зависит коэффициент размножения – a и b – могут быть распределенными. Пусть распределенным является параметр a; это означает, что каждая «a-группа» имеет свою собственную верхнюю границу численности. В этом случае вспомогательные переменные определяются уравнениями

dq/dt=-qbN, q(0)=1,

dp/dt=b, p(0) = 0,

то есть p(t) = bt. Численность неоднородной модели удовлетворяет уравнению

dN/dt = bN(E_ta – N). (5.16)

Это – уравнение Бернулли, и при заданном E_ta решение (5.16) выписывается в явном виде. Существенно, что асимптотика решения этого уравнения заранее известна и равна E_ta. Это свойство уравнения (5.16) позволяет строить неоднородную логистическую модель с желаемой асимптотикой.

Например, N(t)~bt тогда и только тогда, когда начальное распределение параметра a является нормальным. Действительно, в этом случае распределение параметра в любой момент времени будет нормальным со средним E_ta = a₀ + s₀²bt, и

dN/dt = bN[a₀ + s₀²bt -N] (5.17)

Обратное утверждение вытекает из (3.4) и вида ПФМ для нормального распределения.

6. Примеры.

6.1. Модели самоизреживания неоднородных древостоев. Значительная часть известных, внешне совершенно различных полуэмпирических формул изреживания популяций деревьев являются частными случаями модели Мальтуса dN/dt=-aсN с распределенным параметром a (c=const – масштабный множитель). Пусть, например, коэффициент смертности a имеет пуассоновское распределение. Тогда текущая численность N(t) определяется формулой

N(t) = N(0) exp(-a₀(1-exp(-ct)). (6.1)

Это – известная функция Гомпертца, которая совпадает ([4]) с формулой самоизреживания Хильми, выведенной из совершенно других соображений.

Предположение о пуассоновском распределении параметра a соответствует разбиению популяции на счетное число i=0,1,... групп, так что в каждой группе коэффициент смертности равен ic. «Биологический» смысл имеет разбиение лишь на конечное число групп, что отвечает «усеченному» пуассоновскому распределению. Пусть параметр a принимает лишь конечное число значений i=0,1,...k с вероятностями

P₀(a=i) = l₀(i)/N(0) = C₀ (k) (a₀)ⁱ /i!,

где C₀ (k) =[(a₀)ⁱ /i!]^-1 (6.2)

Тогда распределение параметра остается усеченным пуассоновским в любой момент t с заменой a₀ на exp(-ct)a₀ .Численность популяции с начальным распределением (6.2)

N(t)= N(0)C₀ (k)(a a₀ exp(-ct))ⁱ /i! (6.3)

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

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

dN/dt = N (l – N)(N – m), l>m>0. (6.4)

Агрегированное влияние условий местообитания может быть учтено аддитивным параметром a [1].

Рассмотрим модификацию модели (6.4) с отрицательной обратной связью между параметром a и «слишком большой» численностью

dN/dt = N F(N, a)TN(l – N)(N – m) + aN(l – N) (6.5)

Пусть начальное распределение параметра a нормально. Тогда

dN/dt = N (l-N)[N – m + a₀ + s₀²p], (6.6)
dp/dt= l-N, p(0) = 0.

В этом варианте модели (6.6) возникает режим «восстановления численности» исчезающего вида, если только s₀ ¹0, который не существует в исходной модели. Если N(0)<m–a₀ , то численность вида сначала уменьшается до значений, близких к 0, а затем через интервал времени t»(m–a₀)/s₀² быстро увеличивается до стабильного значения N=l. Таким образом, «восстановление» вида происходит после длительного периода очень низкой численности. Единственной его причиной является неоднородность вида (по условиям существования популяций или, в конечном счете, по коэффициенту размножения особей), приводящая к увеличению среднего значения параметра до некоторого критического уровня, после которого быстро восстанавливается стабильный уровень численность даже из очень низкой начальной численности. При этом период восстановления обратно пропорционален дисперсии параметра («разнообразию» составляющих вид популяций).

Аналогичная концептуальная модель (6.4) применяется при моделировании динамики лесных насекомых-фитофагов [1]. В этом случае переход в метастабильное состояние l, возможный лишь при N(0)>m, интерпретируется как вспышка численности. В модели (6.6) с нормально распределенным параметром, характеризующем внешние условия, вспышка численности рано или поздно происходит, если только s₀²¹0. При этом возникает режим «отложенной вспышки», который аналогичен описанному выше и не существует в исходной модели с фиксированным параметром. Если N(0)<m–a₀, то численность популяции сначала уменьшается до значений, близких к 0, а затем через интервал времени t»(m–a₀)/s₀² быстро увеличивается до метастабильного значения N=l. При этом текущие среднее значение и дисперсия параметра также выходят на стационарные значения, зависящие от коэффициентов модели.

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

6.3. Модели глобальной демографии. Рост народонаселения Земли вплоть до настоящего времени с хорошей точностью описывается гиперболической кривой:

N(t)=C/ (T – t)^a , (6.7)

где а =0.99, T =2026 ([5]) либо а = 1, T =2025 ([4]), C » 2*10¹¹ [2, гл.3].

Эта формула предсказывает «демографический взрыв» в момент T. При а = 1 она является решением квадратичной модели роста

dN/dt = N² /C, (6.8)

в которой индивидуальный коэффициент размножения пропорционален общей численности человечества (что, по-видимому, не имеет «физического» смысла). Теория, основанная на модификации уравнений (6.7)-(6.8), позволяющая исследовать важное явление демографического перехода при приближении к критической дате T₁, развивается в настоящее время в работах С.П. Капицы ([2], гл. 3).

Теория неоднородных моделей популяций обосновывает другую интерпретацию эмпирической формулы (6.7), которая оказывается не равносильной закону квадратического роста (см.[7]). Отсюда, как мы увидим, следует, что «проблемы демографического взрыва» (в том смысле, что N(t)®¥ при t ® T ) даже теоретически не существует.

Предполагая, что коэффициент размножения различен для разных групп народонаселения, будем строить модель глобальной демографии в рамках модели (2.1), (2.2). Очевидно, что формула (6.7) совпадает с (4.8) при b=0, s =T₁, k=a, С =sN₀ и q(t)=1. Поэтому гипербола (6.7) является следствием как квадратического закона роста (6.8), так и мальтузианской модели с G-распределенным (k=0.99) либо экспоненциально распределенным коэффициентом размножения (5.5). Заметим, что в мальтузианской модели (5.5)-(5.8) дата демографического взрыва определяется дисперсией начального распределения: T= s=1/s₀ и поэтому может быть отнесена сколь угодно далеко даже при росте среднего коэффициента размножения, если, начиная с некоторого момента, уменьшить его дисперсию. (Иными словами, если направить усилия не на ограничение рождаемости в слаборазвитых странах, а на выравнивание условий существования различных народов).

Учет в модели условия ограниченности коэффициента размножения приводит к тому, что N(t) будет при всех t конечной величиной, хотя и неограниченно растущей. Точнее, если предположить, что начальное распределение коэффициента размножения является экспоненциальным в ограниченном интервале [0,c], то согласно (5.13)

N(t) =N₀ s [exp(sc) – exp(ct)]/[(s- t)(exp(sc) -1)]. (6.9)

При этом в момент «демографического перехода» T=s, согласно (5.14),

N(s) = N₀ s c /[1- exp(sc)], E_sa = c/2 (6.10)

Поэтому «демографический взрыв» является следствием заведомо неверного предположения, (неявно заложенного в модель (6.4) и явно – в (6.5)) о возможности неограниченных значений коэффициента размножения. Отметим, что даже грубая оценка величины границы с по реальным демографическим данным дает с<0.06, при этом график формулы (6.6) практически совпадает с гиперболой (6.4) вплоть до настоящего времени (t =2000). Последующий переход от мальтузианской к неоднородной логистической модели (5.15) с экспоненциально распределенным мальтузианским параметром приводит к модели с ограниченной предельной численностью. При этом удается получить хорошее согласие модели при соответствующих значениях параметров с некоторыми прогнозами ООН.

Литература.

1. Березовская Ф.С., Исаев А.С., Карев Г.П., Хлебопрос Р.Г.. Роль таксиса в динамике лесных насекомых // Доклады РАН (Общая биология). – 1999.– Т.365. – №3.

2. Капица С.П., Курдюмов С.П., Малинецкий Г.Г.. Синергетика и прогнозы будущего. – М., Наука, 1997.

3. Карев Г.П.. Эффекты неоднородности в динамике популяций. // Доклады РАН (в печати).

4. Кофман Г.Б.. Рост и форма растений. – Новосибирск: Наука, 1986.

5. Forster, von H. et al. Doomsday: Friday, 13 November, A.D. 2026. // Science. – 1960. – V.132. – No.1291.

6. Karev G.. Dynamics of a non-uniform population. // Proc. of International Conference EQUADIFF – 99. – Berlin, 1999.

7. Umpleby S.A. The scientific revolution in demography // Population and Environment: a journal of interdisciplinary studies. – 1990. – V.11. – No.3.

---

[1] Работа поддержана грантами РФФИ 98-01-00483 и 99-04-49300.