CC BY
34
О замыкании пространственных моментов в биологической модели, и интегральных уравнениях, к которым оно приводит
А.А. Никитин.
Аннотация—Настоящая статья посвящена математической постановке и численному исследованию пространственной модели стационарных биологических сообществ У. Дикмана и Р. Лоу. Главная идея данной модели состоит в том, чтобы найти «проекцию» симулируемого биологического процесса на некоторые характеристики, динамика которых может быть выписана аналитически. В качестве таких «характеристик» в модели Дикмана и Лоу выступают, так называемые «пространственные моменты» Выписывается система интегро-дифференциальных уравнений, описывающая пространственную динамику в этой модели. Ставится проблема замыкания третьего пространственного момента, избавляющая систему от бесконечного числа уравнений. Выписываются примеры замыканий третьего момента, два из которых, так называемые ассиметричное и симметричное замыкания второй степени, изучаются более подробно. Ставится задача с линейным интегральным уравнением равновесия, которая получается при использовании ассиметричного замыкания. Объясняются недостатки полученной модели, в связи с её биологической интерпретацией, а также с применением этого замыкания при изучении двухвидовой модели. Далее формулируется нелинейная проблема, возникающая при использовании симметричного параметрического замыкания второй степени, проводится её численное исследование. Показывается её правомерность, в том числе и в двухвидовой модели.
Ключевые слова—Математическая биология, individualbased model, интегральные уравнения, численные методы;
I. Введение
В последние годы компьютерное моделирование занимает все большее и большее место в различных научных областях человеческого знания. О его месте в современной науке можно судить уже по тому, что в 2013 году три химика М. Карплус (M. Karplus), М. Левитт (M. Levitt) и А. Уоршел (A. Warshel) получили нобелевскую премию по химии за «создание многоуровневых моделей сложных химических систем». В своей нобелевской лекции [1], Уоршел говорит о том, что математическое моделирование позволяет увидеть, как те или иные фи-
Статья получена 3 сентября 2018.
Работа выполнена при поддержке гранта Российского научного фонда, проект №17-11-01168.
А.А. Никитин, доцент, Московский государственный университет имени М.В. Ломоносова, факультет ВМиК, Москва, Ленинские горы, 119991; с.н.с., Российский университет дружбы народов. 117198, Москва, ул. Миклухо-Маклая, 6. (e-mail: nikitin@cs.msu.ru)
зические факторы влияют на итоговое поведение химической системы. Это позволяет глубже понять зависимости внутри системы и согласовать модель с реальными наблюдениями.
Находит компьютерное моделирование место и в биологии. Например, в 2007 году с помощью имитационной модели было показано, что, так называемы «материнский эффект» - негенетическая передача информации от матери к потомству - может быть одной из причин, обуславливающий популяционную динамику рачков Daph-nia Longispina [2]. Авторы обнаружили, что если убрать из модели «материнский эффект», то популяция неотвратимо вымирает. Существование самого материнского эффекта было экспериментально проверено ранее в работе [3].
В настоящей статье будет рассмотрена математическая модель стационарного биологического сообщества. Данная модель была предложена в работах [4]- [6], и уже изучалась нами в работах [7] - [11]. Её отличительной особенностью является присутствие пространственных неоднородностей. Ранее было обнаружено [5], что использование пространственных моделей, учитывающих поведение каждого индивида в отдельности, приводит к результатам более близким к реальности, чем результаты пространственно-однородных, «хорошо перемешанных» (well-mixed) моделей.
В первом разделе проводится краткое описание модели У. Дикмана и Р. Лоу, показано, что число уравнений динамики больше числа неизвестных функций, участвующих в них, и предложена идея замыкания пространственных моментов, т.е. попытка выразить некоторые неизвестные функции через остальные, исходя из биологического смысла. Во второй - ставится математическая проблема, и показывается, что она связана с вышеупомянутым замыканием пространственных моментов. В третьей - предлагается иное замыкание, и численно показывается правомерность этого действия. В заключении подводятся итоги, ставятся новые задачи и намечаются пути их решения.
II. Individual-based model У. Дикмана и Р. Лоу.
Для упрощения выкладок, рассмотрим одновидовую модель на прямой (n=1). Двухвидовая модель была изучена в работе [10], многомерные - в работах [8] и [11].
Основная идея работ [4]-[6] состоит в том, чтобы найти «проекцию» симулируемого биологического про-
цесса на некоторые характеристики, динамика которых может быть выписана аналитически. В качестве таких «характеристик» в вышеупомянутых работах выступают, так называемые «пространственные моменты» (spatial moments). Идея применения этих характеристик была предложена в работе [12], а в [6] проведена формализация и обобщение этой идеи.
Рассмотрим числовую прямую, в некоторых точках которой находятся биологические индивидуумы. Их расположение описывается функцией плотности р(х) (spatial pattern), представляющий собой сумму Ô-функций, имеющих особенности в точках расположения индивидуумов. Индивидуумы в изучаемой точке могут рождаться и умирать, причём вероятность смерти и рождения зависят от других индивидуумов, соседних с данной точкой. Перемещение индивидуумов происходит только при рождении.
Смерти индивидуумов бывают двух видов - от внешних причин и от конкуренции с другими индивидуумами. Внешние причины предполагаются однородными во всём пространстве, и плотность вероятности смерти от них равна d (d>0) для каждого из индивидуумов. Этот параметр называется «темпом смертности». Смертность от конкуренции описывается функцией d ы(х ~ х ), отвечающую за силу взаимодействия между индивидуумами в точках x и x'. Данную функцию будем называть «ядром взаимодействия», и предполагать её неотрицательной, нормированной и чётной. Параметр d' считаем постоянным.
При рождении индивид «выбрасывает» потомка в произвольную точку прямой, где он сразу же становится новым самостоятельным индивидуумом. Плотность вероятности выброса потомка на расстояние х' от родителя равна т(х'}. Здесь m - функция плотности некоторой случайной величины, называемая «ядром рассеивания» или «ядром рождаемости». Сама плотность вероятности рождения индивидуумом в точке x потомка не зависит от точки и равна b (b>0). Этот параметр называется «темпом рождаемости».
Все перечисленные величины: d, b, т(х) и со(х) мы считаем заданными биологической моделью. Перечислим накладываемые на них ограничения:
lim m(x) = lim ы(х) = 0,
1■ ' IVI
r +00 ç +00
I m(x)dx = J $dx = 1
J -co J —CO
Для введения в рассмотрение пространственных моментов, зафиксируем некоторый момент времени и рассмотрим, введённую выше, функцию р(х), равную рОО = Er»е.\- ~~ х') > где X - множество позиций всех индивидуумов. Смысл данной функции состоит в том,
что интеграл fA p(x)dx определяет число индивидуумов, находящихся в некоторой области A.
Далее введём среднюю плотность популяции (первый пространственный момент) N(p):
N(p) = f pCOdt,
л ■ (1)
а также среднюю плотность пар (второй пространственный момент):
С(х,р) = -Ь- Г pWOCt +ж) - m)dt.
А - (2) По сути, это нормированное на площадь, рассматриваемой области, число пар индивидуумов, расстояние между которыми равно а\ и в которых оба индивида различны.
Наконец, средняя плотность троек популяции (третий пространственный момент) Т(х,х ,р)-
Г<х,*',р)=щ| pC00Ct+*)-ffW)-
■ (p(t + х') - £(>') - 8(х - *')) dt
Биологический смысл плотности троек - нормированное на площадь области число троек индивидуумов, расстояние от первого до второго из которых равно х. а от первого до третьего -х'.
В настоящей работе вместо пространственных моментов N (р)i С (х, р) и Г (х, л', р) будем рассматривать их математические ожидания по паттерну р(х) и обозначать их соответственно С(х)и Т(х, х')
Было доказано, что процесс изменения пространственного паттерна РМ во времени описывается с помощью системы интегро-дифференциальных уравнений (master equations). Полный вывод этих уравнений производится в [6] и достаточно громоздок. Приведём уравнения динамики первых двух пространственных моментов. Интегро-дифференциальное уравнение, описывающее динамику средней плотности популяции N имеет вид:
d Г+со
— N=bN-dN-d'\ СОП^ОЛ^' at -у- (3)
Первое слагаемое представляет собой изменение плотности вида, вызванное рождением новых индивидуумов, второе - изменение вызванное смертью старых индивидуумов из-за внешних причин, третье - изменение плотности вида, вызванное смертями из-за конкуренции.
Уравнение динамики второго момента получается в виде:
jtC(0=b m{-ON + b J_+Jm(f')CCf +
f+ОЭ
-de (f) - d'wcocco - d' oiinm, nd?.
J—со
Первое слагаемое демонстрирует влияние на плотность пар
рождения нового индивидуума. Второе -влияние на рождения нового индивидуума родителем, находящимся на расстоянии ? от исследуемой точки (по всем таким расстояниям производится интегрирование). Третье слагаемое демонстрирует вероятность смерти одного из индивидуумов от внешних причин. Четвёртое - отражает вероятность того, что один из индивидуумов в паре умер из-за конкуренции с другим индивидуумом. А пятое - что один из индивидуумов в паре умер из-за конкуренции в тройке.
Заметим, что уравнение динамики первого момента использует второй момент, уравнение динамики второго момента - третий. Уравнение динамики третьего про-
странственного момента будет использовать четвертый момент и так далее. Поэтому количество неизвестных функций (пространственных моментов) всегда будет превышать число уравнений. Для того, чтобы избавиться от подобных зависимостей было предложено «замкнуть» уравнения динамики, исходя из биологических ограничений на исследуемые функции, [13]. Т.е. выразить некоторый пространственный момент через предыдущие.
Если мы выразим второй момент через первый СО) * N2,
то наша система интегро-дифференциальных уравнений преобразуется к хорошо изученному дифференциальному уравнению Ферхюльста [18]:
— /V = О - d)N - dN2, dt
которое, как известно, лишено пространственной структуры. Замыкание четвертого пространственного момента (через первый, второй и третий) приводит к весьма сложной системе интегро-дифференциальных уравнений. Поэтому, единственно разумным представляется замыкание третьего момента через первые два. Большой анализ и изучение данного вопроса было проведено в статье [13]. В этой работе на функцию замыкания Т » F(N,C)
были наложены три следующих условия (которые вытекают из биологического смысла), которые заведомо должны выполняться: lim Г (л:, х ) =
Если С(х) = N2, то Т(х, х') = N3.
В той же работе [13] был предложен ряд «кандидатов» на замыкание третьего момента:
2N
С(х)С(х')
(6)
Т(х.х')
X ■ (7)
С (ж) С (ж' )с(х-хг)
N3
(8)
ющеи системы уравнении: О = bN — dN — d'
cCfXiW
О = Ъ т{-ОК + Ъ + ГЖ'- (9)
Г+0О
-йсю - «¡'«те - «сггаг)<ц1.
J—со
Выражая первый пространственный момент N из первого уравнения, и подставляя его во второе, получаем следующее выражение:
=
-+СО ^ ^
b-d
Г Т IAJ
d'I toi fOrCf.fOdf.
J— со
которое называется «интегральным уравнением равновесия». В данное уравнение подставляются замыкания из предыдущего раздела. Ранее были широко изучены уравнения, получающиеся при замыканиях (6) и (7). В случае асимметричного замыкания второй степени, (6)
получается лине иное интегральное уравнение:
(Ь + = Ь I т(т)С(г + f)cÎT +
J — со
+
ъ
b-d
+ СО
(11)
ÎTl(^) ci'w(r)C ( r)ii r.
Это интегральное уравнение, которое и было изначально предложено математическому сообществу Ульфом Дикманом в 2005 году было достаточно подробно изучено, см, например, работы [14], [15], [16], [8]. Было доказано, что при некоторых дополнительных условиях, накладываемых на исследуемые функции, решение уравнения (11) существует только при условии d=0. Не существование решения данного уравнения при d>0 можно показать и графически. Ниже изображены решения видоизмененного уравнения (11) (несобственный интеграл по числовой прямой заменен на собственный по отрезку [-А,А])
г -t-A
(Ь + d'<u(0)c(0 = ь I ™(т)с(т+
J-A
(12)
+
b-d
тЮ
С+А
d'
J-А
ù)(r)C(j)dT.
для различных значений переменной интегрирования А.
Замыкания называются по степени второго момента C, которая в них используется. Т.о. замыкание (5) называется замыканием первой степени, а замыкание (8) -третьей. Замыкание второй степени, (6) называется симметричным, а (7) - асимметричным.
Ранее было показано [13], что наиболее близкими к данным компьютерных симуляций являются замыкания второй степени (6) и (7).
III. Уравнение равновесия
В предыдущем разделе была рассмотрена система уравнений динамики первого и второго пространственных моментов. В нашей работе мы рассматриваем данную систему в стационарной точке, т.е. решение следу-
1.045 1.04 I.0Î5 1.03 1.025 1.02 1.015 I 01 1.005 1
Рис. 1 График зависимости второго пространственного
момента, С для различных отрезков интегрирования [-А,А].
Как мы видим, решение интегрального уравнения равновесия сильно зависит от длины отрезка интегрирования, что демонстрирует отсутствие решения этого уравнения. Аналитическое доказательство отсутствия решения при ^0 проводится, используя прямое и обратное преобразование Фурье. Подробности см., например, в работе [8].
Данный факт (отсутствие решения при ^0) существенно снижает ценность решаемой биологической задачи, т.к. случай d=0 соответствует модели, в которой отсутствует «естественная смертность», и индивидуумы могут умирать исключительно от конкуренции друг с другом.
Но, между тем, изучение случая асимметричного замыкания второй степени (7) и интегрального уравнения, к которому оно приводило, было продолжено. Так, были установлены достаточные условия существования решения , [14]-[16]. А также исследованы многомерные случаи в пространствах одной, двух и трёх размерностей (п=1,2,3). Для интегральных ядер т и ю взятых в виде гауссиан:
^0") = —-5—777, га(г) = —-——-
был разработан, так называемый «метод понижения размерности» и получена зависимость средней плотности популяции (первого пространственного момента) от значений интегральных ядер (точнее от значений параметров ат и аш). Ниже приведены графики этих зависимостей в одно-, двух- и трёх- мерных случаях, [8].
Подробное описание, и рассуждения о приведённых графиках см в работе [7]. Изучение аналогичной задачи для ядер-куртозиан, [17] в пространствах различных размерностей проводилось в работе [9].
ю
Рис. 2 График зависимости средней плотности популяции в зависимости от ядер конкуренции и движения в одномерной среде.
2D
Рис.3 График зависимости средней плотности популяции в зависимости от ядер конкуренции и движения в двумерной среде.
3D
Рис. 4 График зависимости средней плотности популяции в зависимости от ядер конкуренции и движения в трёхмерной среде.
Окончательный отказ от асимметричного замыкания второй степени произошёл после работы [10], в которой изучалась двухвидовая модель У. Дикмана и Р. Лоу [4], [5]. Задача, заключающаяся в нахождении первого и второго пространственных моментов для модели двух видов индивидуумов (даже в случае отсутствия смертности от окружающей среды, d=0) приводит к системе из пяти интегральных уравнений, два из которых аналогичны уравнению равновесия (9), в которых коэффициент, стоящий у интеграла в правой части равенства не обращается в единицу, что приводит к отсутствию решения этого уравнения, см. [15], [8].
Исходя из полученного результата, был сделан вывод об отсутствии равновесия в двухвидовой модели при использовании асимметричного замыкания моментов второй степени. Данный факт привёл к необходимости использования другого замыкания третьего момента, которое мы будем исследовать в следующем разделе.
IV. Параметрическое замыкание второй степени.
Проблемы, изложенные выше, привели к необходимости использования другого вида замыкания третьего момента. Нами было использовано параметрическое замыкание, являющееся линейной комбинацией замыканий (6) и (7), которое было предложено и исследовано в работе [13]:
ТП п -■ а, ссрссг - О , сстг - О Л
2\ N N N I
гауссианами, d>0.
+(1 - а)
С(_х)С(х')
. (13)
Данное замыкание при подстановке в систему (9) приводит к системе с нелинейным интегральным уравнением:
d'î^an^'W
N =
b-d
f + œ
(d - d'w(f))C(0 - b m(f)W - b C(f + f) m(f') dÇ' =
J — 00
f+œ
c(0 d'ccf)« (fVf-
J—oo
(f +CQ г +00 \
C(f) d'ccfMOdf + cif) d'C(f+f)ö(Odf -
J — cл J—œ /
1 — a N
-2 -1,5 -1 -0,5
2N
-J^A \ d'otf + - d'N4
Аналитическое изучение второго уравнения пока не представляется возможным. Поэтому, мы подошли к его решению численно, используя метод рядов Неймана [14]. Подробности реализации данного метода могут быть найдены в работе [11].
Ниже приведены графики равновесного второго момента для различных данных модели. Видно, что решение быстро выходит на требуемую асимптоту при аргументе, стремящемся к бесконечности.
1.25
Рис. 6 График равновесной парной плотности плотности с ядрами-гауссианами, d>0.
Хорошо себя зарекомендовало это замыкание и в двухвидовой модели. Задача по её исследованию была сформулирована в виде системы из пяти интегральных уравнений. Предложен алгоритм по решению данной системы, который был применен к входным данным, вытекающим из классических биологических сценариев: Competition-colonization trade-off и Heteromyopia. Подробности могут быть найдены в работе [10].
В завершении мы приведём график, в логарифмической шкале, изменения разности двух соседних (по количеству итераций) вторых моментов. Как видно, убывание к нулю практически линейное, а значит, мы имеем экспоненциальную сходимость.
Рис. 5 График равновесной парной плотности с ядрами-
Рис. 7 График разности старого и нового интеграла С:-Г Ь.;х в зависимости от числа итераций.
Цветом обозначено значение параметра d. Чем больше значение d, тем «краснее» график.
V. Заключение
В настоящей работе была рассмотрена модель, основанная на индивидуальностях (individual-based model) У. Дикмана и Р. Лоу. Была поставлена проблема замыкания пространственных моментов и описаны два примера таких замыканий. Далее было объяснено, почему асимметричное замыкание второй степени (7), которое изначально казалось биологам наиболее удачным и приводило к линейной математической задаче, оказалось несостоятельным. Предложено иное параметрическое замыкание третьего момента (13), в исследовании которого и видятся пути исследований. В дальнейшем требуется провести большую работу над созданием адекватных компьютерных симуляций, которые были бы стабильны на большом множестве параметров модели. После создания симуляций, которые бы выходили на стабильное состояние равновесия, требуется провести работу над математической постановкой. Согласно результатам проверки на корректность различных значений параметра а в (13), наиболее близкие к симуляциям результаты достигаются при а=2/5, см. [13]. Однако, стоит отметить, что используемое замыкание необязательно является наиболее подходящим и задача о поиске наиболее оптимального замыкания является открытой и в известной мере одной из самых актуальных.
Вполне может быть, что адекватные компьютерные симуляции могут поднять вопрос о корректности системы интегро-дифференциальных уравнений, описывающих динамику пространственных моментов (3), (4), или привести к добавлению (или исключению) некоторых слагаемых в них.
Автор благодарит Ульфа Дикмана за постановку задачи, и внимание к нашей работе. Также он благодарит своих коллег-учеников А. Бодрова и А. Савостьянова за полезные обсуждения и работу над приведёнными в настоящей статье иллюстрациями.
Библиография
[1] Warshel A. Multiscale modeling of biological functions: from enzymes to molecular machines (Nobel Lecture). // Angew Chem. Int. Ed. Engl. - 2014. - Sept -Vol. 53, no. 38.
[2] Алексеев В.Р., Казанцева Т.И. Использование индивидуально-ориентированной модели для изучения роли материнского эффекта в смене типов размножения у Cladocera // Журнал общей биологии. -2007. - т.68, №3. - 231-240. (in russian)
[3] Alexeev V., Lampert W. Maternal control of resting-egg production in Daphnia // Nature. - 2001. - Vol. 414 -899-901.
[4] Dieckmann U., Law R. Introduction // The Geometry of Ecological Interactions: Simplifying Spatial Complexity / ed. by U. Dieckmann, R. Law, J. Metz. - Cambridge University Press, 2000. - 1-6.
[5] Dieckmann U., Law R. Relaxation projections and the method of moments // The Geometry of Ecological Interactions: Simplifying Spatial Complexity / ed. by U. Dieckmann, R. Law, J. Metz. - Cambridge University Press, 2000. - 252-270.
[6] Dieckmann U., Law R. Relaxation projections and the method of moments // The Geometry of Ecological Interactions: Simplifying Spatial Complexity / ed. by U. Dieckmann, R. Law, J. Metz. - Cambridge University Press, 2000. - 412-455.
[7] Bodrov A. G., Nikitin A. A. Qualitative and numerical analysis of an integral equation arising in a model of stationary communities // Doklady Mathematics. — 2014. — Vol. 89, no. 2. — P. 210-213.
[8] Bodrov A. G., Nikitin A. A. Examining the biological species steady-state density equation in spaces with different dimensions // Moscow University Computational Mathematics and Cybernetics. — 2015. — Vol. 39, no. 4. — P. 157-162.
[9] Kalistratova A. V., Nikitin A. A. Study of dieckmann's equation with integral kernels having variable kurtosis coefficient // Doklady Mathematics. — 2016. — Vol. 94, no. 2. — P. 574-577.
[10]Nikitin A. A., Savostianov A. S. Nontrivial stationary points of two-species self-structuring communities // Moscow University Computational Mathematics and Cybernetics. — 2017. — Vol. 41, no. 3. — P. 122129.
[11] Nikitin A.A., Nikolaev M.V. Equilibrium Integral Equations with Kurtosian Kernels in Spaces of Various Dimensions // Moscow University Computational Mathematics and Cybernetics. — 2018. — Vol. 42, no. 3. — P. 105-113.
[12] Bolker B., Pacala S. Using moment equations to understand stochastically driven spatial pattern formation in ecological systems // Theoretical Population Biology. -1997. - Vol. 52.
[13] Murrel D.J., Dieckmann U., Law R. On moment closures for population dynamics in continuous space // Journal of Theoretical Biology. - 2004. - Vol. 229. -421-432
[14] Данченко В.И., Давыдов А.А., Никитин А.А. Об интегральном уравнении для стационарных распределений биологических сообществ // Проблемы динамического управления. - 2010. - №3. - 15-29. (in russian)
[15] A.A. Davydov, V.I. Danchenko, M.Yu. Zvyagin, 2009, published in Trudy Matematicheskogo Instituta imeni V.A. Steklova, 2009, Vol. 267, pp. 46-55.
[16] Данченко В.И., Рубай Р.В. Об одном интегральном уравнении стационарного распределения биологических систем. // Современная математика. Фундаментальные направления. Том 36 (2010). 50-60. (in russian)
[17] Seier E., Bonett D. Two families of kurtosis measures // Metrika. 2003. 58. 59-70.
[18] Verhulst P.F. Notice sur la loi que la population poursuit dans son accroissement // Correspondance mathematique et physique. - 1838. - Vol. 10. - 113-121.
Никитин Алексей Антонович родился в г. Москва 14.02.1983, В 2000 году закончил Лицей информационных технологий №»1533, в 2005 году закончил факультет ВМиК МГУ им. М.В. Ломоносова, в 2008 году защитил диссертацию на соискание степени к.ф.-м.н. на тему «Граничное управление третьим краевым условием» под руководством В.А. Ильина. С 1 мая 2008 года работает на кафедре Общей математики факультета ВМиК МГУ в должности ассистента, и доцента с октября 2013 года.
Область научных интересов: математическое моделирование, нели- модернизации и информатизации системы образования, использова-нейные интегральные уравнения, математическая биология, вопросы ние визуальных образов в процессе образования.
On the closure of spatial moments in the biological model, and the integral equations to which
it leads
A.A. Nikitin
Abstract - This article is devoted to the mathematical formulation and numerical study of the spatial model of stationary biological communities by U. Dickman and R. Low. The main idea of this model is to find a "projection" of the simulated biological process on certain characteristics, the dynamics of which can be written analytically. Such "characteristics" in the Dickman' and Low's model are the so-called "spatial moments". A system of integro-differential equations is described describing the spatial dynamics in this model. The problem of moment's closure of the third spatial moment, which relieves the system of an infinite number of equations is posed. Examples of closures of the third moment are written out, two of which, the so-called asymmetric and symmetric closures of the second degree, are studied in more detail. The problem is posed with a linear integral equilibrium's equation, which is obtained by using an asymmetric closure. The shortcomings of the obtained model are explained, in connection with its biological interpretation, and also with the application of this closure in the study of the two-species model. Next, we formulate a nonlinear problem that arises when a symmetric parametric closure of the second degree is used, and its numerical study is carried out. Its legitimacy is shown, including in the two-species model.
Keywords - Mathematical biology, individual-based model, integral equations, numerical methods;
