Исследование уравнения равновесной плотности биологического вида в пространствах различных размерностей Текст научной статьи по специальности «Математика»

УДК 519.642.4

А. Г. Бодров, А. А. Никитин2

ИССЛЕДОВАНИЕ УРАВНЕНИЯ РАВНОВЕСНОЙ ПЛОТНОСТИ БИОЛОГИЧЕСКОГО ВИДА В ПРОСТРАНСТВАХ РАЗЛИЧНЫХ РАЗМЕРНОСТЕЙ*

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

Ключевые слова: интегральные уравнения, математическая биология, численные методы.

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

Ñ = (b^d)N - dj

Rn

(1)

С(0 = Nbm(0 + b I m(V)C(í + V) dV^(d + dw(0)C(0 - 'Ц^- J w(V)C(V) dV,

Rn 1»

где N(t), C(t,£) — неизвестные функции (аргумент t далее опускаем); b, d, d — неотрицательные константы — параметры модели; го(£), т(£) — плотности вероятностей, которые также являются параметрами модели; 1 ^ п ^ 3. Система (1) описывает динамику одновидовой пространственной биологической модели, предложенной У. Дикманом и Р. Лоу в работах [1, 2]. Подробное описание рассматриваемой здесь биологической модели приведено в статье авторов [3], продолжением которой является данная работа.

Для более полного понимания модели приведем описание ее параметров:

• b, d — темпы рождаемости и смертности соответственно. Средние промежутки времени, в которые индивидуум порождает потомство/умирает, b > d ^ 0;

• Ó 0 и w(0 — сила и ядро конкуренции соответственно; dw(£) описывает в единицу времени интенсивность того, что индивидуум на расстоянии £ убьет другого;

• т(£) — ядро рождаемости. Описывает плотность вероятности порождения нового индивидуума, отстоящего от своего родителя на вектор Í,. Функции т(£), го(£) являются неотрицательными, экспоненциально убывающими на бесконечности нормированными плотностями вероятностей, удовлетворяющими следующим условиям:

•»(£)>«, /•»(£)« = ! . »(«>0, /„(£)« = 1.

Rn 1»

В рамках данной статьи будет рассматриваться радиально-симметричный случай (безусловно, наиболее естественный с биологической точки зрения), в котором функции т(£), го(£) предполагаются радиально-симметричными, что в свою очередь приводит к существованию и единственности радиально-симметричной равновесной функции С(£) (более подробно об этом будет сказано далее);

1 НИУ ВШЭ, студент магистратуры, e-mail: drnobQgmail.com

2 Факультет ВМК МГУ, доц., к.ф.-м.н., e-mail: nikitinQcs.msu.su

* Работа выполнена при поддержке гранта Президента РФ для государственной поддержки молодых российских ученых — кандидатов наук МК-6108.2015.9.

• N — средняя плотность популяции, ожидаемое число индивидуумов на единицу площади, N ^ 0;

• С(£) — ожидаемая попарная плотность, т. е. среднее число пар индивидуумов, отстоящих друг от друга на вектор £ на единицу площади.

Предполагается, что плотности индивидуумов в точках, находящихся на больших расстояниях друг от друга, независимы, т.е. lim С(£) = N2. Помимо этого, из очевидного биологического

—s-oo

смысла следует, что С(£) ^ 0. Принимая это во внимание, заметим, что первое уравнение динамики в (1) представляет собой модифицированное уравнение Ферхюльста N = gN — cN2 (см., например, [4]), в котором внутривидовая конкуренция зависит от расположения индивидуумов в пространстве (здесь параметр д играет роль "естественного прироста" популяции, а с описывает силу внутривидовой конкуренции).

Для того чтобы система в дальнейшем имела более удобный вид, перейдем к нормированной ожидаемой попарной плотности С(£) = C{f;)/N2. Действительно, после такого перехода уравнение динамики N примет вид

N = (Ъ- d)N - dYN2, где Y = j w(£)C

Rn

и его схожесть с уравнением Ферхюльста станет еще более заметна. Заметим, что интеграл Y, описывающий вызванные пространственным расположением индивидуумов изменения в силе конкуренции, зависит только от времени. С таким набором договоренностей перейдем к поиску равновесия в системе.

2. Поиск равновесия. Выпишем равновесную систему уравнений, т. е. систему при = 0 и Г(/.0 = 0:

"0 = (Ь-й)ЛГ-<йЧУ2,

о = ^¡р- + ь I т(г])С(£ + г}) йг)- йСЦ) - МОСЮ - ¿УМС(0-

кп

Нетривиальное решение первого уравнения из этой системы можно легко выразить в явном виде: N = (Ь — <])/<]¥. Теперь приведем равновесное уравнение на подставив в него равновесное Ж:

Кп

Далее будем искать пару из функции С* € Ь 1(КП) и числа У*, в которой С* является решением уравнения (2) при У = У*, а У* = ^ Ц£)С(£)

2.1. Решение уравнения равновесия в одномерном случае. Предполагая симметричность функции т(£), после упрощения (2) получим

(¿40 + (.)£■«) = + Ь[т * С]«). (3)

Здесь символом [тФС1({) обозначена свертка ,ти* функций: / т(,)СК-,) % В несколько „ругой

Кп

формулировке исследованием этого уравнения занимались А. А. Давыдов и др. [5, 6]. Ими был получен ряд интересных результатов*, в частности доказаны теоремы существования и единственности решения уравнения (3).

* Несмотря на то что исследования проводились при параметре Л = 0, было предварительно показано, что решение при Л ф 0 не существует.

2.1.1. Отсутствие решения. Наиболее важным результатом А. А. Давыдова, В. И. Данченко и М. Ю. Звягина, имеющим важное прикладное значение, является установление отсутствия решения уравнения (3) при й ф 0.

Для доказательства этого факта возьмем преобразование Фурье от обеих частей уравнения (3) и получим

¿ш(7(/х) + ЬС{ц) = ^Ут(д) + Ып{ц)С{ц).

Здесь /(/х) — образ Фурье функции /(С)- Рассматривая это уравнение в 0 и учитывая, что «(0) = /ш(() * = 1 (аналогично дл, «,({)), 6уде„ и„„ь

Кп

йУ + ЬС{ 0) = + ЬС{ 0).

Таким образом, после элементарных преобразований получаем йй¥ = 0. Данное равенство равносильно одной из двух альтернатив:

1) справедливо йУ = 0, что влечет за собой отсутствие внутривидовой конкуренции, и достаточ-

но просто как с биологической (экспоненциальный рост популяции), так и с математической (аналитическое решение, подразумевающее при й ф 0 выполнение у}(£)С(£) = 0 почти всюду) точки зрения;

2) й = 0, что подразумевает рассмотрение биологического вида, смертность которого обусловлена

исключительно внутривидовой конкуренцией.

Для второго случая В. И. Данченко и М. Ю. Звягин предложили наборы аналитических функций С(£), ии(!;) и т(£), для которых данное уравнение имеет решение. Однако с биологической точки зрения интерес представляет случай, в котором т(£) ж и)(£) являются плотностями нормальных распределений с дисперсиями от и ои! соответственно. И в этом случае найти аналитическое решение не получается, поэтому приоритетной становится задача поиска численного решения данного уравнения. Данная задача уже рассматривалась авторами в статье [3], где был предложен численный метод для решения уравнения (3) в случае одномерной популяции.

2.1.2. Метод Нистрёма. Пусть дано уравнение А/ = Ь, где А — исследуемый нами линейный интегральный оператор, Ь — известная функция, / — неизвестная функция. Тогда можно дискретно аппроксимировать оператор А и функцию Ь на некоторой сетке О матрицей Аа и вектором Ьа- В результате получим систему линейных алгебраических уравнений Аа/а = Ьа, решая которую с использованием, например, Ы1-разложения найдем вектор /а — дискретную аппроксимацию функции /.

сю

Для аппроксимации несобственного интеграла* вида J /(ж) йх на сетке О следует воспользо-

— сю

ваться тем, что он представим в виде

оо Ь

//(ж) йх = Нт / /(ж) йх.

Ь^ос J

-Ь

Далее воспользуемся тем, что решение рассматриваемого уравнения существует (см. статью [5]). Поэтому подынтегральная функция представляет собой произведение функции, сходящейся на бесконечности к константе, на экспоненциально убывающую функцию. Все это подразумевает, что при достаточно больших Ь выполняется неравенство

1^

/(ж) йх — J /(ж) йх

< £

-ь

Этот интеграл мы будем рассматривать в смысле главного значения.

для любого сколь угодно малого е > 0. После "округления" несобственного интеграла

оо Ь

/(ж) (¿ж и У /(ж) dx

—ос —Ь

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

ь

J /(ж) dx И WiJ{xi),

_L Xi£G

где щ — веса точек сетки ж*. Саму сетку в таком случае следует брать из соображений оптимального покрытия отрезка [—Ь,Ь]. Небольшую проблему в данном случае представляет свертка

т(г])С(£ — г]) йг} как интеграл, зависящий от параметра, поскольку приближение

lj

m(,)е(£-,) А, «/m(,)(?(£-,) А,

—оо —Ь

может оказаться неверным при £ близком к Ь, так как если мы берем интеграл по отрезку [—Ь, Ь] то значения т(г])С(£ — г]) за пределами сегмента игнорируются. В силу того что С(—Ь) и 1 и

-ь о

m(i)e(, + i)A,«I/m(,)A,,

мы получаем значительную погрешность округления. Однако в этом случае можно воспользоваться тем, что lim С(£) = 1, и приблизить свертку как

ICI—s-oo

ос L оо

I шШ(( -,)*,»/ mW (CK - ,) - 1) Л, + / mW л,.

—оо — L — оо

При использовании данного метода и квадратурной формулы трапеций была достигнута относительная поточечная ошибка в 0.5% на наборе тестовых функций, предложенных В. И. Данченко и М. Ю. Звягиным [6]. При увеличении порядка квадратур наблюдалось увеличение ошибки.

3. Многомерный случай. Сложность вычисления решения методом Нистрёма есть 0(|G|3), где — число точек на сетке G. Поскольку для получения приближенного решения с приемлемой точностью требуется экспоненциально увеличивать число точек, то для решения многомерной задачи этот метод будет неэффективным. Сложность решения методом Нистрёма по памяти — 0(|G|2), что также затрудняет вычисления. Однако здесь можно воспользоваться описанным в [3]

оо

методом рядов Неймана ^подразумевающим поиск решения в виде / = ^^ АгЬфактически пред-

г= 1

ставляющим собой метод Якоби (простых итераций), неявно примененный к матрице системы Aq.

Напомним, что все функции в задаче радиально-симметричные, а нетривиальное решение единственно [5]. Это позволяет доказать следующее утверждение.

Лемма. Если функция С (ж) является решением (3), то функция С (ж) радиально-симметрична, т. е.

~ 1 I ~ -ч - ПЖП/2 _

С{х) = ^щу f С(ж ) dx , 5n(|x|) = f dx = p^n ^ W\n >

где 5п(|ж|) — площадь поверхности n-мерной сферы радиуса |ж|.

Доказательство. Сначала покажем, что если С(х) является решением, то решением будет и функция

Сг( х) = <£ C(x)dx.

Sn(\x\) J

Действительно, если взять от исходного уравнения интеграл j) dx , то получим уравнение для Сг(х). Тот факт, что lim Cr(x) = 1, также легко проверяем. Поскольку решение единственно,

\х\ —too

то С(ж) = C'i- (х). Следовательно, можно совершить переход к многомерной сферической системе координат в уравнении (3) и получить уравнение

оо сю

{Ь + Лш{г))С{г) = Ь JM(£,r)C(£)d£ + dm(r) J Щ)ш{£)С{£) di, о о

где S(£) = тт(2£)п~1 при п = 2, 3, а Ш(£,г) — функция, полученная в результате преобразования якобиана, домноженного на — r¡).

В двумерном случае, когда т(£) — плотность многомерного нормального распределения с диагональной матрицей ковариации с о2т на диагонали, ядро Ш(£, г) будет иметь вид

2тг

—,„ Г 1 f (г + Icos op)2 + (I sin op)2 \ ,

о

2тг

1 f г2+£2\ f ( £r cos (p \ , 1 ( r2 + £2\n f£r ■ exp--——— / exp---— dep = --exp--——— j0

Зтго"2а/ я/ \ ®и% / Зтго"2о~^¡^ J с

о

Здесь 30 {£г/а%) — модифицированная функция Бесселя первого рода. В трехмерном случае Ш(£, г) будет иметь вид

Шв , У7^ ( г2 + Г2 ^ 2г/

4. Биологические результаты.

4.1. Мотивация. Биологической мотивацией для использования пространственной модели является принцип конкурентного исключения (принцип Вольтерры-Гаузе): "полные конкуренты не могут сосуществовать бесконечно".

Данный принцип можно обосновать, исходя из модели Лотки-Вольтерры для двух конкурирующих видов (с плотностями N и М соответственно):

Л; 1>х N - /¡хм.ММ. М = ЬмМ - /1м х.ММ.

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

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

Одним из решений данной проблемы является переход к рассмотрению пространственной модели, в которой один вид оказывается более "сильным" в конкурентном плане, а второй способным быстрее открывать новые области (так называемый competition-colonization trade-off).

4.2. Результаты. В данной модели наблюдаются два хорошо известных биологам эффекта кластеризация (clustering) и переразброс (overdispersion).

Первый эффект сводится к тому, что при малых значениях а„, дисперсии распределения т(£) и больших значениях aw дисперсии распределения «;(£) наблюдается склонность популяции к жизни в маленьких группах-кластерах, внутри которых наблюдается высокая конкуренция, что приводит к снижению численности популяции N и ярко выраженному "пику" в области нуля у попарной корреляции С(£).

Второй эффект вызван тем, что при а„, 3i> 0, aw и 0 выживают преимущественно потомки на больших расстояниях, что приводит к практически отсутствующей конкуренции с потомками, увеличению численности популяции N и ярко выраженной "яме" в области нуля у попарной корреляции С(£).

На представленном изображении можно заметить, что эффект переразброса становится более ярко выражен при увеличении размерности, в то время как эффект кластеринга проявляет себя по-разному. Для одномерного случая (рисунок, а, г) можно эмпирически заметить склонность к вымиранию вида при увеличении аи; в случае кластеризации, в двумерном случае (рисунок, ¿>, д) можно заметить стремление к асимптоте сверху, в то время как в трехмерном случае (рисунок, в, е) заметно восстановление популяции после резкого ее уменьшения в самом начале.

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

Модель также зависит от значений параметров (1, и Ь, однако можно легко заметить, что после деления уравнения (3) на Ь этот параметр будет входить в (3) только в составе й/Ь, таким образом.

можно рассматривать только d/b. Его увеличение приводит к тому, что вышеописанные эффекты

становятся более ярко выраженными.

СПИСОК ЛИТЕРАТУРЫ

1. Law R., Dieekmann U. Moment approximations of individual-based models // The Geometry of Ecological Interactions: Simplifying Spatial Complexity. Cambridge: Cambridge University Press, 2000. P. 252-270.

2. Law R., Dieekmann U. Relaxation projections and the method of moments // The Geometry of Ecological Interactions: Simplifying Spatial Complexity. Cambridge: Cambridge University Press, 2000. P. 412-455.

3.Бодров А.Г., Никитин А. А. Качественный и численный анализ интегрального уравнения, возникающего в модели стационарных сообществ // Докл. АН. 2014. 455. № 5. С. 507-511.

4. Verhulst P.F. Notice sur la loi que la population poursuit dans son accroissement // Correspondance Mathematique et Physique. 1838. 10. P. 113-121.

5. Данченко В.И., Давыдов А. А., Никитин А. А. Об интегральном уравнении для стационарных распределений биологических сообществ // Проблемы динамического управления: Сборник научных трудов. Вып. 3. М.: Изд. отдел ф-та ВМК МГУ, 2009. С. 15-29.

6. Давыдов А. А., Данченко В. И., Звягин М.Ю. Существование и единственность стационарного распределения биологического сообщества // Тр. МИАН. 2009. 267. С. 46-55.

Поступила в редакцию 16.01.15

STUDYING THE SPECIES EQUILIBRIUM DENSITY EQUATION IN SPACES OF DIFFERENT DIMENSIONS

Bodrow A. G., Nikitin A. A.

This article is devoted to the numerical solution of an integral equation appeared in a spatial biological model. Used numerical methods (Nystrom method and Neumann series method) lead to find an exact approximation in one-dimensional case. In this article we propose a way to reduce multi-dimensional case to one-dimensional and find out explicit formula for converting integral kernels in the case of normal kernels.

Keywords: integral equations, mathematical biology, numerical methods.