Журнал общей биологии, 2020, T. 81, № 4, стр. 257-271

МАТРИЧНАЯ МОДЕЛЬ ЦЕНОПОПУЛЯЦИИ A. ALBANA

Объект. По нашим данным, A. albana в условиях альпийских лишайниковых пустошей проявляет себя как травянистый хамефит, летнезеленый стержнекорневой многолетний монокарпик (Казанцева, 2016; Логофет и др., 2017б). Вид включен в “Красную книгу Краснодарского края” (2007) и “Красную книгу Республики Адыгея” (2012) как редкий категории 3. Биология, экология и онтогенез вида описаны ранее (Казанцева, 2016; Логофет и др., 2017б), и построен граф жизненного цикла по стадиям онтогенеза (рис. 1), которые различимы в поле по морфологии растения. Наличие в ГЖЦ нетривиальных задержек в стадиях im и v, а также ускоренного перехода pl im есть проявление поливариантности онтогенеза A. albana в условиях альпийского пояса Северо-Западного Кавказа (Логофет и др., 2019).

Методы исследования. Исследования проводились в Карачаево-Черкесской Республике, на территории Тебердинского государственного природного биосферного заповедника, на постоянных площадках, заложенных в 2009 г. на лишайниковой пустоши, на горе Малая Хатипара, 2800 м над ур. м. Метод ежегодной переписи всех растений на площадках описан в работе Д.О. Логофета с соавт. (2017б). Мониторинг проводили каждый год в конце августа в течение 11 лет с 2009 по 2019 год включительно.

Собранные данные относятся к типу “идентифицированные особи” (“identified individuals”; Caswell, 2001, p. 134), т.е. содержат сведения о стадии каждого растения A. albana с момента его появления на площадке и до момента цветения с шагом в один год. Поэтому, наряду со структурой ценопопуляции в каждый год наблюдения, т.е. вектором(-столбцом)

компоненты которого суть числа особей, обнаруженных в соответствующих стадиях онтогенеза в год наблюдения t (табл. 1), известны и события переходов между стадиями за прошедший год, а также размеры ежегодных пополнений.

Матричная модель и калибровка матриц. Наблюдаемые изменения вектора x(t) по годам описываются основным уравнением модели в векторно-матричной форме:

где L(t) – годичная ПМП, ассоциированная с ГЖЦ, который представлен на рис. 1. Все 10 ПМП L(t) имеют размер 5 × 5 и общее строение: хотя в отдельные годы некоторые из демографических параметров (vital rates) a, b, …, l, m, оказываются нулями из-за отсутствия соответствующих переходов в данных (табл. 2). Эти данные вкупе с меняющейся структурой популяции позволяют калибровать годичные ПМП абсолютно точно согласно уравнению (1) и смыслу параметров как частот соответствующих событий переходов, зафиксированных в наблюдениях (табл. 3).

Прогноз выживания ценопопуляции дается величиной λ₁(L) – доминантным собственным числом калиброванной ПМП L, которое показывает асимптотическую скорость роста популяции: если L(t) = L остается постоянной во времени, то при t → ∞ структура популяции x(t) сходится к равновесию – собственному вектору x*, отвечающему λ₁(L), – а ее динамика x(t) ~ λ₁(L)^tx*, т.е. стремится к геометрическому росту (или убыли) с показателем λ₁(L) > 0 (Logofet, 1993; Caswell, 2001; Логофет, Уланова, 2018). Поэтому прогноз положительный, если λ₁(L) > 1, и отрицательный, если λ₁(L) < 1. В неавтономной (т.е. с зависимостью параметров от времени) модели вместе с изменением матрицы L(t) в зависимости от t существенно меняется и величина λ₁(L), оставляя результирующий прогноз неопределенным (Логофет и др., 2019, табл. 3).

Прогноз выживания через осреднение годичных ПМП. Чтобы зафиксировать онтогенетические переходы и пополнение популяции, необходимы данные двух последовательных переписей (в моменты t и t + 1), и потому число калиброванных ПМП всегда на единицу меньше числа лет наблюдений (11 лет наблюдений обеспечивают калибровку 10 матриц L(t)). Нет оснований полагать, что они окажутся одинаковыми, а практика калибровки показывает различия не только количественные, но и качественные: в отдельные годы прогноз выживания меняется на противоположный (Логофет и др., 2017а, б, 2019). Итог всего периода мониторинга способна отразить матрица G, найденная методом СМО неотрицательных матриц (Logofet, 2018; Логофет и др., 2019).

Логика СМО такова: согласно калиброванной модели (1) вектор, x(2009), структуры популяции в начальный момент наблюдений t = 2009 преобразуется в x(2019) – вектор структуры в последний момент t = 2019 – путем последовательного умножения на 10 годичных ПМП L(t), т.е. (после переобозначения L(t) = L_j , j = t – 2009)

То же самое должна делать и средняя матрица G, т.е.

Уравнение (4) указывает на мультипликативный характер осреднения, а общее для всех усредняемых матриц строение (2) требует того же и от средней матрицы G. Логично потребовать, кроме того, чтобы все элементы средней матрицы не выходили за границы соответствующих элементов усредняемых матриц.

Данная логика приводит к математической задаче минимизации ошибки приближения в решении матричного уравнения (4) (Логофет и др., 2019, Приложение А). Положение доминантного собственного числа λ₁(G) средней матрицы G относительно 1 позволяет сделать вывод о судьбе ценопопуляции в долговременной перспективе на основании нескольких лет наблюдений: рост, если λ₁(G) > 1; убыль, если λ₁(G) < 1; стабилизация, если λ₁(G) = 1. Десять лет мониторинга ценопопуляции A. albana привели к значению λ₁(G) = = 0.89262 < 1 у средней матрицы G, найденной с погрешностью в сотых долях (Логофет и др., 2019, табл. 3), т.е. к негативному прогнозу выживания ценопопуляции. Отметим, что предельный вектор, y*, матрицы G был весьма близок к вектору x(2019) – относительной структуре популяции в последний год наблюдения: мера несходства Кейфица (Keyfitz, 1968; Логофет и др., 2019) между ними составляла менее 1%. Итог 11 лет мониторинга, т.е. результат СМО 10 годичных ПМП, представлен в настоящей статье.

МЕТОД ОБРАТНОГО ПРОГНОЗА В МАТРИЧНОЙ МОДЕЛИ A. ALBANA

Обратный прогноз призван подтвердить или поставить под сомнение выводы прогноза прямого – в частности, по значению λ₁(G_–) структурно-мультипликативного среднего, G_–, годичных матриц обратного прогноза и предельному вектору, $z_{{{{G}_{ - }}}}^{*}$, средней матрицы. Как отмечено выше, рутинные методы обратного прогноза работают лишь в специальных случаях, а для матричной модели формальное обращение уравнения (1), т.е.

неприемлемо, поскольку матрица L(t)^–1, обратная к ПМП, не обязана даже быть неотрицательной, не говоря уже о ее строении, никак не отражающем ГЖЦ (рис. 1).

Простой альтернативный подход к обратному прогнозу: получать наблюдавшиеся структуры ценопопуляции в обратном порядке, от 2019 к 2009 году, с помощью ПМП, B(t), отвечающих прежнему ГЖЦ (рис. 1), – наталкивается на неоднозначность решения задачи. В прямой модели все переходы между стадиями, состоявшиеся за один год, т.е. все изменения вектора x(t) = [p(t), j(t), im(t), v(t), g(t)]^T во времени, известны из наблюдений, а обратных наблюдений в реальности нет. Основное уравнение обратной модели, которое связывает данные наблюдений в обратном порядке,

имеет при каждом t много разных решений даже в рациональных числах (лишь 5 скалярных уравнений для 10 неизвестных элементов матрицы B). Поэтому нужен специальный подход к реализации идеи обратного прогноза.

Идея обратного прогноза через обращение вспять оси времени в наблюдениях была успешно реализована для дискретной марковской цепи (Logofet, 2019a), и в следующем разделе мы применяем ее для построения обратного прогноза в матричной модели популяции.

Обратный прогноз через “обращение” годичной ПМП. Как отмечено выше, обращение годичной ПМП в смысле вычисления обратной матрицы неприемлемо для обратного прогноза, и потому “обращение” приобретает в нашем случае специальный смысл. Данные ежегодных наблюдений содержат сведения не только о структуре ценопопуляции (табл. 1), но и обо всех переходах по стадиям онтогенеза, состоявшихся за прошедший год (табл. 2) согласно ГЖЦ (рис. 1). Представим эти данные в форме целочисленных таблиц T_L(t), имеющих такое же строение, как и сами годичные ПМП, т.е. в соответствии с ГЖЦ (табл. 2), но показывавших не частоты или коэффициенты, а фактические числа переходов и пополнений. Тогда строчные суммы элементов T_L(t) будут закономерно совпадать с размерами соответствующих стадийных групп (технические подробности в Приложении А).

Обратить время в этих данных значит изменить направления всех переходов жизненного цикла на противоположные или сделать то же самое с дугами ГЖЦ (рис. 2). При этом само количество переходов не изменяется, и в терминах таблицы T_L(t) обращение означает транспонирование (т.е. перенос каждого элемента симметрично относительно главной диагонали) или, что эквивалентно, замену строк таблицы ее столбцами. Полученную таким образом таблицу T_L(t – 1)^T нужно превратить в матрицу B(t) обратного прогноза так, чтобы выполнялось условие калибровки (6). Для этого достаточно положить

и найти положительные калибрующие множители d₁, d₂, …, d₅ из уравнения (6) обратного прогноза: (технические подробности в Приложении Б).

Таким образом, каждая из 10 матриц (7) обратного прогноза неотрицательна, и ее строение соответствует сильно связному орграфу (рис. 2), т.е. матрица неразложима (irreducible; Horn, Johnson, 1990), и по теореме Перрона–Фробениуса (Гантмахер, 1967) она имеет доминантное собственное число, λ₁(B(t)) > 0, и соответствующий ему положительный собственный вектор, y* > 0, который представляет собой устойчивое предельное распределение (Caswell, 2001; Логофет, Уланова, 2018) в модели (6) обратного прогноза. После калибровки матрицы B(t) все ее собственные числа и соответствующие им собственные векторы вычисляем стандартной процедурой, из результатов которой выбираем доминантное собственное число λ₁ и нормируем соответствующий ему вектор (Приложение Б). По смыслу обратного прогноза ожидается, что λ₁(B(t)) > 1 для тех лет, для которых λ₁(L(t – 1)) < 1, т.е. результирующий обратный прогноз следует получать, как и прямой прогноз, через осреднение годичных матриц (Приложение В).

Обратный прогноз через осреднение годичных матриц. Логика осреднения годичных матриц обратного прогноза такая же, как и в прямой модели (Логофет и др., 2019): структура популяции, зафиксированная в первый год наблюдения, должна переходить в структуру последнего года путем ее умножения на последовательность годичных матриц в хронологическом порядке (от 2019 к 2010 году в данном случае), и такой же результат должна давать средняя матрица G_– в степени, равной числу шагов. Формально

т.е. структурно-мультипликативную среднюю матрицу G_– следовало бы вычислять как матричный корень 9-й степени из произведения 9 годичных матриц. Поскольку такое произведение представляет собой положительную матрицу (все элементы больше 0), матричный корень в лучшем случае будет также положительным (McDonald et al., 2014; Politi, Popolizio, 2015). Но у средней матрицы строение должно быть таким же, как и у матриц усредняемых, т.е. а значения ее элементов не должны выходить из границ, в которых находятся элементы усредняемых матриц.

Условие (10) приводит к переопределенности (GUFO.ME, 2020) матричного уравнения (9) как системы из 25 скалярных уравнений для 10 неизвестных a_–, b_–, …, l_–, m_–, т.е. к отсутствию точного решения. Таким образом, поиск средней матрицы G_– сводится к приближенному решению уравнения (9), т.е. к задаче ограниченной минимизации ошибки приближения (Приложение В).

РЕЗУЛЬТАТЫ

Калибровка годичных ПМП. Результаты калибровки матриц L(t) представлены в табл. 3. К девяти матрицам, калиброванным ранее (Логофет и др., 2019, табл. 3), последний, 2019-й, год наблюдений добавил матрицу L(2018) с доминантным собственным числом, меньшим 1 (последняя строка табл. 3).

Калиброванные матрицы обратного прогноза. Результаты вычисления годичных матриц B(t) по формулам (7), (8) представлены в табл. 4. Вариации меры приспособленности локальной популяции (доминантного собственного числа λ₁(B(t)) по годам весьма значительны, как и в прямой модели. Как и ожидалось, λ₁(B(t)) > 1 во всех тех случаях, когда прямая модель давала значение λ₁(L(t – 1)) < 1 (табл. 3), и наоборот. Явной закономерности в предельных распределениях (векторах y*(t)) не наблюдается: в разные годы доминирование перемещается между стадиями j, im и v, и только при наименьшем значении λ₁(t = 2011) доминирующей оказалась стадия p.

Осреднение годичных матриц прямого и обратного прогноза. Структурно-мультипликативное среднее 10 годичных ПРМ L(t), полученное как приближенное решение уравнения (4), показано в первой строке табл. 5. Добавление 10-й матрицы с λ₁(L(2018)) < 1 (табл. 3) к набору из 9 усредняемых ПМП (Логофет и др., 2019, табл. 3) закономерно снизило значение λ₁(G) (в сотых долях), т.е. 10-летний мониторинг дает прогноз с более быстрым убыванием ценопопуляции, нежели 9‑летний. Ошибка приближения снизилась почти на порядок. Предельное распределение стадийных групп сохранило двугорбый характер с доминированием в стадиях j и v, и оно по-прежнему мало отличается от распределения x(2019) в последний год наблюдения (4.5%, см. последний столбец табл. 5).

Приближенное решение уравнения (9) для 10 матриц из табл. 4, т.е. СМО годичных матриц обратного прогноза, с приемлемой точностью приближения дает матрица G_–, представленная в последней строке табл. 5. Строение G_–, очевидно, соответствует (10), т.е. обращенному ГЖЦ (рис. 2), а значения ее элементов не выходят из границ, заданных соответствующими элементами усредняемых матриц (Приложение В, условие (В.2)). Доминантное собственное число λ₁(G_–) матрицы обратного прогноза оказалось больше 1, подтвердив ожидания, а соответствующий доминантный собственный вектор z*, т.е. предельное распределение численности стадийных групп при движении назад, оказался весьма близким к фактическому в начальный год наблюдений: показатель несходства между ними по Кейфицу около 2% (табл. 5).

ОБСУЖДЕНИЕ И ЗАКЛЮЧЕНИЕ

Ежегодные наблюдения за ценопопуляцией проломника Androsace albana на постоянных площадках продолжаются с 2009 г. (Логофет и др., 2017б, 2019), и каждый очередной год добавляет новую матрицу в набор годичных ПМП, расширяя спектр вариаций условий среды, опосредованно отраженных в матрицах. Этот набор служит исходным пунктом в оценке жизнеспособности ценопопуляции по каждому из двух дальнейших альтернативных направлений. Оценка стохастической скорости роста, λ_S, популяции в случайно меняющейся среде (Tuljapurkar, 1986, 1990; Caswell, 2001) преобладает в англоязычной литературе (Caswell, 2010), а вычисление доминантного собственного числа λ₁(G) средней матрицы G следует оригинальной концепции структурно-геометрического (структурно-мультипликативного в данной статье) осреднения неотрицательных матриц (Logofet, 2018, 2019b; Логофет и др., 2019). Оба подхода дают качественно одинаковый – отрицательный – прогноз выживания ценопопуляции A. albana (Логофет и др., 2019), хотя количественные результаты, т.е. скорости ежегодного снижения численности, заметно различаются (Логофет и др., 2019, табл. 3, 4).

Поскольку результат СМО годичных ПМП представляет собой приближенное решение уравнения (9) (Приложение В), причину вышеназванного различия можно усмотреть в неточности приближения, так как λ₁ любой неотрицательной неразложимой матрицы весьма чувствительно к вариациям ее элементов (монотонно возрастает при увеличении любого элемента и наоборот; Гантмахер, 1967, §XIII.3, Теорема 6; Логофет, Уланова, 2018, утверждение (П7.1–2)). Это классическое положение имеет качественный характер и нуждается в дальнейшем теоретическом уточнении в связи с задачей СМО годичных ПМП. Несколько математических результатов были получены в этой связи (Protasov, 2019; Протасов, Широков, 2019), однако до практических рекомендаций дело еще не дошло.

В подобной ситуации возрастает роль альтернативных подходов к задаче прогноза, способных подтвердить выводы, полученные методом СМО, и обратный прогноз выступает одним из таких методов. Годичные матрицы, B(t), обратного прогноза (табл. 4) обнаружили не меньшее разнообразие характеристик, чем ПМП прямой модели (табл. 3), в зависимости от года наблюдения t. Примечательно, что мера, λ₁(B(t)), “обратной” жизнеспособности ценопопуляции оказалась больше 1 во все те годы, когда “прямая” мера на предыдущем шаге, λ₁(L(t – 1), была меньше 1, и наоборот, – подтверждая формально очевидное положение здравого смысла.

Как и в прямой модели, итоговый прогноз жизнеспособности дает структурно-мультипликативная средняя матрица, G_–, ряда годичных матриц своим доминантным собственным числом λ₁(G_–). Оно закономерно оказалось больше 1 (табл. 5), указывая на рост ценопопуляции в обратном направлении времени, а предельное распределение, z*, размеров стадийных групп практически совпадает с финальным распределением, x(2009), при движении назад (различие менее 2%; табл. 5).

Предложенная методика обратного прогноза оказалась эффективной в матричной модели популяции A. albana со стадийной структурой, и она не специфична для A. albana: единственное принципиальное условие – неразложимость годичных матриц, т.е. сильная связность исходного ГЖЦ (рис. 1), – обеспечено непрерывностью онтогенеза, начальная стадия которого пополняется в результате репродукции. Поэтому методика обратного прогноза применима во всех случаях, где имеется временнóй ряд данных типа “идентифицированные особи” (“identified individuals”; Caswell, 2001, p. 134). В приложении к данным многолетнего мониторинга ценопопуляции A. albana данная методика подтвердила справедливость “неутешительного прогноза выживания ценопопуляции” на альпийской пустоши (Логофет и др., 2019).

Сбор и обработка полевого материала осуществлены при финансовой поддержке РФФИ (грант № 14-04-00214). Модельная часть исследований поддержана РФФИ, грант № 19-04-01227. Вычисления проделаны в среде MATLAB, версия R2019b.