CC BY
5
Известия высших учебных заведений. Прикладная нелинейная динамика. 2021. Т. 29, № 5 Izvestiya Vysshikh Uchebnykh Zavedeniy. Applied Nonlinear Dynamics. 2021;29(5)
Научная статья УДК 575.174, 517.925
DOI: 10.18500/0869-6632-2021-29-5-706-726
О генетической дивергенции в системе двух смежных популяций, обитающих на однородном ареале
Е. Я. Фрисман, М.П. Кулакови
Институт комплексного анализа региональных проблем ДВО РАН, Биробиджан, Россия E-mail: [email protected], [email protected] Поступила в редакцию 08.04.2021, принята к публикации 27.04.2021, опубликована 30.09.2021
Аннотация. Цель работы - исследование механизмов, приводящих к возникновению генетической дивергенции -устойчивых генетических различий между двумя смежными популяциями, связанными миграцией особей. Рассматривается ситуация, когда приспособленность особей жестко определяется генетически единственным диаллельным локусом с аллелями А и а, популяция панмиктична с менделевскими правилами наследования. Динамическая модель содержит три фазовые переменные: концентрации аллеля А в каждой популяции, а также доля (вес) одной из популяций в общей численности. При этом предполагается, что численности изменяются либо независимо со скоростями, определяемыми средними значениями коэффициентов отбора (мальтузианскими параметрами), которые зависят от концентрации аллелей и приспособленностей гомо- и гетерозигот каждой из популяций, либо строго синхронно. Методы. Для исследования модели использовались качественные методы исследования обыкновенных дифференциальных уравнений, включающие построение параметрических и фазовых портретов, бассейнов притяжения и бифуркационных диаграмм. Исследуются бифуркации, обеспечивающие принципиальную возможность генетической дивергенции. Результаты. Если у гетерозигот приспособленность выше, чем у гомозигот, то обе популяции оказываются полиморфными с одинаковой концентрацией гомологичных аллелей. В случае пониженной приспособленности гетерозигот и независимо от изменения численности популяций со временем в популяциях установится одинаковый мономорфизм по одному из аллелей. Динамика при этом оказывается бистабильной. Показано, что дивергенция в такой системе - результат субкритической бифуркации вил неустойчивого полиморфного состояния. В этом случае дивергентное состояние неустойчиво и проявляется в динамике как часть переходного процесса при движении к одному из мономорфных состояний. Заключение. Устойчивой генетическая дивергенция оказывается лишь для популяций, сохраняющих определенным образом соотношение численностей. В этом случае дивергенции предшествует седлоуз-ловая бифуркация, а динамика оказывается квадростабильной - в зависимости от начальных условий возможен либо мономорфизм, либо дивергенция.
Ключевые слова: генетическая дивергенция, дифференциальные уравнения, динамика, седлоузловая бифуркация, би-стабильность, квадростабильность.
Благодарности. Работа выполнена в рамках государственного задания Института комплексного анализа региональных проблем ДВО РАН.
Для цитирования: Фрисман Е.Я., Кулаков М.П. О генетической дивергенции в системе двух смежных популяций, обитающих на однородном ареале // Известия вузов. ПНД. 2021. T. 29, № 5. С. 706-726. DOI: 10.18500/0869-6632-2021-29-5-706-726
Статья опубликована на условиях Creative Commons Attribution License (CC-BY 4.0).
706
© Фрисман Е. Я., Кулаков М.П., 2021
Article
DOI: 10.18500/0869-6632-2021-29-5-706-726
On the genetic divergence of two adjacent populations living in a homogeneous habitat
E. Ya. Frisman, M. P. Kulakovи
Institute for Complex Analysis of Regional Problems, Far Eastern Branch, Birobidzhan, Russia E-mail: [email protected], [email protected] Received 08.04.2021, accepted 27.04.2021, published 30.09.2021
Abstract. The purpose is to study the mechanisms leading to the genetic divergence, i.e. stable genetic differences between two adjacent populations coupled by migration of individuals. We considered the case when the fitness of individuals is strictly determined genetically by a single diallelic locus with alleles A and a, the population is panmictic and Mendel's laws of inheritance hold. The dynamic model contains three phase variables: concentration of allele A in each population and fraction (weight) of the first population in the total population size. We assume that the numbers of coupled populations change independently or strictly synchronously. In the first case, the growth rates are determined by fitness of homo- and heterozygotes, the mean fitness of the each population and the initial concentrations of alleles. In the second case, the growth rates are the same. Methods. To study the model, we used the qualitative theory of differential equations studies, including the construction of parametric and phase portraits, basins of attraction and bifurcation diagrams. We studied local bifurcations that provide the fundamental possibility of genetic divergence. Results. If heterozygote fitness is higher than homozygotes, then both populations are polymorphic with the same concentration of homologous alleles. If the heterozygotes fitness is reduced, then over time the populations will have the same monomorphism in one allele, regardless of the type of population changes. In this case, the dynamics is bistable. We showed that the divergence in the model is a result of subcritical pitchfork bifurcation of an unstable polymorphic state. As a result, the genetic divergent state is unstable and exists as part of the transient process to one of monomorphic state. Conclusion. Divergence is stable only for populations that maintain a population ratio in a certain way. In this case, it is preceded by a saddle-node bifurcation and dynamics is quad-stable, i.e. depending on the initial conditions, two types of stable monomorphism and divergence are possible simultaneously.
Keywords: genetic divergence, differential equations, dynamics, saddle-node bifurcation, bi- and quad-stability.
Acknowledgements. This work was carried out within the framework of the state targets of the Institute for Complex Analysis of Regional Problem of the Far Eastern Branch of the Russian Academy of Sciences.
For citation: Frisman EYa, Kulakov MP. On the genetic divergence of two adjacent populations living in a homogeneous habitat. Izvestiya VUZ. Applied Nonlinear Dynamics. 2021;29(5):706-726. DOI: 10.18500/0869-6632-2021-29-5-706-726
This is an open access article distributed under the terms of Creative Commons Attribution License (CC-BY 4.0).
Введение
Одной из принципиальных и интересных задач математической популяционной генетики и биофизики является описание механизмов и определение условий возникновения первичной генетической дивергенции - устойчивого различия генетических структур в системе изначально однородных популяций. Математическому моделированию первичной генетической дивергенции много внимания было уделено в 70-80 годы прошлого века и тогда казалось, что все основные идеи здесь уже высказаны и все основные результаты уже получены [1-5]. Однако за последние 40 лет произошло бурное развитие модельных подходов и класса решаемых задач, а также методов анализа нелинейных динамических систем [6,7]. Создано множество программных средств, позволяющих вывести исследования на качественно новый уровень, быстро и точно провести этот анализ, ярко и наглядно представить его результаты (см., например, [8]). Как и прежде, современные исследователи для изучения движущих факторов возникновения генетической дивергенции используют как системы обыкновенных дифференциальных уравнений [9-11], так и разностные уравнения (отображения) [10, 12-17]. Развитие вычислительной
техники сделало возможным выполнение множества численных экспериментов с использованием соответствующих имитационных моделей генетического состава популяций [15,18]. Вместе с тем пришло понимание того, что отбор и эволюция могут идти с сопоставимыми экологическим процессам скоростями. В результате активно стали развиваться модели, учитывающие генетический дрейф и миграцию особей между популяциями, заселяющими разветвленные сети связанных участков [16]; модели, учитывающие частотно-зависимый отбор, материнский генетический эффект [17], плотностно-зависимые факторы [13,19], а также хищничество [11,20] и паразитизм [12]. В связи с этим можно указать несколько интересных обзоров популяционно-экологических и популяционно-генетических моделей, учитывающих расселение и другие факторы, приводящие, в частности, к генетической дивергенции [10,19,21,22].
Поэтому нам представляется, что настало время воспользоваться новыми открывшимися возможностями и пересмотреть модели первичной генетической дивергенции, применив современный арсенал методов и средств анализа динамических систем. В настоящей работе рассматривается математическая модель, основанная на дифференциальных уравнениях и описывающая изменение частот аллелей и соотношения численностей в системе двух смежных панмектичных популяций, связанных изотропной миграцией. Методами бифуркационного анализа исследуются возможности, условия и механизмы формирования устойчивого различия генетических структур рассматриваемых популяций.
Для понимания основных закономерностей микроэволюции популяции диплоидных организмов под действием естественного отбора ограничимся подробным рассмотрением простой модельной ситуации, когда все адаптивное разнообразие в популяции определяется одним диал-лельным локусом с аллеломорфами А и а, причем фенотип особи жестко определяется ее генотипом; популяция панмиктична и в ней действуют менделевские правила наследования. В этом случае действие естественного отбора можно описать количественно, поставив в соответствие каждому из генотипических классов особей АА, Аа и аа по одному коэффициенту. Для разных моделей действие отбора характеризуют либо величинами waa, WAa и waa, называемыми приспособленностью (fitness) особей соответствующего генотипа, либо величинами saa = WAA — 1, SAa = WAa — 1 и saa = waa — 1, называемыми коэффициентами отбора. Для непрерывно и равномерно размножающейся популяции коэффициенты отбора Sij равны сумме среднего числа выживших потомков, произведенных одной особью данного генотипического класса в единицу времени и доли выживших особей этого генотипа за то же время [2, 10]. Сделанных предположений достаточно, чтобы получить следующее дифференциальное уравнение, описывающее динамику частоты q аллеля А:
В математической популяционной генетике это уравнение, описывающее динамику генетической структуры, изначально рассматривалось независимо от динамики численности популяции [1,9]. Более того, численность популяции считалась постоянной и независимой от генетической структуры. Исследователей интересовало только изменение генетической структуры, которое, как они считали, описывается уравнением (1) при неизменной численности популяции.
Такой подход требует некоторых дополнительных обсуждений. Если исходить из приведенного определения коэффициентов отбора, то можно получить следующее уравнение для динамики численности популяции
1. Краткое описание модели
dq dt
q(1 — q)(sAa — Saa + (sAa + saa — 2sAa)q).
(1)
(2)
где 8 = 8ааЧ2 + 2вАаЧ(1 — д) + 8аа(1 — д)2 - среднее значение коэффициента отбора в популяции в данный момент времени. Заметим, что
^ = 2((вАА + 8аа — 2вАа)Я + ЯЛа — Яаа),
поэтому
¿д д(1 — д) ¿8 2 йд
йз йз йд д(1 — д) / 4 2
\Лд)
М йд М 2 \ йд
Таким образом, при фиксированных коэффициентах отбора численность популяции изменяется по экспоненциальному закону с переменным коэффициентом роста (мальтузианским параметром), который в результате естественного отбора может лишь возрастать, стремясь к постоянному значению стационарных положениях генетической структуры. Для того чтобы численность оставалось неизменной, требуется, чтобы мальтузианский параметр был равен нулю в каждый момент времени. Добиться этого можно, если постоянные значения коэффициентов отбора 8^ заменить переменными значениями, равными их отклонениям от средних: = 8^ — 8. В этом случае уравнение (1) не изменится, а уравнение (2) преобразуется к виду (Ш/(И = 0.
Рассмотрим теперь две панмиктичные менделевские однолокусные диаллельные популяции, расположенные на однородном участке ареала и связанные миграционными потоками. Предположим, что особи каждой из этих популяций, имеющие одинаковый генотип, не отличаются по приспособленностям, но частоты аллелей в популяциях могут различаться. Предположим далее, что популяции обмениваются мигрантами и интенсивности миграционных потоков пропорциональны численности той популяции, оттуда эти мигранты проистекают. Обозначим коэффициент миграции (доли мигрантов в единицу времени) через т, ^ 0.
Будем описывать эволюцию такой системы популяций с помощью следующих динамических переменных: д\ и д2 - концентрации (частоты) одного из аллелей (для определенности аллеля А) в первой и второй популяции, соответственно (0 ^ д\ ^ 1, 0 ^ д2 ^ 1), а также веса одной из популяций (например, первой) р = ^/(N1 + N2), где N1 и N2 - численности популяций (0 < р < 1).
Сформулированных предположений и допущений вполне достаточно, чтобы получить следующую систему дифференциальных уравнений, описывающих процесс эволюции двух популяций, обитающих на однородном по отбору ареале и связанных миграциями:
^ = 41(1 — — 8аа + (вАА + — 28Аа)Я1) + Ш --- (д2 — Яг),
аЪ р
= 42(1 — д2)(вАа — ваа + ^АА + 8аа — 28Аа)д2) + (д1 — д2), (3)
йъ 1 — р
= р(1 — р)(^1 — 8) + т(1 — 2р),
где 81 = 8ААд{ + 28Аад1(1 — д\) + 8аа(1 — д1)2 и 82 = 8ААд2 + 28Аад2(1 — Ы + ваа(1 — Я2)2 -средние значения коэффициента отбора в первой и второй популяции, соответственно. Более подробный процесс вывода уравнений (3) можно найти в монографии [5].
Конкретизируем тип отбора. Нас интересует возможность и условия формирования генетической дивергенции (устойчивого различия генетических структур) в системе двух смежных
и
популяций, обитающих на однородном ареале. Из предположения об однородности ареала следует, что значения приспособленностей и, соответственно, коэффициентов отборов генотипов (вАА = ^аа — 1, <ЪАа = г^Аа — 1 и заа = таа — 1) в рассматриваемых популяциях не различаются. Возникает вопрос: какой тип отбора может привести к генетической дивергенции? Движущий отбор, когда вдА > вАа > $аа или вдА < вАа < $аа, здесь явно не подходит. Значит, надо разбирать ситуацию, когда приспособленность гетерозиготы лежит вне диапазона приспо-собленностей гомозигот. Ограничимся рассмотрением «симметричного» случая 'Маа = и положим вАА = 8аа = 0, ЗАа = 5. Следует отметить, что именно такая ситуация исследуется в большинстве современных работ; обычно полагают, что в < 0 [14-16].
Подставляя эти значения коэффициентов отбора в систему (3), получаем:
^ = ,вд1(1 — дг)(1 — 2дх) + т -—- (д2 — дг), М р
^ = зд2(1 — д2)(1 — 2д2) + (дх — д2), (4)
йъ 1 — р
— = 2вр(1 — р)(д1 — д2){1 — дх — д2) + т(1 — 2р).
Из биологических соображений имеет смысл говорить только о тех решениях системы (4), которые целиком лежат в единичном кубе (0 ^ дх ^ 1, 0 ^ д2 ^ 1 и 0 < р < 1), поскольку фазовые переменные - это доли от целого. Однако для полного понимания особенностей динамики необходимо изучить, на первый взгляд, биологически незначимые типы решений, так как их изменение сказывается качественно на значимых решениях.
2. Исследование устойчивости
Остановимся подробнее на исследовании режимов динамики данной модели. Система (4) имеет следующие стационарные состояния:
1. Ео(0, 0,1/2) - обе популяции представлены только особями с генотипом аа (отсутствуют генотипы АА или Аа - мономорфная популяция);
2. Е\(1/2,1/2,1/2) - обе популяции включают особей всех трех генотипов при одинаковой концентрации каждого из двух аллелей (полиморфная популяция);
3. Е2(1,1,1/2) - обе популяции состоят только из особей с генотипом АА (мономорфная популяция);
„ ^ /1 , Vs2 + 4ms 1 Vs2 + 4ms А „
4. пара точек Ь3,4 - ±-, - ^-, - , которая существует при s < 0 и
\ 2 2s 2 2s 2 I
s + 4т < 0, либо при s > 0 и s + 4т > 0. Однако лишь в части двух этих диапазонов точки равновесия будут иметь биологически значимые координаты (лежать в единичном кубе). Так, в первом случае пара Е3,4 лежит в единичном кубе, только если т ^ 0 (область 4 на рис. 1, a), а во втором - если т < 0 (область 8 на рис. 1, a). Соответственно несложно показать, что на линии т = —s/4 рождается (или исчезает) пара точек Е3,4 вследствие бифуркации вил (pitchfork bifurcation). В этом случае точка Е\ «расщепляется» на две дополнительные точки Ез и Е4. Каждая из этой пары точек соответствует дивергентному состоянию системы популяций: в одной из полиморфных популяций преобладает один из аллелей: например, достаточно высокая концентрация аллеля А и низкая - аллеля а, в то время как во второй популяции преобладает другой аллель: высокая концентрация аллеля а и низкая - аллеля А.
Для подробного исследования системы (4) рассмотрим величину Q(t) = д1(Ь)р(Ь) + + д2(Ъ)(1 — р(Ъ)), равную концентрации аллеля А во всей системе популяций. Дифференцируя ее, получаем:
dQ йд 1 ^2 , , Лр
= Р~Ж + (1 — р) + 1 — ^ м-
Подставляя в это выражение значения производных из системы (4), получаем следующее дифференциальное уравнение:
^ = — 1)(С — Q), (5)
где С = р(д2 — д?2) + дИз (5) следует, что на поверхности Q = 1/2 производная общей концентрации dQ/dt = 0. Это, в частности, означает, что ни одна из модельных траекторий не может пересечь эту поверхность, и они всегда остаются расположенными по одной из ее сторон. Таким образом, Q = 1/2 - это сепаратрисная поверхность, разделяющая режимы с разными соотношениями концентрации аллеля А. Выше этой поверхности Q > 1/2, ниже Q < 1/2. Поверхность Q = С, обладающая схожими свойствами, расположена вне единичного куба и лишь касается части его ребер (д 1 = 0,1, д2 = 0,1, р е К).
В тоже время из (5) следует, что из любой начальной точки (в том числе вне единичного куба), удовлетворяющей соотношению Q = 1/2, формируется траектория, целиком лежащая на этой поверхности, а сама поверхность является центральным многообразием. Это означает, что на ней динамика системы (4) в точности до малых может быть описана системой двух дифференциальных уравнений, то есть можно выполнить процедуру понижения размерности системы. Обозначим эту поверхность символом Ш, а также будем добавлять к нему верхний индекс 5 или SS для обозначения устойчивости траекторий, лежащих на ней (для которых верно dQ/dt = 0). Соответственно, Ш3 - устойчивое многообразие - такое, что слабые возмущения выбивают траекторию с этой поверхности; Ш33 - безусловно устойчивое многообразие - такое, что даже сильно возмущенная траектория вернется на Ш33 и будет двигаться по ней, пока не достигнет устойчивой точки (аттрактора). Очевидно, что характер устойчивости этого многообразия определяется типом устойчивости особых точек, лежащих на ней, и размерностью их собственных подпространств. Несложно заметить, что точки Е1, Е3 и Е4 лежат на поверхности Q = 1/2.
Для исследования устойчивости состояний равновесия системы (4) будем вычислять собственные значения стандартным образом. В результате в плоскости параметров 8 ит, несложно выделить несколько областей, соответствующих разному типу устойчивости каждой особой точки. Опишем их.
Для начала отметим, что при т ^ 0, то есть в областях, указанных на рис. 1 под номерами 1-4, существуют устойчивые стационарные состояния, лежащие в единичном кубе. При этом единичный куб оказывается притягивающим, то есть любая траектория, которая имеет начало даже за пределами куба (из некоторой области, ограниченной поверхностью Q = С), стягивается к одной из стационарных точек Е± (г = 0,1,... 4). В областях 5-8 это не так. Там даже если существуют устойчивые состояния равновесия, то лежат они за пределами единичного куба (области 5 и 6). В результате траектория, рожденная в кубе, покинет его, а система (4) потеряет биологический смысл.
В областях под номерами 1-2 точка Е1 - устойчивый узел. Поэтому любая траектория из единичного куба (даже некоторой области притяжения вне него) сходится к этой точке. Однако, за счет того, что области 1-3 отличаются количеством и типом устойчивости других особых точек, качественно модельные траектории оказываются достаточно разнообразными в разных областях. Так, в области 1 точки Ео и Е2 - седловые с одномерным устойчивым многообразием, которое совпадает с ребром куба д1 = д2 = 1 или д1 = д2 = 0 (рис. 1, Ь). В то время как пара
л
S = 2, /77 = 0.5
1© 0.5- \m s/2 © ' © I 1 s
-2®-1 © a -0.5-1 . w -s/4
=0.5
/
W \ s=-1, m = 0.(
Рис. 1. а - Параметрический портрет системы (4) и Ь-€ - фазовые портреты системы (4) в биологически значимых областях 1-4
Fig. 1. a - Parametric and b-e - phase portraits of model (4) in biologically significant domains 1-4
1
d
e
точек - неустойчивые узлы, из которых исходят две пары ветвей - неустойчивое многообразие Wи (одномерная сепаратриса), которое втекает в точку Е\. В результате любая траектория, лежащая между двумя ветвями Wи (в некоторой от них окрестности), стремится к устойчивому узлу , даже если начальная точка лежит вне единичного куба. Напротив, если начальная точка лежит выше или ниже соответствующей ветви и вне единичного куба, то траектория оказывается неограниченной. В области 2 при пересечении линии m = s/2 в окрестности седловых точек Е0 и Е2 устойчивое многообразие оказывается двумерным, и оно «отсекает» значительно большую часть фазового пространства, чем в области 1, и «добавляет» к области притяжения точки Е\ большую область вне единичного куба (рис. 1, c). Кроме того, неустойчивое многообразие Wи с началом в точках круто изгибается вверх или вниз (верхняя и нижняя ветви) и уходит в бесконечность при удалении от линии m = s/2. Таким образом, область притяжения полиморфной точки Е\ в области 2 значительно больше, чем в области 1.
При переходе в область 3 точек не существует, точка Е\ теряет устойчивость (становится седловой), а устойчивость приобретают точки Е0 и Е2 (рис. 1, d). В этом случае динамика системы становится бистабильной. В результате для начальных условий, лежащих по разные стороны от сепаратрисного многообразия Ws, траектория стремится к одной из мономорфных точек: при Q > 1/2 к Е2, при Q < 1/2 к Е0. Точки Е0 и Е2 в этом случае оказываются глобально устойчивыми и притягивающими для любых начальных условий. В области 4 система также би-стабильна с тем отличием, что при пересечении линии m = -s/4 рождается пара седловых точек
Е^^, целиком лежащих в единичном квадрате (рис. 1, в). Устойчивые собственные подпространства этих точек строго лежат на поверхности Q = 1/2, а неустойчивое порождает неустойчивое многообразие Ши, которое изгибается вверх или вниз для точек Е3 и Е4 в зависимости от значения Q (см. рис. 1, e). Наличие в двух этих областях параметров такой топологической структуры фазового пространства, определяемое размерностью собственных подпространств, означает, что в ряде случаев формируются довольно причудливые немонотонные модельные режимы динамики. Например, в области 3, если начальная точка расположена на поверхности многообразия Ш^ или достаточно близка к ней, то траектория медленно движется по ней, пока не достигнет седла Е\. В этом случае слабое возмущение быстро «вытолкнет» траекторию в сторону Е0 или Е2, которых она достигнет достаточно быстро. Концентрации ql и д2, а также вес р при этом будут изменяться довольно плавно.
В области 4 траектории будут немного сложнее. Например, траектория, лежащая на поверхности многообразия Ш^, будет двигаться по ней, пока не достигнет точки Е3 или Е4. Точка Е\ в этом случае порождает сепаратрису на Ш^ такую, что, если начальная точка лежит правее или левее от нее, то траектория стремится к Е3 или Е4. После чего слабые возмущения «вынуждают» траекторию двигаться к одному из мономорфных состояний Е0 или Е2. Траектория при этом оказывается расположена достаточно близко к неустойчивому многообразию Ши (см. рис. 1, e). Наличие таких траекторий (да и самого Ши) означает, что должен быть отмечен скачок значения веса р. Другими словами, по мере приближения траектории к одному из мономорфных состояний происходит стремительный рост (или падение) численности одной из субпопуляций, а после вес падает (растет), пока не достигнет стационарного значения 1/2 в точке Е0 или Е2. Рис. 2 демонстрирует примеры таких режимов динамики.
На рис. 2 показаны четыре разные траектории при различных незначительно отличающихся начальных условиях (координаты стартовых точек приведены на рис. 2, a). Для первой начальной точки строго выполняется Q = 1/2, то есть она лежит на многообразии Ш^. Формируемая при этих начальных условиях траектория с точностью до погрешности численного интегрирования движется по поверхности Ш^. Далее она ненадолго (снова же, из-за
WS E3 WU
1
41 А E2
/ 1
0.5 ( Р 42 А ; / /
0.5
50
150 200
41 Р 1 ---- E2
42
41 e2
/\ у "Г s Р
"q'2""
10 20 t 30 40 50
(b)
(c)
(d) (e)
(qi(0), <fc(0), p(0)) (0.6, 0.2, 0.75) (0.61, 0.21, 0.751) (0.601, 0.205, 0.751) (0.599, 0.201, 0.751)
10 20 t 30 40 50
10 20 t 30 40 50
Рис. 2. Примеры режимов динамики в области 4 при s = -0.9, т = 0.15 и начальных условиях (^1 (0), Я2(0),р(0)), указанных на рисунке
Fig. 2. Examples of dynamics regimes in domain 4 at s = -0.9, rn = 0.15 and initial conditions (qi (0),^2(0),p(0)) shown on the figure
0
t
b
1
0
0
c
e
погрешности интегрирования) задерживается в окрестности точки Ез, пока возмущения не выбьют траекторию в окрестность неустойчивого многообразия Ши. В этот момент величина р возрастает, и по мере приближения к устойчивой точке Е2 ее рост прекращается и сменяется падением до значения 1/2 (рис. 2, Ь). Для большей наглядности на рис. 2, Ь выделены перечисленные этапы динамики, соответствующие движению вдоль Ш3 и Ши и стационарной динамике в мономорфном состоянии Е2 или Ео (рис. 2, в). Две другие стартовые точки отличаются от первой на сотые. Это приводит к тому, что траектория проходит лишь в окрестности многообразий Ш^ и Ши и намного быстрее достигает мономорфных состояний Е2 или Ео, однако небольшой скачок величины р по-прежнему наблюдается (рис. 2, с).
На бифуркационной линии т = 0, то есть на границе, отделяющей биологически значимые области параметров 1 и 4 от незначимых областей 8 и 5, система (4) вырождается: первые два уравнения не зависят друг от друга и от третьего уравнения. Очевидно, что стационарные значения для первого и второго уравнения системы (4) равны 0, 1/2 и 1. При 8 > 0 устойчива точка д = 1/2, при 8 < 0 - точка д = 0 и д = 1. Вместе с тем третье уравнение зависит от первых двух следующим образом. Если д\ и д2 достигли одного из своих стационарных значений ((1/2,1/2) при § > 0, (0, 0), (0,1), (1, 0) или (1,1) при § < 0), то третье уравнение обращается в нуль независимо от значения р. Как следствие, величина р не достигает своего стационарного значения - 0 или 1 (это становится возможным, только если р(0) = 0 или 1). В результате при т=0 стационарное значение р (соотношение численностей на обеих территориях в асимптотическом случае) может быть любым и зависит от начальных значений д\(0), д2(0) и р(0). Это означает, что в системе (4) при т = 0 и в > 0 существует прямая устойчивых особых точек: (1/2,1/2, р) (рис. 3, а). Соответственно при < 0 таких прямых уже четыре: (0, 0,р), (0,1,р), (1, 0,р) и (1,1,р), где р - любое действительное число (рис. 3, Ь). В последнем случае область допустимых начальных значений переменных (единичный куб) разбивается на четыре равные бассейна притяжения, ограниченные плоскостями д\ = 1/2, д2 = 1/2 и гранями куба. При 0 < д1(0) < 0.5 и 0 < д2(0) < 0.5 переменные д\ и д2 стремятся к 0, при 0.5 < дг(0) < 1 и 0.5 < д2(0) < 1 - к 1. При 0 < дг(0) < 0.5 и 0.5 < д2(0) < 1 частота аллеля А в первой популяции - д\ стремится к 0, а во второй - д2 - к 1, при 0.5 < дг(0) < 1 и 0 < д2(0) < 0.5 частота аллеля А в первой популяции стремится к 1, а во второй - к 0. Две последние ситуации можно интерпретировать как генетическую дивергенцию популяций. Причем все это никак не зависит от начального значения р(0).
Наконец, в биологически незначимых областях 5 и 6 на рис. 1, а устойчивыми оказываются точки Ез и Е4, но они лежат за пределом единичного куба, а концентрации д\ или д2 оказываются отрицательными или больше единицы (области 5 и 6 отличаются размерностью устойчивого многообразия в окрестности точек Ео и Е2). В областях 7 и 8 все особые точки лежат в единичном кубе, но среди них нет устойчивых (области 7 и 8 отличаются количеством особых точек).
Обсудим содержательные биологические результаты, которые следуют из проведенного исследования. При этом нас, в первую очередь, естественно, интересуют области значений параметров, соответствующие биологическому содержанию задачи: 8 > —1 и 0 < т < 1. Проведенное исследование показывает, что при этих значениях параметров в рассматриваемой системе двух популяций, связанных миграциями, никакой устойчивой генетической дивергенции быть не может. При положительных значениях параметра 8, то есть в случае повышенной приспособленности гетерозигот (области 1 и 2 на рис. 1, а), в каждой из популяций устанавливается, естественно, одинаковый устойчивый полиморфизм с равными концентрациями аллелей А и а. При отрицательном значении параметра 8, то есть в случае пониженной приспособленности ге-терозигот (области 3 и 4 на рис. 1, а), в каждой из популяций устанавливается одинаковый устойчивый мономорфизм либо по генотипу АА (при Q > 1 /2), либо по генотипу аа (при Q < 1 /2).
(0.5, 0.5, p)
-0.5
(qi(0), q2(0), p(0))
1 - (0.7, 0.99, 0.2)
2- (0.001, 0.3, 0.8)
\Я2
(qi(0),q2(0), p(0)) (0.49, 0.89, 0.98) (0.49, 0.89, 0.95)
401 0.5f->
/Р
0.5
0
a
0
10 t 20
30
0
b 0
10 20 t 30 40 50
Рис. 3. Примеры динамики системы (4) при т = 0, s = -0.5 (a) и s = 1 (b)
Fig. 3. Examples of dynamics regimes of system (4) at m = 0, s = -0.5 (a) and s = 1 (b)
0
1
1
2
При относительно малых значениях коэффициента миграции (0<т< —в/4, область 4 на рис. 1, a) в системе появляются положения равновесия , соответствующие генетической дивергенции популяций, когда в одной из полиморфных смежных популяций преобладает один из аллелей, а во второй - другой, но эти положения равновесия оказываются неустойчивыми для любых значений параметров модели, при которых они существуют (за исключением вырожденного случая т = 0). Как уже указывалось, эти точки лежат на поверхности Ш^ (при Q = 1/2), и любые малые возмущения, выводящие из этой поверхности, приводят к стремлению к мономорфному равновесию всей системы в целом. Содержательно это можно интерпретировать следующим образом. Если начальные условия таковы, что Q(0) = 1/2, то отбор против гетерозигот при одинаковой приспособленности гетерозигот никак не сможет нарушить это соотношение. Однако, если при этом начальные частоты аллелей в каждой из популяций не равны 0.5, то тот же отбор против гетерозигот будет способствовать увеличению доли преобладающего аллеля в данной популяции. Этому процессу будет противодействовать миграция, которая сглаживает различия в генетических структурах популяций. Если интенсивность миграции достаточно велика ( т > — в/4), различия между популяциями полностью сглаживаются и частоты аллелей в каждой из популяций выравниваются (стремятся к 0.5). При меньших интенсивностях миграций (т < — в/4) в системе превалирует дизруптивный отбор и устанавливается одно из двух возможных дивергентных
равновесий Ез или Однако такое равновесие оказывается возможным, только если ф(0)=1 /2, и оно неустойчиво по отношению к любым возмущениям, нарушающим это равенство. Выход из поверхности Ш^ приводит к тому, что либо в каждой из связанных популяций будет преобладать один и тот же аллель (и в условиях бистабильности он будет стремиться к фиксации), либо в популяциях преобладают разные аллели, но частоты преобладающих аллелей различаются (поскольку = 1/2). Но тогда популяция, имеющая большую частоту преобладающего аллеля, будет иметь и большую среднюю приспособленность. В силу этого она будет преобладать по скорости роста, а затем и по численности (весу), и по потоку мигрантов («эмигрантов»). В результате преобладающий в одной из популяций аллель со временем окажется преобладающим во всей системе популяций. В итоге это приведет к мономорфизму и выравниванию численностей (весов) за счет миграций.
Заметим, что при полной изоляции (то есть при т = 0) генетическая дивергенция популяций, естественно, возможна при пониженной приспособленности гетерозигот Аа (в < 0). В этом случае каждая отдельно взятая популяция всегда оказывается мономорфной (присутствует только генотип АА или аа). Генотип на разных территориях при этом может полностью совпадать (наблюдается одинаковый устойчивый мономорфизм) или быть различным (дивергенция). Вместе с тем, из-за различий в начальных значениях средних приспособленностей особей, финитное отношение численностей популяций оказывается любым (0 < р < 1) в зависимости от начальных значений частот генов и соотношения численностей (веса ). Примеры динамики на рис. 3 демонстрируют тот факт, что в асимптотическом случае концентрации аллеля А достигают значения 0, 1/2 или 1, в то время как величина веса р каким-то образом зависит от начальных условий и редко достигает значения 1/2.
3. Частный случай системы (4) при р = 1/2
Несложно заметить, что в невырожденном случае (m = 0) любое из стационарных состояний Е0-Е4 характеризуется значением веса р = 1 /2, то есть обе популяции в асимптотическом случае оказываются равными по численностям независимо от того, какие частоты аллелей А и а соответствуют этому состоянию. В силу этого интересно проанализировать, как будет меняться концентрация аллелей в популяциях с постоянными численностями или с постоянными соотношениями численностей, то есть при р = const. Для начала рассмотрим простейшую ситуацию: р = 1/2 - численности равны на обеих территориях. В этом случае система (4) имеет вид [2,3]:
^ = sqi(1 - qi)(1 - 2qi) + m(q2 - q 1),
dl (6) — = sq2(1 - q2)(1 - 2q2) + m(q 1 - q2).
Система (6) имеет тот же набор особых точек, что и (4): Е0-Е4. Однако к ним добавляются еще четыре точки Е5-Е8, которые находятся как корни полинома 8s2q4-16s2qf+2s(5s-2т)q\ + + 2s(2m-s)qi + m(2m-1), азначения q2 находятся по формуле: q2 = qi(-2sq2+3sqi+m-s)/m. Описание этих особых точек достаточно громоздко, поэтому не будем приводить здесь их явный вид. Однако несложно показать, что этот полином имеет строго четыре действительных корня при s < 0, если s/2 < m < -s/6, а при s > 0, если -s/6 < m < s/2. На границах этого интервала, то есть при m = ±s/2 или ±s/6, полином имеет только два корня второй кратности, совпадающих либо с точками Ео и Е2 при s > 0, либо с Е3 и Е4 при s < 0. Соответственно в этих диапазонах параметров точки Е5-Е8 существуют, а появляются они в результате бифуркации вил из пары Ео,2 или Е3 4.
m= 0.25, 5= - 0.25
1
42 0.5
m
m=0.25, 5=0.25
0 0.5 q1 1 m= 0.16, 5= - 0.8
1
q2 0.5
3 2
q2 0
-0.5 m= - 0.2, 5= - 0.1
0 u 0.5 q1 1 m= 0.1, 5= - 0.9
-2-о
41 -
-2 0 q1 2 3 m= - 0.2, 5= - 0.7
1.5 1
q2 0 0.5
-1 0 q1 1 2 m=0.1, 5=0.35
m= - 0.18, 5=0.9
-0.5 0 q1 1 1.5 m=- 0.1, 5=0.9
1
q2 0.5
/J®
inf-
1
q2 0.5
0 0.5 q1 1
-1 0 q1 1 2 0 0.5 q1 1
0 0.5 q1 1
Рис. 4. Параметрический портрет системы (6), а также фазовые портреты, совмещенные с бассейнами притяжения устойчивых стационарных состояний. Серые области (inf) на бассейнах под номерами 1, 1', 2 и 8 соответствуют неограниченным модельным траекториям. В области 7 все траектории неограниченные (уходят на бесконечность). В остальных случаях белые и серые области - бассейны сосуществующих устойчивых состояний равновесия
Fig. 4. Parametric and phase portraits of the system (6) which are combined with the basins of attraction of stable equilibrium. Gray areas (inf) with numbers 1, 1', 2 and 8 correspond to unbounded model trajectories. In area 7 all trajectories are unbounded (seek to infinity). In other cases, white and gray areas are basins of coexisting stable equilibrium
0
0
0
0
Области существования и устойчивости особых точек Е0—Е4 совпадают с аналогичными областями системы (4), а для точек Е5-Е% они оказываются «вложены» в них. В верхней части рис. 4 в плоскости параметров показаны эти области. Рассмотрим их подробнее.
В области под номером 1, а также в области 1', не имеющей биологического смысла (по отрицательному значению параметра т, а не значениям фазовых переменных), существует девять точек, Ео-Е$. В области 2 и 8 - пять, Е0-Е5. Устойчивым является полиморфное состояние Е1. Интересно, что при положительных значениях параметра т некоторые особые точки расположены за пределом первого квадранта. Поэтому область притяжения полиморфной точки Е1 включает как весь единичный квадрат, так и зоны, расположенные за его пределами. Однако при отрицательных значениях т все стационарные точки и бассейн притяжения Е1 расположены в единичном квадрате.
В области 3, а также в биологически незначимой области 6 динамика системы (6) становится бистабильной. Здесь существуют две особые устойчивые точки - мономорфные Ео и Е2 с областями притяжения, разделенными прямой д1 = 1 — д2 (в области 3) и д1 = д2 (в области 6).
В области 4 появляется пара точек Ез, которые кардинально не влияют на области притяжения мономорфных состояний равновесия, но вносят в динамику особенности, аналогичные тем, что были в области 4 трехмерной системы (4). В этом случае модельные траектории содержат немонотонные участки динамики. Если сумма начальных концентраций д 1(0) + (¡2(0) = 1 (или достаточно близко к 1), то сначала траектория движется вдоль сепаратрисы д\ = 1 — д2, пока не достигнет окрестности точки Е3 или Е4 (соответствующих дивергенции, где она может «задержаться»), а затем вследствие малых возмущений она устремляется к точке Ео или Е2. В результате вначале величина будет расти, а 2 падать или наоборот. Затем и 2 будут одновременно расти при движении к Е2 или падать при движении к Ео.
Наконец, в областях 5 и 5' от каждой седловой точки Е% и Е4 отщепляется по паре Е§,6 и Е7,8, а полиморфные точки Ез и Е4 приобретают устойчивость. В результате в этих областях помимо устойчивых мономорфных состояний существуют устойчивые состояния равновесия, соответствующие полиморфной популяции с высокой концентрацией аллеля А и низкой аллеля а или наоборот. Таким образом, динамика системы (6) оказывается квадростабильной, то есть в зависимости от начальных условий величины и 2 с течением временем (при достаточно большом ) стремятся к четырем принципиально разным значениям: либо мономорфным состояниям с нулевой концентрацией аллеля А или а, либо к полиморфным с неравными концентрациями (с преобладанием того или иного аллеля).
Наконец, в биологически незначимой области 7 нет устойчивых особых точек, а все модельные траектории оказываются неограниченными.
Таким образом, генетическая дивергенция в системе миграционно связанных популяций -наличие устойчивых состояний равновесия, соответствующих полиморфной популяции с высокой концентрацией аллеля А и низкой аллеля а, или наоборот - оказывается возможной в модели, соответствующей системе (6) при пониженной приспособленности гетерозигот ( < 0) и достаточно малом коэффициенте миграций: т < — в/6 (область 5 на рис. 4).
Предположим теперь, что численности на разных территориях различные, но в целом сохраняется отношение численностей р = N\/(N\ + N2) = const, такое что 0 < р < 1. Будем считать величину р параметром, который если и изменяется, то достаточно медленно. В этом случае систему (4) перепишем в виде:
Как и прежде, система (7) имеет от трех до девяти особых точек, которые в этом случае зависят еще от параметра р, за исключением точек Ео(0, 0), Е1 (1/2,1/2) и Е2(1,1). Остальные шесть точек Е3-Е8 при р = 1/2, очевидно, совпадают с описанными в предыдущем разделе. В данном случае модель (7) - трехпараметрическая. Поэтому изменение характера устойчивости стационарных состояний необходимо выполнить при изменении всех трех параметров.
Для начала несложно показать, как изменяются координаты устойчивых и неустойчивых особых точек, например, при изменении параметра 8 и фиксированных значениях р и т. Рис. 5 демонстрирует это, а также показывает механизм рождения всех особых точек. Очевидно, что простые точки Ео, Е\ и Е2 существуют при любых значениях параметров (области 1-У на рис. 5), в то время пара точек Е3,4 с координатами вне единичного квадрата - лишь при 5>0
4. Обобщение модели (6) для любых р
(7)
pf*
е
En
-2 E
4e3
E£r.
m=0.15,p=0.5
®
En
e6
E
6
e5
2 -2 /
E8
qi
■ , ______2 ,
E
m=0.15, p=0.49
H^qi^1 I' m=0.05, p=0.35 -__ез_________ 1
tc e
4 eg
■4 J-I о
--------E7____
E6 ek
'"""ЕГ1
qi
Е^_ Е2 ч Ei E3
:-Ез sn 1 pf sn Ei \
6
0.5 Еп , , /5, ^ЕГ i
-2 Е' - Е8 V -1 IV 0 III II
m=0.1, p=0.51
Е2 \ 4-el_
е5 sn 1 Ез ei
е8 Е 4 pf Е6 Е
■2 е© -1 IV е0 0 III — II Е4 i
6
d
Рис. 5. Зависимость координат особых точек системы (7) от параметра s (бифуркационная диаграмма) с отмеченными точками бифуркации (PF - бифуркация вил, SN - седло-узловая и TC - транскритическая бифуркация). Сплошная линия q1 (s) - полностью устойчивая особая точка, пунктир - неустойчивая
Fig. 5. Coordinates of equilibrium points of system (7) depending on the parameter s (bifurcation diagram), where bifurcation points are marked. PF, SN and TC is pitchfork, saddle-node and transcritical bifurcation respectively. Solid and dotted lines q1 (s) are stable and unstable equilibriums respectively
a
c
s
(область I и II). При < 0 пара с теми же обозначениями Е34 рождается путем расщепления полиморфного состояния Е (1/2,1/2) на два, то есть в результате бифуркации вил (область IV на рис. 5, a). При в > 0 пары Е5,6 и Е7,8 появляются в результате расщепления мономорфных состояний равновесия Ео и Е2 (область I на рис. 5, a). Как видно из рисунка, часть из них лежит вне единичного квадрата. Однако, как было показано ранее, появление таких, на первый взгляд, биологически незначимых стационарных точек сигнализирует об изменениях характера динамики определенного рода, а не о потере смысла математической модели. В данном случае при больших значениях устойчивой является точка Е в центре единичного квадрата. Кроме нее на фазовой плоскости существует пять или девять состояний равновесия, которые соединены устойчивыми и неустойчивыми многообразиями (сепаратрисами). Это указывает на существование изолированных траекторий определенного рода. В бассейнах притяжения даже единственной устойчивой полиморфной точки Е (область I, I', II на рис. 6) несложно выделить подобласти, в каждой из которых существуют траектории, которые не способны покинуть ее, так как они ограничены сепаратрисами.
При < 0 пары точек, обозначенные как Е§,6 и Е7 8, полностью лежат в единичном квадрате (область V на рис. 5, a-d), а система также содержит изолированные траектории. Но более интересно, что механизм рождения этих точек зависит от значения веса р. Строго при р = 1/2 ив < 0 пары Е5,6 и Е7,8 - результат субкритической бифуркации вил. Это означает, что правее точки РБо на рис. 5, a неустойчивое состояние Е3 (Е4) - это седло. Левее РБо, когда появляется пара Е5,6 (Е7,8), оно приобретает устойчивость. То, как устроено фазовое пространство до и после этой бифуркации, можно проследить на рис. 4 (область 4 и 5).
m=0.1, 5=0.35, p=0.3 5
m=0.5, 5=2, p=0.1
-0.5
a
-0.5 0 1 1.5
m=0.1, 5=-1, p=0.43
b
m=0.25, 5=0.25, p=0.3
1
42 0.5
0-
d
E3
2
I
m=0.16, 5=-0.8, p=0.35 m=0.25, 5=-0.25, p=0.7
у! 1
0.5 q1 1
III l\v //E2
- EjVf
ж ,
0.5 q1 1
f
0.5 q1 1
Рис. 6. Фазовые портреты режимов динамики системы (7) из разных областей на бифуркационной диаграмме, совмещенные с бассейнами притяжения устойчивых стационарных состояний
Fig. 6. Phase portraits and basins of attraction of the system (7) for the parameters from different domains on the bifurcation diagram
При незначительной вариации параметра р характер бифуркации кардинально меняется (рис. 5, b). С одной стороны, при s < 0 графики зависимостей координат точек Е3 (Е4) и Е5,в (Е7,8) от параметра s перестают пересекаться. Это означает, что на фазовой плоскости точки Е5,6 (Е7,8) лежат левее (правее) от Е3 (Е4) (рис. 6, d), а не расположены от нее по обе стороны (рис. 6, a). В результате фазовые портреты и соответствующие бассейны притяжения оказываются менее симметричными (см. рис. 6), чем в случае р = 1/2. С другой стороны, появление пар точек Е§,6 и Е7,8 не меняет тип устойчивости Е3 и Е4 - они как были седлами, так седлами и остаются. Вместе с тем эта пара состоит из устойчивого узла ( Е5 и Е$) и седловой точки ( Еб и Е7). Это означает, что при р = 1/2 эти пары точек появляются в результате седло-узловой бифуркации (SN), а не бифуркации вил (PF). В результате при низких значениях s(s < 0, 0 < т < — s/6) система (7) оказывается квадростабильной (существуют четыре устойчивые точки).
При s > 0 наблюдается аналогичное изменение механизма рождения пар стационарных точек Е§,6 и Е78 с тем отличием, что они всегда неустойчивы, как и точки Ео и Е2, от которых они «отщепляются». В данном случае небольшое отклонение параметра р от значения 1 /2 приводит к тому, что пары Е§,6 и Е7 8, как и в случае s < 0, появляются в результате седло-узловой бифуркации, а не вил, и оказываются расположенными несимметрично относительно точек Ео и Е2 на соответствующем фазовом портрете (рис. 6, b). Поэтому в момент появления этих пар их координаты всегда оказываются вне единичного квадрата. Область значений параметров, соответствующих такой конфигурации стационарных состояний на бифуркационных диаграммах, обозначим через I' (рис. 5, b-d и рис. 7). По мере роста параметра s стационарные значения концентраций q\ и q2 растут таким образом, что точки Еб и Е7 сближаются с мономорфными состояниями Ео и Е2, проходят сквозь друг друга, испытывая транскритическую бифуркацию (TC), и вновь попадают в единичный квадрат (область I на рис. 5, b-d). В результате пары точек Е§,6 и Е78 вновь лежат по обе стороны от точек Ео и Е2, как в случае р = 1/2.
0
0
e
Описанные закономерности формирования устойчивых, а также би- и квадростабильных режимов динамики модели (7) можно изобразить в виде двумерной бифуркационной диаграммы, приведенной на рис. 7, a. Для обозначения областей параметров, лежащих между разными бифуркационными линиями на рис. 5-7, используются как арабские, так и римские цифры, которые, в общем-то, обозначают одно и то же - это области с разным устройством фазового пространства и режимами динамики. Арабские цифры в точности соответствуют предыдущему случаю при р = 1/2, а также, в определенной мере, областям параметров исходной системы (4). Римские цифры соответствуют более «общему» случаю 0 <р < 1.
Как показано выше, при разных значениях р полиморфные состояния равновесия Е^-Е^ формируются по-разному. На линии р = 1/2 они отщепляются от соответствующего мономорф-ного (Ео и Е2 при > 0) или полиморфного состояния (Ез и Е4 при < 0) и лежат достаточно симметрично относительно Ео, Е2, Ез или Е4. Это происходит в точках РБо (рис. 7, a). При р = 1/2 полиморфные состояния равновесия Е5—Е8 появляются в стороне от этих точек, когда параметры пересекают одну из ветвей бифуркационной линии Как видно из рис. 7, a,
1
0.8
Р
0.5
0.2
0 -2
m=0.2 SN I' TC
PF III II I
V IvX ...... " " ф '
J N I' TC
-1
ш PF°
SN
I' TC
V (I') Р> 1/2
E E5 (E0)
V (I) p=1/2
(E0)
V (I') p< 1/2
E5 E6 У E3 (E0)
0
1
2
I
0
1
2
b
a
Рис. 7. а - Бифуркационная диаграмма системы (7) при вариации параметров р, s и фиксированных значениях т. b - Схематическое описание механизма рождения полиморфных состояний равновесия Е5 и Е6 при переходе параметров из области IV в V (из II в I')
Fig. 7. а - The bifurcation diagram of the system (7) with variations in the parameters p, s and fixed values of m. b - Schematic of the birth of polymorphic equilibrium states E5 and E6 during the transition of parameters from region IV to V (from II to I')
эти ветви сходятся в точке PFo на прямой р = 1/2. В свою очередь при s > 0 в точке PFo, дополнительно к линиям SN, сходятся линии транскритической бифуркации TC. На рис. 7, b довольно схематично с помощью одномерного фазового портрета показан механизм рождения рассматриваемых пар полиморфных состояний Е5,6 при разных значениях р. Другая пара Е7,8 рождается аналогично, но с противоположной от точки Е4 (Е2) стороны. Вместе с тем неустойчивые полиморфные точки Е3 и Е4 при s < 0 отщепляются от точки Е2 при пересечении линии PF. Это означает, что в системе (7) имеют место бифуркации вил коразмерности 1 и 2 в зависимости от параметра .
Таким образом, в рассматриваемой системе генетическая дивергенция возможна, но только для популяций, которые, как минимум, сохраняют соотношение своих численностей (р = Ni/(Ni + N2) = const). Это является довольно сильным предположением. В любом случае постоянное отношение численностей возможно в случае их абсолютно синхронного роста или падения. В реальности синхронизация динамики связанных элементов, как правило, наблюдается при достаточно сильной миграционной связи, либо если начальные их состояния достаточно близкие. Однако проведенное исследование показывает, что дивергенция наблюдается лишь при небольших значениях коэффициента миграции (помимо условия пониженной приспособленности геретерозигот, то есть s < 0). При — s/6 < т < — s/4 (область 4 на рис. 4 или IV на рис. 7) она оказывается асимптотически неустойчивой и проявляется лишь как часть переходной динамики для специально подобранных начальных условий. Устойчивой дивергенция оказывается при достаточно малых значениях коэффициента миграции, а именно при 0 < т < — s/6 (область 5 на рис. 4 или V на рис. 7). В последнем случае система популяций оказывается квадростабильной, когда одновременно возможен мономорфизм или дивергенция в зависимости от начального соотношения концентраций зигот с геномом А или а на разных территориях. При одних соотношениях на обеих территориях устанавливается мономорфизм по генотипу А или a (рис. 8 в центре). При других на смежных территориях преобладают разные аллели -
Рис. 8. Бассейны притяжения (слева) решений системы (7) при s = -0.75, гп = 0.05, р = 0.66041 (а) и р = 0.33959 (b) (достаточно близко к границе SN). Примеры динамики, демонстрирующие переход к мономорфизму (посередине) и дивергенции (справа). Стартовые точки выбраны достаточно близко к границам бассейнов
Fig. 8. Basins of attraction (on the left) of system (7) at s = -0.75, rn = 0.05, p = 0.66041 (а) and p = 0.33959 (b) (close to the SN boundary). Examples of dynamics demonstrating the transition to monomorphism (in the middle) and genetic divergence (on the right)
на одной больше генотипа А, на другой практически мало генотипа а (рис. 8 справа) или наоборот. Размер областей притяжения состояний равновесия, соответствующих дивергенции, практически не зависит от соотношения численностей р. Изменение величины р лишь сдвигает границы бассейнов притяжения (сепаратрис) следующим образом. Если р > 1/2 (первая популяция более многочисленна), то граница смещается в сторону мономорфизма по генотипу А, то есть к точке Е2 (рис. 8, a). Если р < 1/2 (первая популяция малочисленна), то в сторону мономорфизма по генотипу а, то есть к точке Ео (рис. 8, b).
Заключение
Итак, в данной работе мы еще раз рассмотрели простейшие модели первичной генетической дивергенции и исследовали их с помощью современных методов нелинейного анализа и компьютерной визуализации полученных результатов. Были построены параметрические портреты рассматриваемых систем, фазовые портреты, совмещенные с бассейнами притяжения устойчивых стационарных состояний, а также интересные, с нашей точки зрения, бифуркационные диаграммы и примеры динамики. Все это позволило существенно расширить первоначальный анализ, проведенный в ранних работах по генетической дивергенции, сделать его полнее и значительно нагляднее.
Что же касается содержательных результатов, то основной вывод связан здесь с необходимостью ограничения численности популяции: только для систем популяций, численности которых как-то ограничены (например, соотношения численностей не меняется во времени), оказывается возможна генетическая дивергенция. Это приводит к естественному предположению, что в возникновении генетической дивергенции кроме популяционно-генетических факторов существенную роль играет и экологическое лимитирование роста численности. Прямое введение экологического лимитирования в популяционно-генетические модели дивергенции и подробный современный анализ этих эколого-генетических моделей сулят заманчивые перспективы будущих исследований.
Список литературы
1. Любич Ю. И. Основные понятия и теоремы эволюционной генетики свободных популяций // УМН. 1971. Т. 26, № 5(161). С. 51-116.
2. Базыкин А. Д. Пониженная приспособленность гетерозигот в системе двух смежных популяций // Генетика. 1972. Т. 8, № 11. С. 155-161.
3. Базыкин А. Д. Отбор и генетическая дивергенция в системах локальных популяций и популяциях с непрерывным ареалом (математическая модель) // Проблемы эволюции. 1973. Т. 3. С. 231-241.
4. Фрисман Е. Я., Шапиро А. П.Избранные математические модели дивергентной эволюции популяций. М.: Наука, 1977. 152 с.
5. Фрисман Е. Я. Первичная генетическая дивергенция (Теоретический анализ и моделирование). Владивосток: ДВНЦ АН СССР, 1986. 160 с.
6. Арнольд В. И., Афраймович В. С., Ильяшенко Ю. С., Шильников Л. П. Теория бифуркаций (Динамические системы - 5) // Итоги науки и техн. Сер. Соврем. пробл. мат. Фундам. направления. Т. 5. М.: ВИНИТИ, 1986. С. 5-218.
7. Шильников Л. П., Шильников А. Л., Тураев Д. В., Чуа Л.Методы качественной теории в нелинейной динамике. Москва-Ижевск: Институт компьютерных исследований, 2003. 428 с.
8. DhoogeA., Govaerts W., Kuznetsov Y.A., MeijerH. G. E., Sautois B. New features of the software MatCont for bifurcation analysis of dynamical systems // Mathematical and Computer Modelling of Dynamical Systems. 2008. Vol. 14, no. 2. P. 147-175. DOI: 10.1080/13873950701742754.
9. Stewart I., Elmhirst T., Cohen /.Symmetry-Breaking as an Origin of Species // In: Buescu J., Castro S. B. S. D., da Silva Dias A. P., Labouriau I. S. (eds) Bifurcation, Symmetry and Patterns. Trends in Mathematics. Basel: Birkhauser, 2003. P. 3-54. DOI: 10.1007/978-3-0348-7982-8_1.
10. Btirger R. A survey of migration-selection models in population genetics // Discrete & Continuous Dynamical Systems - B. 2014. Vol. 19, no. 4. P. 883-959. DOI: 10.3934/dcdsb.2014.19.883.
11. Yamamichi M., Ellner S. P. Antagonistic coevolution between quantitative and Mendelian traits // Proc. R. Soc. B. 2016. Vol. 283, no. 1827. P. 20152926. DOI: 10.1098/rspb.2015.2926.
12. Telschow A., Hammerstein P., Werren J. H. The effect of Wolbachia on genetic divergence between populations: Models with two-way migration // The American Naturalist. 2002. Vol. 160, no. S4. P. S54-S66. DOI: 10.1086/342153.
13. Жданова О. Л., Фрисман Е.Я. Динамические режимы в модели однолокусного плотност-но-зависимого отбора//Генетика. 2005. Т. 41, № 11. С. 1575-1884.
14. Altrock P. M., Traulsen A., Reeves R. G., Reed F. A. Using underdominance to bi-stably transform local populations // J. Theor. Biol. 2010. Vol. 267, no. 1. P. 62-75.
DOI: 10.1016/j.jtbi.2010.08.004.
15. Yeaman S., Otto S.P. Establishment and maintenance of adaptive genetic divergence under migration, selection, and drift // Evolution. 2011. Vol. 65, no. 7. P. 2123-2129.
DOI: 10.1111/j.1558-5646.2011.01277.x.
16. Laruson A. J., Reed F.A. Stability of underdominant genetic polymorphisms in population networks // J. Theor. Biol. 2016. Vol. 390. P. 156-163. DOI: 10.1016/j.jtbi.2015.11.023.
17. Yamamichi M., Hoso M.Roles of maternal effects in maintaining genetic variation: Maternal storage effect//Evolution. 2017. Vol. 71, no. 2. P. 449-457. DOI: 10.1111/evo.13118.
18. Wakeley J. The effects of subdivision on the genetic divergence of populations and species // Evolution. 2000. Vol. 54, no. 4. P. 1092-1101. DOI: 10.1111/j.0014-3820.2000.tb00545.x.
19. Фрисман Е.Я., Жданова О. Л., Кулаков М.П., Неверова Г. П., Ревуцкая О. Л. Математическое моделирование популяционной динамики на основе рекуррентных уравнений: результаты и перспективы. Ч. II // Известия РАН. Серия биологическая. 2021. № 3. С. 227-240. DOI: 10.31857/S000233292103005X.
20. Neverova G. P., Zhdanova O. L., Frisman E. Y. Effects of natural selection by fertility on the evolution of the dynamic modes of population number: bistability and multistability // Nonlinear Dyn. 2020. Vol. 101, no. 1. P. 687-709. DOI: 10.1007/s11071-020-05745-w.
21. Shea K., Metaxas A., Young C.R., Fisher C.R. Processes and Interactions in Macrofaunal Assemblages at Hydrothermal Vents: A Modeling Perspective // In: Lowell R. P., Seewald J. S., Metaxas A., Perfit M. R. (eds) Magma to Microbe: Modeling Hydrothermal Processes at Ocean Spreading Centers. Vol. 178 of Geophysical Monograph Series. Washington, DC: Blackwell Publishing Ltd, 2008. P. 259-274. DOI: 10.1029/178GM13.
22. Sundqvist L., Keenan K., Zackrisson M., Prodohl P., Kleinhans D. Directional genetic differentiation and relative migration // Ecol. Evol. 2016. Vol. 6, no. 11. P. 3461-3475.
DOI: 10.1002/ece3.2096.
References
1. Lyubich YI. Basic concepts and theorems of the evolutionary genetics of free populations. Russian Math. Surveys. 1971;26(5):51-123. DOI: 10.1070/RM1971v026n05ABEH003829.
2. Bazykin AD. Disadvantages of heterozygotes in a system of two adjacent populations. Sov. Genet. 1974;8(11):1453-1457.
3. Bazykin AD. Selection and genetic divergence in local systems populations and populations with a continuous areal (mathematical model). Evolution Problems. 1973;3:231-241 (in Russian).
4. Frisman EY, Shapiro AP. Selected Mathematical Models of Divergent Evolution of Populations. Moscow: Nauka; 1977. 152 p. (in Russian).
5. Frisman EY. Primary Genetic Divergence (Theoretical Analysis and Modeling). Vladivostok: Far East Scientific Center of the Academy of Sciences of the USSR; 1986. 160 p. (in Russian).
6. Arnold VI, Afraimovich VS, Ilyashenko YS, Shilnikov LP. Bifurcation Theory (Dynamical Systems - 5). In: Results of Science and Technology. Series Modern Problems of Mathematics. Fundamental Directions. Vol. 5. Moscow: VINITI; 1986. P. 5-218 (in Russian).
7. Shilnikov LP, Shilnikov AL, Turaev DV, Chua LO. Methods of Qualitative Theory in Nonlinear Dynamics. Part I. Singapore: World Scientific; 1998. 416 p. DOI: 10.1142/3707.
8. Dhooge A, Govaerts W, Kuznetsov YA, Meijer HGE, Sautois B. New features of the software MatCont for bifurcation analysis of dynamical systems. Mathematical and Computer Modelling of Dynamical Systems. 2008;14(2):147-175. DOI: 10.1080/13873950701742754.
9. Stewart I, Elmhirst T, Cohen J. Symmetry-Breaking as an Origin of Species. In: Buescu J, Castro SBSD, da Silva Dias AP, Labouriau IS, editors. Bifurcation, Symmetry and Patterns. Trends in Mathematics. Basel: Birkhauser; 2003. P. 3-54. DOI: 10.1007/978-3-0348-7982-8_1.
10. Burger R. A survey of migration-selection models in population genetics. Discrete & Continuous Dynamical Systems - B. 2014;19(4):883-959. DOI: 10.3934/dcdsb.2014.19.883.
11. Yamamichi M, Ellner SP. Antagonistic coevolution between quantitative and Mendelian traits. Proc. R. Soc. B. 2016;283(1827):20152926. DOI: 10.1098/rspb.2015.2926.
12. Telschow A, Hammerstein P, Werren JH. The effect of Wolbachia on genetic divergence between populations: Models with two-way migration. The American Naturalist. 2002;160(S4):S54-S66. DOI: 10.1086/342153.
13. Zhdanova OL, Frisman EY. Dynamic regimes in a model of single-locus density-dependent selection. Russ. J. Genet. 2005;41(11):1302-1310. DOI: 10.1007/s11177-005-0233-3.
14. Altrock PM, Traulsen A, Reeves RG, Reed FA. Using underdominance to bi-stably transform local populations. J. Theor. Biol. 2010;267(1):62-75. DOI: 10.1016/j.jtbi.2010.08.004.
15. Yeaman S, Otto SP. Establishment and maintenance of adaptive genetic divergence under migration, selection, and drift. Evolution. 2011;65(7):2123-2129.
DOI: 10.1111/j.1558-5646.2011.01277.x.
16. Laruson Aj, Reed FA. Stability of underdominant genetic polymorphisms in population networks. J. Theor. Biol. 2016;390:156-163. DOI: 10.1016/j.jtbi.2015.11.023.
17. Yamamichi M, Hoso M. Roles of maternal effects in maintaining genetic variation: Maternal storage effect. Evolution. 2016;71(2):449-457. DOI: 10.1111/evo.13118.
18. Wakeley J. The effects of subdivision on the genetic divergence of populations and species. Evolution. 2000;54(4):1092-1101. DOI: 10.1111/j.0014-3820.2000.tb00545.x.
19. Frisman EY, Zhdanova OL, Kulakov MP, Neverova GP, Revutskaya OL. Mathematical modeling of population dynamics based on recurrent equations: Results and prospects. Part II. Biol. Bull. Russ. Acad. Sci. 2021;48(3):227-240. DOI: 10.1134/S1062359021030055.
20. Neverova GP, Zhdanova OL, Frisman EY. Effects of natural selection by fertility on the evolution of the dynamic modes of population number: bistability and multistability. Nonlinear Dyn. 2020;101(1):687-709. DOI: 10.1007/s11071-020-05745-w.
21. Shea K, Metaxas A, Young CR, Fisher CR. Processes and Interactions in Macrofaunal Assemblages at Hydrothermal Vents: A Modeling Perspective. In: Lowell RP, Seewald JS, Metaxas A, Perfit MR, editors. Magma to Microbe: Modeling Hydrothermal Processes at Ocean Spreading Centers. Vol. 178 of Geophysical Monograph Series. Washington, DC: Blackwell Publishing Ltd; 2008. P. 259-274. DOI: 10.1029/178GM13.
22. Sundqvist L, Keenan K, Zackrisson M, Prodohl P, Kleinhans D. Directional genetic differentiation and relative migration. Ecol. Evol. 2016;6(11):3461-3475. DOI: 10.1002/ece3.2096.
Фрисман Ефим Яковлевич - родился в Сталинабаде (ныне Душанбе, 1948). Окончил Новосибирский государственный университет (1971). Защитил диссертации на соискание ученой степени кандидата биологических наук по специальности «генетика» (1982) и доктора биологических наук по специальности «биофизика» (1989). В 2011 году избран членом-корреспондентом Российской академии наук по специальности «общая биология». После окончания университета работал в Институте автоматики и процессов управления ДВО РАН: занимал должности от стажера-исследователя до заведующего лабораторией. С 2002 по 2018 - директор Института комплексного анализа региональных проблем ДВО РАН. Ныне - научный руководитель этого института. Автор более 200 публикаций по общим вопросам математического моделирования и по конкретным моделям экологических и по-пуляционных систем. Научные интересы связаны с моделированием динамики популя-ционных и экологических систем, математической популяционной генетикой и задачами оптимального управления.
Россия, 679016 Биробиджан, ул. Шолом-Алейхема, 4
Институт комплексного анализа региональных проблем ДВО РАН
E-mail: [email protected]
ORCID: 0000-0003-1629-2610
AuthorID: 1986
Кулаков Матвей Павлович - родился в Биробиджане (1982). Окончил Биробиджанский государственный педагогический институт (2004). Защитил диссертацию на соискание ученой степени кандидата физико-математических наук по специальности «биофизика» (2018). Старший научный сотрудник лаборатории математического моделирования популя-ционных и экологических систем Института комплексного анализа региональных проблем ДВО РАН. Научные интересы связаны с использованием методов нелинейной динамики, качественной теории дифференциальных уравнений, бифуркаций, хаоса и синхронизации для изучения особенностей функционирования пространственно распределенных живых систем: популяций животных и биологических сообществ.
Россия, 679016 Биробиджан, ул. Шолом-Алейхема, 4
Институт комплексного анализа региональных проблем ДВО РАН
E-mail: [email protected]
ORCID: 0000-0002-7060-2731
AuthorID: 170285
