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

УДК 519.642.4

А. А. Никитин, А. С. Савостьянов2

НЕТРИВИАЛЬНЫЕ СТАЦИОНАРНЫЕ ТОЧКИ ДВУХВИДОВЫХ САМОСТРУКТУРИРУЮЩИХСЯ СООБЩЕСТВ*

В статье рассматривается двухвидовая модель самоструктурирующихся стационарных биологических сообществ, предложенная У. Дикманом и Р. Лоу. Разработан метод изучения системы интегро-дифференциальных уравнений, описывающей положение равновесия в модели; найдены нетривиальные стационарные точки; исследованы ограничения на пространство параметров модели, приводящие к подобным стационарным точкам. Полученный метод применен к ряду общеизвестных биологических сценариев.

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

1. Введение. В работе рассматривается стационарное положение биологических сообществ двухвидовой популяции, конкурирующей за ресурсы (одновидовые работы изучались ранее в [1, 2]). Данный вопрос является одним из самых естественных и в то же время содержательных вопросов в данной области. Ввиду этого, описанные ранее модели прошли существенную эволюцию, рассмотренную в [3], от аналитически простых, исходящих из пространственной однородности биологического сообщества, до моделей пространственно неоднородных популяций (Individual Based Models, IBM), требующих большой вычислительной мощности. В данной работе рассмотрена модель, предлагающая сравнительно неосложненный аналитический аппарат для изучения пространственных неоднородностей, возникающих в процессе самоструктуризации популяции; изучены стационарные положения системы на предлагаемом моделью пространстве параметров, а также влияние размерности области обитания сообщества на найденные стационарные положения.

2. Пространственная модель стационарных сообществ. Рассматриваемая в работе модель была предложена У. Дикманом и Р. Лоу в [4, 5].

2.1. Метод пространственных моментов. Для характеристики состояния сообщества необходимо рассматривать не только количественные характеристики каждого вида, но и их пространственную структуру, а также события, происходящие с каждым индивидом из рассматриваемой популяции. Пусть сообщество находится в ограниченной области А объема

Конкретное распределение особей вида г в области А будем называть паттерном г-го вида

где 6(х) — дельта-функция Дирака, а X С А — множество точек, в которых присутствует особь г-го вида. Совокупность паттернов всех видов полностью описывает состояние сообщества в конкретном распределении особей.

Для каждого паттерна введем следующий набор статистик:

• Ni(p) = - Pi(x) dx — плотность индивидов в области при данном распределении;

• = т-77 / pi(x)pj(x + £) dx — плотность пар индивидов видов {i — j), где индивид вида

J

А

j удален от индивида вида г на расстояние

1 Факультет ВМК МГУ, доц.; РУДН, ст. науч. сотр., к.ф.-м.н., e-mail: nikitinQcs.msu.su

2 Факультет компьютерных наук НИУ ВШЭ, студ., e-mail: a.s.savostyanov@gmail.com

* Результаты разделов 2, 4, 6 получены первым автором при поддержке Российского научного фонда, проект

х'£Х

А

1

№ 17-11-01168.

• Tyfc(£', = т-jr f pi(x)pj(x + С)Рк{% + С') dx — плотность троек индивидов видов {i — j — к),

\А\ J

А

где индивиды видов j и к удалены от индивида вида г на расстояния £ и соответственно. Аналогичным образом можно ввести плотности пространственных структур и более высокого порядка. В данной работе рассматривается двухвидовая модель, т. е. i,j, к = 1,2.

Пространственными моментами первого, второго и третьего порядков будем называть математические ожидания по пространству паттернов введенных выше статистик

Ni(t)=EpNi(p), Cij(t,t)=EpCij(t,P), Tijk(W,t)=¥vTijk(W,p),

где вероятность реализации каждого паттерна зависит от времени t.

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

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

• Вероятность рождения потомка вида г в точке 'J от родителя, находящегося в точке обозначим Bi(= т— £), где функция nii(x) — ядро рождения (dispersal kernel),

J nii{x) dx = bi — темп рождаемости, 0 < bi < 1.

Rn

• Вероятность смерти конкретного индивида г-го вида, находящегося в точке обозначим

А(с,р) = di + f wa (с - о мо -

j А

где di — вероятность смерти от влияния окружающей среды, 0 < di < 1; ядро конкуренции Wij(x) — плотность вероятности смерти индивида г-го вида от конкуренции с индивидом j-го

вида, находящимся на расстоянии х, J Wij(x) dx = d'i;j — сила конкуренции, 0 < d\j < 1.

Rn

В рамках данной работы будем считать, что величины rrii(x)/bi и Wij(x)/d'ij распределены нормально с нулевым математическим ожиданием (nii ~ N(Q,afl) Wij ~ iV(0, ст^)), причем обе функции являются радиально-симметричными, т.е. тДж) = тД|ж|) и Wij(x) = Wij(|ж|), исходя из биологических соображений.

2.3. Динамика пространственных моментов. При помощи введенных выше величин можно описать динамику предложенного набора статистик, формальный вывод которой проведен в [4]. Приведем результирующие уравнения динамики для моментов первого и второго порядка.

Для моментов первого порядка имеем

Jt Щ = (bi - di)Ni - ¿; J WijiOCijiO dt (1)

3 = %п

Несложно заметить, что все слагаемые в правой части имеют простую биологическую интерпретацию, которая приведена в статье [4]. Аналогично может быть выписана динамика второго момента:

d_ dt

CiM) = Sijnii^ONi + Jmi{OCiM + C') di,' - diCijiO-

V / wik(aTijk(t,£)d^ WijiOCijiO + (2)

k=i

Здесь также можно охарактеризовать каждое слагаемое. Член (г,^, £ —> г, —£) означает, что далее следует повторить рассуждения, но теперь уже для "второй" особи из пары с поправкой на изменение ориентации расстояния.

2.4. Замыкания пространственных моментов. Как видно выше, динамика первого момента Щ зависит от второго момента Су (£). Таким же образом динамика второго момента Су (£) зависит от третьего момента Tyfc(£,£'). Подобная зависимость объясняется тем фактом, что отрицательная динамика пространственных моментов зависит от конкурентной гибели участников пространственной структуры, т. е. зависит от плотности пространственных структур следующего порядка.

Для разрешения возникающей иерархии зависимостей используются замыкания пространственных моментов, т. е. аппроксимация пространственных структур более высокого порядка через пространственные структуры низших порядков. В частности замыкание Су (£) = N^Nj приводит к известной динамике Лотки-Вольтерры. В настоящей работе будут рассматриваться замыкания третьих моментов через моменты первого и второго порядка.

Заметим, что если пространственная структура стала однородной на низшем уровне, то она должна остаться такой и на более высоких, — это означает, что необходимо наложить на рассматриваемые замыкания некоторые требования:

1) если Су = NiNj, то Tijk = NiN:jNk;

2) lim Су (С) = ^Nj — отсутствие пространственной структуры у бесконечно удаленных точек;

3) lim £') = NiCjkÜ' — С) — бесконечно удаленная точка уничтожает пространственную

структуру, приводя к пространственному моменту низшего порядка;

4) аналогично введенному выше lim Tyfc(£, £') = Су(£)Л^.

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

т ,, « fCij(0Cik(?) , CijiQCj^'-Q , Cik(?)Cjk(? - Q \ 4Cy(g)Cifc(f)

Tiikiu } = 2 V-Ж-+-Щ-+-Щ--V ( }-Ni-'

(3)

В работах [6, 7] было выбрано а = 4/5 как наиболее подходящее под проведенные биологические симуляции.

3. Стационарные точки системы. Проведем замыкание (3) в выписанных ранее уравнениях динамики (1) и (2). Далее будем искать стационарные точки системы:

= <цу2(*) = <1Си(Ы) = = ^ йС22{^1) = ^ (4)

ей ' ей ' ей ' м ' м '

где = 1,2. В результате получим систему из 5 интегральных уравнений для неизвестных N1, N2, Си(С), С12(С) и СггСС)- Прежде чем построить эту систему, приведем несколько обозначений и соглашений.

• Будем опускать зависимость пространственных моментов от времени, поскольку ищется стационарная точка.

• Второй и третий моменты нормируем и центрируем на их пределы на бесконечности:

= тМ = АЖ) = СУ(0-1.

Таким образом, будем иметь Нт -Оу(£) = 0.

• Введем общий суммарный уровень конкуренции между видами % и ] с учетом их пространственной структуры:

Уц = j WijiOCijiO

Обозначим через [/ * д] = / + £') с^' свертку двух функций / и д.

3.1. Результирующая система. Сворачивая уравнения (4), получим систему = [т 1*С12Ш) + [т2 *С21]Ы) - [<11 + из12(0}С12(0 - |»и«')^!«, С)

-1 w22(C)т122((,C)d(l^[d2 + w21Ы)}C2lЫ)^ I у>21 (Оад-^'К-

— <¿22 I

В нормированном случае она будет иметь вид

= [(тщ + т2) * С12] (С) - [йг + ё2 + у)12(0 + го21(£)]С12(£) - ^ |»21(ОТ2ц(-{, С')

— -^2 I у>22(£')Т212Ы,?)<%'-М1 I ыиI «;12(т22(Ш<е

Выпишем нормированное замыкание

ГУ Г 1

Тг^Ы') = 2 [^'(С)^(С) + - С') + Сгк(?)С,к(£ - С') - 1] + (1 - а)Су(£)Сг*(£')

и запишем систему уравнений для первых моментов

Ьг-(11= N^11 + Ж2У12, Ь2 - й2 = ЛГ1У21 + ^22/22,

^1 + N2 I ы22(ОТ212Ы^')^' + ^ I >т11Ц')Т121{^0<Ц'+

+ Ж2 У и112((,')Т122((,,(,') = N1С12У21 + ^2С\2у22 + Л^С^уп + М2С12у12 =

= (Ъг + Ъ2 - (¿1 - й2)С12($) = (¿1 + Ь2 - 1 - Й2)£>12($) + (¿1 + Ьг ^ 1 - й2)(!-!)>

где через обозначена подстановка £') = Су (^)С^(^'). В итоге получим

О! \ О! 1

1 ~ 2 / ^ + ^ + + ^ + 1012 + 1)12 = + * ~ 1012 ~

а Г , ,

- - ]В12 (Ж1 [ади *В21] + Ж2 [ила *£>22] + Ж1 [ад21 *£>ц] + Н2 [ш22 *1>12] + Н2й'12 + Н1 ё21+Ы2ё22 +Ы^) +2Л;,[7Г,, + + */.> 11 ] + 2Л7^[<гТ2*I) 12] + /V1 1)и\ +Ъ-^г-,-,I)-,-,* 1)^]

+ ЛГ, [/Г], />! , * В12] + * В22]

Преобразуя это соотношение, найдем

]в12 =

= [(гп1 + т2) * 1)\-2\ - ад 12 -(£>12 + 2) ([гои * />!•_.] + [го21 * £>п]) + [го21£>12 * Вп] + [®цОц * -

(£>12 + 2) ([и; 12 * -С22] + [го22 * В12}) + [ш22В22 * В12] + [п)12В12 * В22] . (5) Аналогично для двух оставшихся вторых моментов будем иметь

(X \ (X

1 - 2 ) (к + ь2) + - ((¿1 + й2 + й'11Н1 + й'12Н2 + й'21Н1 + й'22Н2) + «712 +

2^

(l - |)bi + |(di + Nld>n + N2d!l2) + wn

(Du + 2)[wn * Dn] + [wnDn * Du]

Dи = — + [mi * Du] - wir iVi

« i „ -<! JVi

No

(Dn + 2) [«712 * D12] + [wi2Di2 * D12]

1 - 77 + ^((¿2 + JVi4i + N2d'22) + w22 D22 = ^ + [m2 * D22] - W22~ Z / Z J iV2

а

N,

(D22 + 2)[ад22 * D22] + [«722D22 * D22]

N1

(D22 + 2) [«721 * -C12] + N21D12 * DI2]

Решая систему уравнений для первых моментов методом Крамера, получим

JVi =

(h - di)y22 - (b2 - (¿2)1/12

N0 =

(b2 - d2)yu - (h - (¿1)2/21

У11У22 - У12У21 У11У22 - У12У21

Таким образом, система уравнений (5)-(8) описывает искомую стационарную точку.

(6)

(7)

(8)

4. Численное решение системы. Полученная система уравнений (5) (8) представляет собой систему нелинейных интегральных уравнений и может быть записана в виде

N1 = ^(уц,у12,у21,у22), N2 = Ь2(уц,у12,у21,у22),

Сц = Кп(Мъ N2,Cn, С12), С12 = Ki2(Ni, N2,Cn, С12, С22), С22 = K22(Ni, N2, С21, C22),

где Li линейные интегральные операторы, вклк,,ак,щие /„!уК)Сц (0 а К„ интеграль-

Кп

ные операторы с нелинейностями вида Сц~((,) • [и)ц~ * Сц~}((,) и [(тц. • Сц~) * Сц~}((,)-

4.1. Метод рядов Неймана. Для решения построенной системы используем метод рядов Неймана. Пусть дано уравнение вида

/=д + к1 (9)

где К — некоторый оператор, д — известная, а / — искомая вектор-функции.

Теорема 1. Если ряд д + Кд + К2д + ... + Кпд + ... сходится, то он сходится к решению уравнения (9).

Теорема 2. Если оператор К — сжимающий, т. е. выполнены условия

За, 0 < а < 1: Ух, у => р(Кх, К у) < а • р(х, у),

то ряд д + Кд + К2д + ... + Кпд + ... сходится.

Для доказательства существования и единственности решения системы (5)-(8) можно использовать технику, разработанную в [8, 9], и сформулированные теоремы 1, 2. Вопросы, связанные с исследованием сходимости использованных численных методов, в работе не рассматриваются.

4.2. Ускорение вычислений в двумерном случае. Обратимся теперь к изучению двумерного случая. Наиболее сложным с вычислительной точки зрения вопросом является работа с двумерной сверткой. Для нее верно свойство преобразования Фурье:

[/«.(/I , ПП.П-Пи\\-

После дискретизации на квадратной сетке размера К х К прямое и обратное преобразования выглядят следующим образом:

К-1К-1 К-1К-1

п — _ V"' V"' /; рц^(тпи+пь)

"ни — / у / J -0тпс ; ■ь-пт — / J / J

п=1 т=1 п=1 т=1

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

1 К—1 г 1 к-1

I V—\ I V—\ '¿гтгъ

Guv — к к ж™»е к

К ^ К

п=1 L m= 1

_ 2iirто

e K

что дает алгоритмическую сложность 0(K3log(K)), поскольку необходимо вычислить быстрое преобразование Фурье (К1о%(К)) в каждом узле сетки. Данную асимптотику можно улучшить, применив преобразование Ханкеля [10].

Перейдем в преобразовании Фурье к полярным координатам

+ ОС +ОС +ОС 7Г

Р[/\(шх,шу)= ! I ¡(х,у)е~^х+шУу) = ^ I f(r,в)e-irpcos^-в)rdrdв.

—оо —сю 0 — 7Г

Используя радиальную симметричность исследуемых функций, получаем соотношение

+ СЮ 7Г +СЮ

¿?[/](р, Ф)= I г ¡(г) йг I е-ггрсоз^-в) м = 2ж I г/(г) 70(гр) dr,

0 —7Г 0

которое известно как преобразование Ханкеля 0-го порядка, где

ж

J0(x) = -L J e~ixcoSTdT

— функция Бесселя нулевого порядка. Если дополнительно сделать экспоненциальную замену переменных г = г0еж, р = Pqcv , то вместо преобразования Фурье получим свертку двух функций, алгоритмическая сложность которой 0(Klog(K)), т.е. было произведено ускорение в К2 раз.

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

Рассматривались два механизма. Первый механизм — это heteromyopia, в котором особи внутри вида конкурируют на большем расстоянии, чем особи разных видов (в рассматриваемой модели механизм реализован за счет изменения дисперсии ядер интравидовой aft = о22 и интервидой конкуренции af2 = (721). Параметры модели были выбраны так: ст^ = о22 и 0Т2 = a2i приняты равными;

Ьг = Ь2 = 0.4, dx = d2 = 0.2, d'n = d22 = d21 = d'l2 = 0.001, ст™ = ст™ = 0.06.

Зависимость первого пространственного момента от параметров первой модели приведена для одномерной (А С К1) (см. рис. 1, а) и для двумерной (А С К2) (см. рис. 1, б) области обитания.

Второй механизм — это competition-colonization trade-off, в котором один вид сильнее конкурирует, а другой имеет большую дисперсию распространения (в рассматриваемой модели механизм реализован за счет изменения дисперсии ядра рождения одного вида ог2 и силы конкуренции другого d'12). Оставшиеся параметры модели были фиксированы и выбраны следующим образом:

bi = b2 = 0.4, d1 = d2 = 0.2, d'n = d22 = d21 = 0.001,

„ТО n (\Л „W „W „W „W n ПЛ

ol = 0.04, au = a12 = a21 = a22 = 0.04.

Зависимость первого пространственного момента от параметров второй модели приведена для одномерной (А С К1) (см. рис. 2, а) и для двумерной (А С К2) (см. рис. 2, б) области обитания.

Вычисления проводились на суперкомпьютере "Ломоносов" в рамках проекта "Моделирование биоценозов при помощи суперкомпьютера".

Отметим, что полученные результаты в одномерном случае не противоречат результатам, полученным в исходной работе [11], за исключением переходных областей, что, видимо, вызвано лучшей точностью нашего численного метода. Результаты в двумерном случае были нам ранее

Рис. 2

неизвестны. Возможно, отсутствие данного случая в работе [11] было связано с тем, что сложность метода Рунге Кутты, используемого ими, возрастает квадратично при переходе в Ж2 .

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

Работа выполнена в ходе проведения исследования (16 05 0069) в рамках Программы «Научный фонд Национального исследовательского университета «Высшая школа экономики» (НИУ ВШЭ)» в 2016 2017 гг. и с использованием средств субсидии на государственную поддержку ведущих университетов Российской Федерации в целях повышения их конкурентоспособности среди ведущих мировых научно-образовательных центров, выделенной НИУ ВШЭ.

Авторы благодарят У. Дикмана за постановку задачи и внимание к работе, а также А. Бодрова и С. Гаджиева за помощь в подготовке рукописи и ценные обсуждения.

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

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

2. Бодров А. Г., Никитин А. А. Исследование интегрального уравнения плотности биологического вида в пространствах различных размерностей // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2015. № 4. С. 7-13.

3. Plank М. J., Law R. Spatial point processes and moment dynamics in the life sciences: a parsimonious derivation and some extensions // Bull. Math. Biol. 2015. 77. P. 586-613.

4. Dieckmann U., Law R. Relaxation projections and the method of moments // The Geometry of Ecological Interactions: Simplifying Spatial Complexity. Cambridge University Press, 2000. P. 252-270.

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

6. Law R., MurrellD. J., Dieckmann U. Population growth in space and time: spatial logistic equations // Ecology. 2003. 84. N 1. P. 252-262.

7. Murrell D. J., Dieckmann U. On moment closures for population dynamics in continuous space // J. Theor. Biology. 2004. 229. P. 421-432.

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

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

10. В add о иг N. Operational and convolution properties of two-dimensional Fourier transforms in polar coordinates //J. Opt. Soc. Amer. 2009. 26. P. 1767-1777.

11. Murrell J., LawR. Heteromyopia and the spatial coexistence of similar competitors // Ecology Letters. 2003. 6. P. 48-59.

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

УДК 519.6

В. И. Пагурова1

О ПРЕДЕЛЬНОМ МНОГОМЕРНОМ РАСПРЕДЕЛЕНИИ ПРОМЕЖУТОЧНЫХ ПОРЯДКОВЫХ СТАТИСТИК

Исследуется совместное асимптотическое при п оо распределение промежуточных порядковых статистик, построенных по выборке объема п.

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

Пусть Xi,..., Хп — независимые и одинаково абсолютно непрерывно распределенные величины, копии случайной величины X, F(x) = ~Р{Х ^ х}, f(x) = F'(x), Х^ — порядковая статистика в вариационном ряду X^ ^ ... ^ построенном по величинам Х\,..., Хп. Обозначим щ = [tina], г = 1,..., т, 0 < t\ < ... < tm, где ti не зависит от п, 0 < а < 1, [ж] означает целую часть числа х. Для различных классов распределений F исследуется асимптотическое при п ^ оо совместное распределение т промежуточных порядковых статистик ■ ■ ■

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

1 Факультет ВМК МГУ, ст. науч. сотр., к.ф.-м.н., e-mail: pagurovaQyandex.ru