CC BY
113
Прикладные задачи
нелинейной теории колеоании и волн
УДК 517.9, 621.372 https://doi.org/10.18500/0869-6632-2019-27-2-5-20
SIRS-модель распространения инфекций с динамическим регулированием численности популяции: Исследование методом вероятностных клеточных автоматов
А.В. Шабунин
Саратовский государственный университет им. Н.Г. Чернышевского Россия, 410012 Саратов, Астраханская, 83 E-mail: shabuninav@info.sgu.ru Поступила в редакцию 14.12.2018, принята к публикации 26.02.2019
Цель. Построение модели распространения инфекции в виде решетки стохастических клеточных автоматов, в которой могут существовать нетривиальные динамические режимы; исследование ее динамики и сопоставление с моделью среднего поля. Метод. численное моделирования квадратной решетки клеточных автоматов по методу Монте-Карло, теоретическое и численное исследование устройства фазового пространства системы обыкновенных дифференциальных уравнений модели среднего поля. Результаты. Построена модифицированная SIRS-модель распространения эпидемий в виде решетки стохастических клеточных автоматов. В модели используется динамическое регулирование численности населения с ограничением максимального числа особей популяции и влиянием заболевания на процессы воспроизводства. Обнаружено, что в зависимости от управляющих параметров она демонстрирует четыре разных установившихся режима: (а) вымирание популяции, (б) стационарный ход заболевания, (в) полное излечение популяции и (г) самоподдерживающиеся колебания числа инфицированных, сопровождающиеся колебаниями общей численности популяции. Последний режим проявляется вблизи границы зоны полного излечения и характеризуется нерегулярными колебаниями числа заболевших с выраженной периодической составляющей. Показано, что при периодическом изменении параметров модель демонстрирует зашумленные периодические или квазипериодические колебания. Обсуждение. Поведение модели клеточных автоматов в целом соответствует модели среднего поля, однако имеются как количественные, так и качественные расхождения. Количественные расхождения заключаются в небольшом смещении значений средней концентрации заболевших особей относительно теоретически предсказанных величин. Расхождения существенно уменьшаются, если в популяции присутствует сильная миграция особей, что можно объяснить большей однородностью популяции в этом случае. В отличие от уравнений среднего поля модель клеточных автоматов демонстрирует необходимость отличного от нуля порога заболеваемости для поддержания эпидемии, а также наличие колебательного режима. Предполагаемая причина колебаний кроется в вероятностном характере работы клеточных автоматов, что приводит к случайным флуктуациям, которые могут усиливаться, подобно тому, как это происходит в возбудимых системах, находящихся под действием внешнего источника шума.
Ключевые слова: SIRS-модель, решетки клеточных автоматов, популяционная динамика.
Образец цитирования: Шабунин А.В. SIRS-модель распространения инфекций с динамическим регулированием численности популяции: Исследование методом вероятностных клеточных автоматов // Изв. вузов. ПНД. 2019. T. 27, № 2. С. 5-20. https://doi.org/10.18500/0869-6632-2019-27-2-5-20
https://doi.org/10.18500/0869-6632-2019-27-2-5-20
SIRS-model with dynamic regulation of the population: Probabilistic cellular automata approach
A.V. Shabunin
Saratov State University 83, Astrakhanskaya str., 410012 Saratov, Russia E-mail: shabuninav@info.sgu.ru Received 14.12.2018, accepted for publication 26.02.2019
Aim. Construction a model of infection spread in the form of a lattice of stochastic cellular automata which can demonstrate nontrivial oscillating regimes; investigation of its dynamics and comparison with the mean-field model. Method. Numerical simulation of the square lattice of cellular automata by the Monte Carlo approach, theoretical and numerical study of the structure of the phase space of its mean-field model. Results. A modified SIRS-model of epidemic propagation has been proposed in the form of a lattice of stochastic cellular automata. It uses dynamic regulation of the population with limit on the maximum number of individuals and takes into account the influence of the disease on reproduction processes. The study has discovered that, depending on the control parameters, the model demonstrates four different steady-state regimes: (a) total dying of the population, (b) stationary course of the disease, (c) complete recovery of the deases and (d) self-sustaining fluctuations of the infected, which accompanied by fluctuations of the population. The last regime is manifested near the border of the zone of complete recovery and is characterized by irregular bursts in the number of ill with a clearly distinguishable periodic component. At periodic modulation of the parameters the model demonstrates noisy periodic or quasi-periodic oscillations. Discussion. In general, the cellular automata model behaves similarly to the mean field equations, however there are quantitative and qualitative differences.The first ones consist in small shifts in the average concentrations compared to theoretically predicted values. These shifts are significantly reduced when individuals in the population can migrate. This effect can be explained by the greater homogeneity of populations in this case. The qualitative difference is that, the model of cellular automata demonstrates the effects, which are absent in the mean-field equations: a threshold of ill individuals necessary for the epidemic spread, and the oscillatory regime. A suggested cause of the last effect is the probabilistic nature of cellular automata, which lead to random fluctuations of the concentrations. The last ones can be amplified just as it happens in excitable systems under external noise.
Key words: SIRS-model, lattice of cellular automata, population dynamics.
Reference: Shabunin A.V. SIRS-model with dynamic regulation of the population: Probabilistic cellular automata approach. Izvestiya VUZ, Applied Nonlinear Dynamics, 2019, vol. 27, no. 2, pp. 5-20. https://doi.org/10.18500/0869-6632-2019-27-2-5-20
Введение
Методы нелинейной динамики позволяют анализировать многие явления реальной жизни путем построения простых (так называемых качественных) моделей. Эти модели позволяют понять характер наблюдаемых явлений и дают возможность управлять их поведением, то есть показывают, за счет каких параметров или воздействий можно получить желаемое изменение тех или иных характеристик рассматриваемой системы. Естественной областью применения такого подхода является исследование распространения эпидемий инфекционных заболеваний в человеческих сообществах или природных популяциях [1-3].
Математическое моделирование эпидемий является одним из разделов популяционной динамики, используемой как в биологии, так и в других науках, например, социологии. Классические модели распространения инфекций представляют собой системы обыкновенных дифференциальных уравнений (ОДУ). Наиболее известной из них является модель SIR, предложенная в 1920-х годах Кермаком и МакКендриком [4], а также ее расширения SIRS-модель. В SIRS-модели популяцию разбивают на группы здоровых и восприимчивых (S - Susceptible), больных (I - Infectious) и выздоровевших (R - Recovered) особей и строят системы ОДУ, определяющие закон изменения относительного числа особей в каждой из групп, опираясь на предположение о случайном и равномерном распределении особей в популяции. Такую систему называют
моделью среднего поля. Преимуществом моделей среднего поля является возможность их аналитического исследования. Однако сделанное предположение ограничивает их применимость, поскольку в реальных системах наблюдается, напротив, сильная неоднородность в распределении инфицированных особей (очаговый характер заболевания) [2].
Альтернативой ОДУ является моделирование методами решеток клеточных автоматов [5-7]. Клеточным автоматом (КА) называется система, обладающая конечным числом состояний, переключения между которыми происходят дискретно во времени по заданному закону. Если закон изменения состояний клетки представляет собой стохастический марковский процесс, то такой КА называется вероятностным. В противном случае клеточный автомат является детерминированным. Клеточные автоматы, как правило, объединяют в сети (решетки), так что динамика каждого КА определяется как его текущим состоянием, так и состоянием соседей. Решетки клеточных автоматов (РКА) являются мощным средством моделирования физических процессов в распределенных системах, позволяя получить их временную и пространственную динамику [10]. В них наблюдаются характерные для динамических систем колебательные и волновые явления: периодические, квазипериодические и хаотические колебания [11, 12], распространение волн и волновых фронтов [13], взаимная и вынужденная синхронизация [14], фазовая синхронизация [15].
При использовании РКА для моделирования инфекционных заболеваний узлы решетки имитируют взаимодействие особей популяции в точках пространства, так что на каждом шаге инфекция может с определенной вероятностью перекинуться от больной особи к соседней, тем самым имитируя процессы распространения заболевания [8]. В рамках такого подхода можно представить также процессы воспроизводства и миграции особей [9].
В данной работе предлагается модифицированная 81Я8-модель распространения эпидемий в виде квадратной решетки вероятностных клеточных автоматов, в которой процессы распространения инфекции дополнены процессами динамического регулирования численности популяции. Особенностью предлагаемой модели является использование предложенного Ферхюль-стом [16] механизма ограничения максимальной численности за счет ограниченности общего ресурса, адаптированного с учетом влияния заболевания на процессы воспроизводства.
1. SIRS-модель распространения инфекционных заболеваний
Несмотря на различия, обусловленные спецификой того или иного инфекционного заболевания, большинство из них имеют общие черты: (а) заболевание распространяется через заражение, происходящее при непосредственном (локальном) контакте двух особей: восприимчивой к заболеванию S и инфицированной I, в результате чего восприимчивая особь заболевает и сама становится источником последующих заражений; (б) заболевшая особь со временем излечивается, приобретая при этом иммунитет к последующим заражениям (R); (в) наконец, с течением времени приобретенный иммунитет утрачивается, возвращая особь к исходному состоянию S. Таким образом, мы наблюдаем циклическую цепочку превращений между дискретным и конечным набором состояний S ^ I ^ R ^ S. Отсюда и название данной модели - SIRS. События (а)-(в) не влияют на число особей в популяции, а лишь перераспределяют их количества в каждом из состояний. Помимо них в популяции протекают события, меняющие ее численность -рождение и гибель особей. Гибель особи можно описать как переход из какого-либо состояния X £ {S, I, R} в вакантное состояние V с вероятностью Pd. Рождение новой особи - наоборот, переход из вакантного состояния V в состояние S с вероятностью P,.1 Все события являются
1 Предполагается, что все вновь рожденные являются здоровыми и восприимчивыми, то есть состояния инфицирования и иммунитета не наследуются.
случайными и характеризуются своими значениями вероятностей (Рк), если полагать их мгновенными, или плотностями вероятностей (рк = (Рк/(£), если рассматривать их как непрерывные во времени. Общая схема превращений показана на рис. 1.
При рассмотрении процессов распространения эпидемии интерес представляет не динамика индивидуальных особей, а изменение их количеств (чисел заполнения) (к Е {5, I, Я, V}). Постулируем, что общая численность популяции Nх = Ns + N1 + Ыд ограничена по величине некоторым максимально допустимым значением N, которое назовем емкостью популяции. Величина N обусловлена некоторым общим ресурсом (количеством пищи, площадью территории, е1:с.), сдерживающим неограниченный рост численности. Разность между максимальной и актуальной численностями Nу = N — Nх - есть число «вакансий» в популяции, которые потенциально могут быть заполнены новыми особями.
Помимо превращений, показанных на схеме (см. рис. 1), значительную роль в распространении эпидемий могут играть процессы миграции. Под миграцией понимаем случайные изменения особями своих пространственных координат (подобно диффузии броуновских частиц), в результате чего каждая из особей, независимо от своего состояния {5, I, Я}, перемещается в произвольном направлении с некоторой средней скоростью у. Миграция не влияет напрямую на величины Nk, однако может влиять на них косвенно, через изменение пространственного распределения особей разного вида.
1.1. Схема распространения заболевания без учета изменения численности популяции. Рассмотрим простую модель, отражающую указанные выше свойства, без учета процессов рождения и гибели, которые рассмотрим позднее. Процессы, связанные с инфицированием, заболеванием и излечением, можно интерпретировать как переходы между тремя возможными состояниями {5, I, Я} некоторой гипотетической частицы, проходящей цепочку превращений. Схема превращений, произошедших в течение малого интервала времени ЛЬ, может быть записана в виде системы реакций
5 +1 А 21,
I А я, (1)
Я А 5,
где буквы над стрелками обозначают вероятности соответствующих переходов. Первая реакция описывает инфицирование, вторая - выздоровление, третья - потерю иммунитета. Как видно из схемы (1), все процессы однонаправленны и порядок смены состояний жестко детерминирован, в то время как сами превращения, а следовательно, и моменты времени, в которые они происходят, являются случайными. Динамика данной модели определяется как параметрами Р^, так и численностью популяции, которая является для системы (1) интегралом движения.
1.2. Схема динамического регулирования численности популяции. Для учета динамического характера численности популяции необходимо добавить к системе реакций (1) про-
Рис. 1. Схема превращений: сплошные линии со стрелками показывают направления переходов, штриховые линии показывают влияющие на реакции компоненты; буквы рядом со стрелками обозначают вероятности переходов
Fig. 1. Transformation scheme: solid lines with arrows show the direction of transitions, dashed lines show the components supporting the reaction; letters near arrows indicate the probability of transitions
цессы рождений и гибели особей, которые регулируют ее величину. Как уже указывалось выше (см. рис. 1), гибель представляет собой переход из текущего состояния X £ (5, I, Я} в вакантное: X V, а рождение - переход из вакантного состояния в состояние 5 при наличии репродуктивной особи: V + X р^ 5 + X. Необходимость «вакансии» для рождения новой особи составляет суть механизма ограничения рождаемости, обусловленного ограниченностью ресурсов.
При моделировании динамического регулирования численности популяции необходимо принять во внимание влияние заболевания на процессы воспроизводства. Оно может реализоваться как в виде дополнительной смертности заболевших особей: I ^ V (Р^ - вероятность летального исхода при заболевании), так и в ограничении рождаемости больными (I) и только что переболевшими (Я) особями. В предлагаемой модели реализуется жесткий вариант последнего ограничения, постулирующий, что все заболевшие и выздоровевшие особи не являются репродуктивными. Учет динамического регулирования численности населения приводит к дополнению схемы (1) следующими реакциями:
5 + V РЬ 25, 5, Я ^ V, (2)
I Р£+Р<И V.
1.3. Миграция особей. При распространении эпидемий существенную роль играет миграция - перемещение особей в пространстве. Перемещение может быть представлено как набор реакций, в каждой из которых соседние особи меняются местами
X о Г. (3)
Здесь X и Г - две соседние особи произвольного вида, включая «вакансию». В ходе трансформации (3) частица X встает на место частицы Г, и наоборот. В данном исследовании миграция предполагается случайной и представляет собой диффузионный процесс, характеризуемый коэффициентом диффузии у.
2. Построение модели в форме сети клеточных автоматов
Смоделируем реакции (1), (2), (3) при помощи квадратной решетки вероятностных клеточных автоматов, отображаемой Ь х Ь матрицей М. Положение автомата в сети идентифицируется двумя пространственными индексами: г - номер строки и ] - номер столбца. Соответственно, он имеет четырех соседей, которыми являются узлы с индексами (г — 1, _/), (г, ] + 1), (г + 1, и (г, ] — 1), где сложение и вычитание осуществляется по модулю Ь, то есть используются периодические граничные условия2. Каждый из узлов решетки может находиться в одном из четырех состояний (5, I, Я, V}, которые ассоциированы с описанными выше состояниями особей популяции (5, I, Я}, либо с отсутствием особи V. Полное число клеток естественным образом задает емкость популяции (Ж): N = Ь2.
Моделирование начинается с задания начального состояния решетки М(0), которое затем меняется в соответствии с отображением последования
М(п + 1)= Р (М(п)), (4)
где п = 0,1, 2,... - дискретное время, Р - случайная функция, выбранная таким образом, чтобы имитировать реакции (1), (2), (3). Соответственно, действие Р разбивается на три последовательных этапа: (1) - заболевание, (2) - воспроизводство и (3) - миграция.
2Периодические граничные условия моделируют поведение подгруппы особей, находящихся в окружении популяции неограниченно большого размера.
2.1. Описание работы алгоритма РКА.
(1) Этап заболевания. Выбираются координаты случайной клетки решетки (i, j) и определяется ее текущее состояние M(i, j) G {S, I, R, V}.
(а) При M = S производится случайный выбор одной из четырех соседних клеток:
Mn G {M(i + 1, j), M(i - 1, j),M(i, j + 1),M(i, j - 1)};
если Mn = I, то состояние исходной клетки с вероятностью Pi трансформируется в I: S ^ I.
(б) При M = I состояние исходной клетки с вероятностью P2 трансформируется в R: I ^ R.
(в) При M = R состояние исходной клетки с вероятностью P3 трансформируется в S: R S.
(2) Этап воспроизводства. Выбираются координаты случайной клетки решетки (i, j) и определяется текущее состояние данной клетки M(i, j) G {S, I, R, V}.
(г) При M = V производится случайный выбор одной из четырех соседних клеток:
Mn G {M(i + 1, j), M(i - 1, j),M(i, j + 1),M(i, j - 1)};
если Mn = S, то состояние исходной клетки с вероятностью Рь трансформируется в S:
V ^ S.
(д) При M = S или M = R состояние исходной клетки с вероятностью Pd трансформируется в V: S ^ V или R ^ V.
(е) При M = I состояние исходной клетки с вероятностью Pd + Pdi трансформируется в V: I ^ V.
(3) Этап миграции. Выбирается случайная клетка решетки M и одна из ее четырех соседних клеток Mn, после чего состояния клеток меняются местами: M о Mn.
Трансформации (а)-(е) представляют собой случайные события, характеризуемые соответствующими вероятностями Pk G [0:1]. Для их моделирования используется генератор случайных чисел с равномерным распределением, генерирующий вещественное число в диапазоне [0:1]. Трансформация происходит, если сгенерированное число попадает в интервал [0 : Pk].
На каждом шаге алгоритма Монте-Карло этапы (1) и (2) реализуются N раз, давая возможность каждой клетке решетки принять участие в реакциях, а этап (3) - yN раз, давая возможность каждой частице в среднем совершить у перемещений. Полученная в ходе отображения
(4) последовательность значений M(n) - есть временная эволюция решетки клеточных автоматов. Эта эволюция является случайным процессом, однако, при N ^ те средние концентрации состояний решетки s = Ns/N, i = N//N, r = Nr/N, v = Nv /N будут представлять собой математические ожидания соответствующих состояний, иными словами, будут стремиться к детерминированным величинам.
3. Анализ модели методом среднего поля
В приближении среднего поля для схемы (1,2) можно составить систему уравнений, управляющих изменением числа особей в каждом из состояний М^ (к £ (5, /, Л, V}). Полагая М^ большими, можно составить уравнения для изменений чисел заполнения за малый интервал
( Мт Мт Мя Мт Му Мт V
ДМт " ДТ М + РздТ + "МТ ту - * дт; ( Мт Мт МТ МТ МТ Ч Л +
ллг ( М/ Мд Мд Ч Д
Переходя к пределу ЛЬ —> 0, и учтя, что в + г + г + V = 1, получим систему обыкновенных дифференциальных уравнений для концентраций
в = — ргвг + рзг + рьв (1 — в — г — г) — Рав,
г = Р1вг — Р2г — Раг — Раг, (5)
Г = р2г — р3г — Раг.
Сложив все три переменные, можно также записать уравнение для полной концентрации особей популяции
с = рьв (1 — с) — Рас — да
где с = в + г + г. При отсутствии заболевания (г = 0) мы получаем логистическое уравнение, введенное Ферхюльстом для ограничения численности популяции [16]. В нашем случае это ограничение получается как естественное следствие работы решетки клеточных автоматов и введенного в разделе 1.2. условия, регулирующего процессы воспроизводства.
Проанализируем решения системы (5). Из условий задачи все переменные и параметры являются положительными; кроме того, переменные суммарно ограничены сверху: в + г + г < 1. Поскольку рг входят в уравнения (5) как множители, их одновременное пропорциональное изменение сохраняет качественное поведение системы и меняет лишь временной масштаб. Поэтому, не снижая общности, можно зафиксировать значение одного из параметров и рассматривать изменения остальных по отношению к нему.
Предельными траекториями системы (5) являются три состояния равновесия: Ео = (0, 0, 0), Е = (1 — рь/ра, 0,0) и Е = (в2, ¿2, Г2), где
Р2 + Ра + Ра
в2 =
р1
г2 = в2(рз + Ра)(Ра — Рь + в2Рь) (6)
2 Р2(Р3 — в2Рь) — в2(р1 + Рь)(Рз + Ра),
Г2
в2 Р2(Ра — Рь + в2Рь)
Р2(Р3 — в2Рь) — в2(Р1 + Рь)(Рз + Ра)'
Первое из них (Е0) соответствует случаю полного вымирания популяции. Оно асимптотически устойчиво при Р^ > Рь. Точка Е соответствует полному излечению популяции со временем. Это решение асимптотически устойчиво при
Ра < Рь,
Р2 + Ра + Ра > 1 — Ра —
Р1 Рь.
Наконец, точка Е2 соответствует установлению динамического равновесия между процессами заболевания и излечения, при котором достигается некоторый средний уровень числа заболевших («2), не изменяющийся в дальнейшем. Данное состояние устойчиво, если заменить знак второго неравенства в (7) на противоположный. Как видно из условий устойчивости, при любых значениях параметров Рг > 0 устойчивым является лишь одно из трех состояний равновесия, то есть в фазовом пространстве отсутствует мультистабильность. Соответственно, установившийся в системе ОДУ режим не зависит от выбора начальных условий из области указанных выше допустимых значений.
Таким образом, дифференциальные уравнения (5) позволяют получить общую картину того, как ведет себя моделируемая система в зависимости от управляющих параметров. Насколько эта картина соответствует более реалистичной модели клеточных автоматов, рассматривается далее.
4. Численные исследования модели клеточных автоматов
Проведем численные исследования решетки клеточных автоматов (4) при отсутствии миграции (у = 0) и сопоставим результаты с предсказаниями модели среднего поля (5). Для этого выберем решетку размером 100 х 100 элементов и зададим начальные условия, «разбросав» случайным образом частицы с начальными концентрациями з(0) = 0.5 и ¿(0) = 0.1 так, как это показано на рис. 2, a. Затем, выбрав значения вероятностей реакций, запустим решетку клеточных автоматов и получим временную эволюцию матрицы М(п) (п = 0,1, 2, 3,...).
100
80 60
• ^
40
20 0
a
20
40 60 i
80
100
Рис. 2. Распределение здоровых (черный цвет), больных (темно-серый) и иммунных (светло-серый цвет) особей по решетке: a - начальное состояние при n = 0; b - состояние при n = 1000; значения параметров:Pi = 1, P2 = 0.1, Рз = 0.3, Pb = 0.1, Pd = 0.02, Pdi = 0.02
Fig. 2. Distribution of healthy (black), ill (dark grey) and immune (light grey) individuals on the lattice: a - initial state at n = 0; b - state at n = 1000; parameter values are: Pi = 1, P2 = 0.1, P3 = 0.3, Pb = 0.1, Pd = 0.02, Pdi = 0.02
В результате на решетке устанавливается некоторое динамическое распределение компонентов реакций, пример которого показан на рис. 2, Ь.
4.1. Динамика решетки клеточных автоматов. Посредством усреднения по ансамблю клеток определим для каждого момента времени п текущие значения концентраций з, г и г, получив таким образом зависимость концентраций от времени. Вид траекторий з(п), ¿(п) и г(п) для рассмотренного случая изображен на рис. 3, a (кривые 1 и 1г, соответственно). На том же рисунке штриховыми и штрихпунктирными линиями построены траектории, полученные интегрированием модели среднего поля (5). Как видно из сопоставления графиков, траектории, полученные при имитационном моделировании, качественно схожи с зашумленными траекториями модели среднего поля. Наличие случайных флуктуаций траекторий имитационной модели обусловлено вероятностным характером автомата. Данные флуктуации нельзя «убрать» из модели, но их интенсивность может быть сделана сколь угодно малой за счет усреднения по большому числу клеток ансамбля, так как величина среднеквадратичного отклонения уменьшается с ростом числа клеток, приблизительно как 1/\/Ж (см. рис. 3, Ь).
Как видно из рис. 3, a, между установившимися значениями концентраций имитационной модели и предсказанием среднего поля наблюдается количественное несоответствие, например, величины 3 и оказываются существенно смещенными друг относительно друга. Причиной несоответствия является невыполнение условий равномерности в распределении компонентов, которые образуют на деле кластерную структуру (рис. 2, Ь). Можно показать, что если кластерная структура разрушается, то обе модели дают близкие стационарные значения концентраций. Для этого рассмотрим, как изменится динамика решетки клеточных автоматов при разрушении
200 400 600 800 1000 n
0
a
0.0312 1 | 1 1 | 1 1 1 1 | 1 1 •
0.0156 г
м 0.00781 О :
0.00391 0.00195 - • V» - %W V 1 1 1 1 1 1 1 1 1 1 1 1 1
32 b 64 128 256 512 L
0.125-
'< 0.0625
0.03125
0.125 0.25 0.5 1 2 4 У
Рис. 3. Динамика имитационной модели: a - временные реализации концентраций; b - зависимость среднеквадратичного отклонения as от размера решетки L; c - зависимость разности между s и s2 от у
Fig. 3. Dynamics of the simulation model: a - time series of concentrations; b - dependence of the standard deviation as on the lattice size L; c - dependence of the difference between s and s2 on у
c
кластерной структуры за счет «перемешивания» взаимодействующих частиц. Предположим, что в популяции, помимо описанных в схеме трансформаций (1, 2), присутствует процесс миграции. Миграция не затрагивает количество особей и их распределение по состояниям, а лишь перемешивает компоненты реакций по решетке, разрушая формирующиеся кластеры. Проведенные расчеты подтверждают что, миграция приводит к сближению средних значений концентраций с предсказанием среднего поля. Например, на рис. 3, c построен график зависимости смещения Дз = в — §2 от величины у (в двойном логарифмическом масштабе). Как видно из рисунка, расхождение между моделями монотонно уменьшается с ростом у; при этом величина смещения хорошо аппроксимируется степенной зависимостью Дв ~ у-0'42, которая отображена на графике штриховой линией. Таким образом, при сильном перемешивании, когда компоненты реакций равномерно распределяются по решетке, имитационная модель дает те же значения средних концентраций, что и модель среднего поля. Однако для реальных процессов характерен именно неоднородный, очаговый характер заболеваний, поэтому модель среднего поля оказывается недостаточно точной.
4.2. Поведение модели в зависимости от параметров. Рассмотрим поведение решетки клеточных автоматов в зависимости от управляющих параметров. Как было отмечено выше, в приближении среднего поля рассматриваемая система демонстрирует три установившихся типа поведения: вымирание популяции (Ео), полное излечение популяции (Е), режим хронического заболевания (Е2), границы между которыми определяются соотношениями Ра = Рь и (Р2 + Р^ + Р(0/Р1 + Ра/Рь = 1, соответственно. Все эти режимы наблюдаются и в модели клеточных автоматов. Проведем численные исследования решетки клеточных автоматов и построим данные границы на плоскости параметров, сопоставив их с предсказаниями среднего поля.
Область вымирания популяции определяется условием превышения смертности над рождаемостью: Ра > Рь; соответственно, ее граница задается медианой Рь = Ра при произвольных
значениях остальных параметров. Проведенные численные исследования с моделью клеточных автоматов также дают линейную взаимосвязь между граничными значениями этих параметров (рис. 4, a), однако граница лежит ниже теоретической: Ра ~ аРь, где множитель а принимает значение от 0.6 при отсутствии миграции, до 0.8 при сильной миграции. То есть в рассматриваемой модели для выживания популяции недостаточно равенства смертности и рождаемости, а необходим некоторый «перевес» рождаемости над смертностью.
Область полного излечения задается соотношением: (Р2 + Ра1 + Ра)/Р1 > 1 — Ра/Рь. Рассмотрим границу между зонами заболевания и излечения в зависимости от величин вероятностей выздоровления Р2 и летального исхода болезни Раь для чего проведем численные исследования при Р1 = 1, Рь = 0.1 и Ра = 0.02, и построим график Р2 (РаО, разграничивающий обе зоны (рис. 4, Ь).
На том же графике штриховой линией отметим предсказание модели среднего поля. Как видно из сопоставления графиков, имитационная модель дает линейно убывающую зависимость: Р2 = С — аРаь качественно подобную теоретической, но с множителем а < 1 и дополнительным смещением. В зависимости от скорости миграции, а меняется от 0.5 при у = 0 до 0.9 при у = 4, а величина смещения - от 0.65 при у = 0 до 0.31 при у = 4. Таким образом, экспериментальная граница располагается существенно ниже линии, предсказанной моделью среднего поля. Наличие процесса миграции приводит к сближению обеих границ, но даже при сильной миграции, полностью разрушающей кластеризацию, (у = 4) наблюдается значительное расхождение в предсказаниях обеих моделей.
Таким образом, в модели клеточных автоматов границы между областями с качественно разным поведением смещены относительно линий, предсказанных моделью среднего поля, причем это расхождение не связано с процессами кластеризации, происходящими в пространстве реакций, а имеют, видимо, иную природу. Можно предположить, что причиной данного расхождения является дискретная природа клеточного автомата и вероятностный характер его динамики, проявляющийся в наличии флуктуации концентраций относительно установившихся значений. Вследствие этих факторов, для модели клеточных автоматов должно быть характерно наличие минимального порога концентрации заболевших особей, необходимого для поддержания распространения инфекции в популяции и минимального числа особей для воспроизводства населения. Данное явление известно в популяционной динамике как явление Олли [17,18]. Такой
1 0.8 0.6 0.4 0.2 0
b
-Ч
А.
Ч\
0
0.2
0.4 0.6 р,.
» у=0 = g=1 . g=4
0.8
Рис. 4. Разбиение области параметров на зоны качественно различного поведения: a - разграничение зон вымирания и выживания популяции на плоскости Pb — Pd при разных у; b - разграничение зон полного излечения и заболевания на плоскости Pdi — P2; штриховые линии задают границы, предсказываемые моделью среднего поля
Fig. 4. Division of the parameter area into zones of qualitatively different behavior: a - dying and survival zones on the plane Pb — Pd at different y; b - zones of complete cure and disease on the plane Pdi — P2; dashed lines define the boundaries predicted by the mean-field model
1
a
порог, очевидно, отсутствует в модели среднего поля, где траектория притягивается к устойчивому состоянию равновесия со сколь угодно малого начального значения.
Для того чтобы проверить данную гипотезу рассмотрим, как меняется в ансамбле клеточных автоматов концентрация больных особей г в ходе приближения к границе зоны излечения. Рассчитаем зависимости ¿(Р^) при фиксированном значении Р2 = 0.1 и разных скоростях миграции и сравним их с теоретической зависимостью (6). Вид соответствующих графиков приведен на рис. 5. Из графиков следует, что концентрация заболевших особей в модели клеточных автоматов меньше теоретически предсказанной, но может быть приближена к последней, если в популяции присутствует сильная миграция. Однако, сближение этих зависимостей не ведет к устранению расхождения в расположении бифуркационных линий, поскольку, как видно из графиков, при уменьшении г до критического значения гс (отмечено пунктирной линией) происходит скачкообразное «схлопывание» концентрации заболевших до нуля, после чего процессы, связанные с заболеванием прекращаются. Величина данного порога практически не зависит от скорости миграции и составляет, для рассматриваемого случая ~ 0.015. При этом, в силу того, что скорость изменения г вблизи нуля мала, наличие данного порога приводит к существенной разнице между бифуркационным значениями параметра Раг и его теоретическим предсказанием по формуле (7).
4.3. Режим самоподдерживающихся колебаний. Помимо стационарных режимов, предсказанных моделью среднего поля, в решетке клеточных автоматов наблюдается режим самоподдерживающихся колебаний, отсутствующих в системе дифференциальных уравнений (5). Данный тип поведения возникает на границы зоны излечения, когда значение концентрации заболевших особей приближается к своему пороговому значению гс. Пример таких колебаний приведен на рис. 6, где построены временные реализации г(п), а также - общей населенности
Рис. 5. Зависимость концентрации заболевших особей i от параметра Pdi при разных скоростях миграции у; штриховая линия отмечает предсказание среднего поля
Fig. 5. Dependence of the concentration of diseased individuals i on the parameter Pdi at different migration rates y; dashed line marks prediction of the mean field equations
Рис. 6. Колебания общей численности популяции и относительного числа заболевших особей: a - временные реализации и b - спектр мощности; значения параметров: P1 = 1, P2 = 0.05, P3 = 0.002, Pb = 0.1, Pd = 0.02, Р<ц = 0.02
Fig. 6. Fluctuations of the population and of the relative number of ill individuals: a - time-series and b - power spectrum; parameter values are: Pi = 1, P2 = 0.05, P3 = 0.002, Pb = 0.1, Pd = 0.02, Pdi = 0.02
Рис. 7. Распределение здоровых (черный цвет) и больных (серый цвет) особей по решетке в моменты времени, соответствующие: a - минимуму и b - максимуму заболевания; значения параметров те же, что и на рис. 6
Fig. 7. The distribution of healthy (black) and ill (grey) individuals on the lattice at different time moments corresponding to: a - minimum and b - maximum of the disease; the values of the parameters are the same as those in Fig. 6
x(n) = s(n) + i(n) + r(n) (рис. 6, a); на том же рисунке построен спектр мощности колебаний (рис. 6, b). Как видно из рисунка, колебания численности популяции носят случайный характер, однако, в них присутствует характерный временной масштаб T ~ 300, которому соответствует выраженный пик на нормированной частоте /о ~ 0.0033, хорошо различимый в спектре (см. рис. 6, b). Общая населенность популяции x(n) и концентрация заболевших i(n) колеблются в противофазе: рост заболевания сопровождается снижением численности популяции и наоборот; при этом относительно небольшое изменение числа заболевших особей влечет за собой весьма значительное изменение общей численности.
Колебания численности популяции, происходящие в системе и характеризуемые длительным временным масштабом должны отражать некоторые глобальные процессы, происходящие в популяции в целом. Рассмотрим два «снимка» решетки, соответствующих соседним минимуму и максимуму колебаний, которые приведены на рис. 7, a и 7, b, соответственно. Как видно из сопоставления рисунков, пространственные распределения на максимуме и минимуме населенности являются обратными по отношению друг к другу: области максимальной заселенности на одном из них соответствуют областям минимальной заселенности на другом, и наоборот. Визуально один из рисунков выглядит как негатив второго. Таким образом, в режиме колебаний популяция ведет себя как единый организм, циклически меняющий свое состояние от почти полного выздоровления до последующего рецидива заболевания, происходящего через характерный временной масштаб.
Режим автоколебаний, как уже отмечалось выше, не наблюдается в детерминированной модели (5). При заданных значениях параметров последняя демонстрирует состояние, соответствующее точке равновесия E2. Однако анализ этой точки показывает, что она (а) является устойчивым фокусом и (б) располагается вблизи инвариантного многообразия i = 0. Вследствие этого, в динамической системе в рассматриваемой области под действием шума будет наблюдаться режим возбудимых колебаний. Подобный режим мы наблюдаем и в решетке клеточных автоматов, с тем отличием, что данная модель является автономной и «двигателем» ее колебаний является вероятностный характер ее динамики.
4.4. Колебания численности, вызванные сезонными факторами. В предыдущих разделах мы рассматривали автономную модель распространения инфекции, где все параметры -постоянные числа. Однако, вероятности процессов, связанных с инфицированием и воспроиз-
водством, как правило, не являются стационарными, а меняются вследствие сезонных и иных факторов. Периодические изменения окружающих условий, вызванные сменой времен года, приводят к периодической модуляции параметров Рг вокруг некоторых средних значений [19]. Эта модуляция ведет к сезонным колебаниям заболевания. Рассмотрим эти процессы применительно к рассматриваемой модели.
Предположим, что параметр, отвечающий за вероятность заражения Р1, модулируется периодической (гармонической) функцией времени вокруг среднего значения Р1 по закону: Р1(п) = Р1 (1 + ^ео8(2л/рп)), где /р - частота, а ^ - индекс модуляции. Выберем в качестве среднего значения Р1 = 0.5 и рассмотрим два характерных случая:
1. стационарный режим, изображенный на рис. 3;
2. режим самоподдерживающихся колебаний, изображенный на рис. 6.
В первом случае модуляция параметра приводит к возникновению колебаний, форма которых близка к периодической (см. рис. 8, a). Соответствующий спектр мощности типичен для периодических негармонических колебаний: он содержит пик на частоте модуляции /р, а также - на гармониках 2/р, 3/р, ... (рис. 8, Ь). Отметим, что небольшие по амплитуде колебания относительного числа инфицированных особей ведут к значительным колебаниям численности популяции в целом, которая на рис. 8, a отмечена штриховой линией. При этом период модулирующего сигнала существенным образом определяет амплитуду возникших колебаний. При очень быстрой модуляции (Т ~ 1) колебания практически незаметны. С ростом периода амплитуда колебаний растет, вплоть до некоторой резонансной частоты /о ~ 0.0029, после чего наблюдается ее уменьшение. Таким образом, амплитудно-частотная характеристика системы обнаруживает резонансный характер (рис. 8, с).
Если выбрать значения параметров из области, где наблюдаются собственные колебания, то возникают биения между двумя колебательными процессами: собственными и параметриче-
1
0.9 0.8 0.7 0.6 ¿0.5 0.4 0.3 0.2 0.1 0
0
2000
4000 6000 n
8000 10000
0.002
0.004 _ 0.006 f
0.008 0.01
0.005
0.01
f
p
0.015
0.02
Рис. 8. Параметрические колебания при периодической модуляции параметра P1 (n) = = 0.5 (1 + 0.4 cos (2я/рп)), где fp = 1/600: a - временные реализации; b - спектр мощности; c - зависимость амплитуды колебаний от частоты воздействия
Fig. 8. Parametric oscillations with periodic modulation of parameter P1 (n) = 0.5 (1 + 0.4 cos (2nfpn)), where fp = 1/600: a - time-series; b - power spectrum; c - dependence of oscillation amplitude on influence frequency
b
a
0
c
Рис. 9. Параметрические колебания при периодической модуляции параметра Pi(n) = 0.5 (1 + 0.2 cos (2nfpn)), где fp = 1/300: a - временные реализации и b - спектр мощности
Fig. 9. Parametric oscillations with periodic modulation of parameter P1(n) = 0.5 (1 + 0.2 cos (2nfpn)), where fp = 1/300: a - time-series and b - power spectrum
скими колебаниями. В результате мы наблюдаем почти квазипериодические колебания с двумя характерными масштабами, как это показано на рис. 9, a. В спектре таких колебаний присутствует пик как на частоте модуляции /р, так и на частоте собственных колебаний (рис. 9, b). Вследствие этого, помимо «быстрых» осцилляций с временным масштабом, близким к периоду сезонных изменений, наблюдаются «медленные» колебания огибающей с разностной частотой А/ = I/o — fp|. На практике эти колебания проявляют себя как почти периодические вспышки инфекции, происходящие без видимых причин, с очень длительным интервалом между ними. Временной масштаб, связанный с изменением огибающей, может намного превышать годичный цикл вспышек заболевания, в результате чего раз в несколько лет (или десятков лет) может наблюдаться значительный рост заболевания.
Заключение
Предложенная SIRS модель клеточных автоматов с динамическим регулированием численности популяции демонстрирует как стационарные режимы, соответствующие устойчивым состояниям равновесия в модели среднего поля, так и отсутствующий в ней режим самоподдерживающихся колебаний с выраженной периодической составляющей. Данный режим аналогичен колебаниям, наблюдаемых в возбудимых системах ОДУ, находящихся под действием внешнего шума.
Рассматривая система демонстрирует также бистабильность, отсутствующую в моделях ОДУ: при слишком малой концентрации заболевших особей система с высокой вероятностью демонстрирует быстрое прекращение эпидемии; в то же время, если эта концентрация превышает некоторый «порог», распространение заболевания продолжается неограниченно долго. Данный эффект является следствием дискретного характера системы клеточных автоматов. Он также хорошо известен из экспериментальных наблюдений за распространением эпидемий в природных популяциях (см. [3]).
При рассмотрении неавтономной SIRS-модели, когда вероятность инфицирования периодически модулируется вследствие сезонных факторов, наблюдаются периодические или квазипериодические колебания численности популяции значительной амплитуды. Зависимость амплитуды колебаний от периода модуляции носит выраженный резонансный характер.
Библиографический список
1. Бейли Н. Математика в биологии и медицине. Москва: Мир, 1970.
2. Hethcote H.W. The mathematics of infectious diseases // SIAM Review. 2000. Vol. 42, № 4. P. 599-653.
3. Андерсон Р., Мэй Р. Инфекционные болезни человека. Динамика и контроль. Москва: Мир, 2004.
4. Kermack W., McKendrick A. A contribution to the mathematical theory of epidemics // Proc. R. Soc. 1927. Vol. A115. P. 700-721.
5. Кобринский Н.Е., Трахтенберг Б.А. Введение в теорию конечных автоматов. Москва: Физ-матгиз, 1962.
6. Цетлин М.Л. Исследования по теории автоматов и моделированию биологических систем. Москва: Наука, 1969.
7. Тоффоли Т., Марголус Н.Машины клеточных автоматов. Москва: Мир, 1991.
8. Boccara N., Cheong K. Automata network SIR models for the spread of infectious diseases in populations of moving individuals // Journal of Physics A: Mathematical and General. 1992. Vol. 25, № 9. P. 2447.
9. Sirakoulis G.C., Karafyllidis I., Thanailakis A. A cellular automaton model for the effects of population movement and vaccination on epidemic propagation // Ecological Modelling. 2000. Vol. 133, № 3. P. 209-223.
10. Ванаг В.К. Исследование пространственно распределенных динамических систем методами вероятностного клеточного автомата // УФН. 1999. Т. 169, № 5. С. 481-505.
11. Provata A., Nicolis G., Baras F. Oscillatory dynamics in low-dimensional supports: A lattice Lotka-Volterra model // J of Chem. Phys. 1999. Vol. 110. P. 8361-8368.
12. Shabunin A., Baras F., Provata A. Oscillatory reactive dynamics on surfaces: a lattice limit cycle model // Physical Review E. 2002. V. 66, no 3. P. 036219.
13. Tsekouras G., Provata A., Baras F. Waves and their interactions in the lattice Lotka-Volterra model // Известия высших учебных заведений. Прикладная нелинейная динамика. 2003. Т. 11, №2. С. 63-71.
14. Wood K., Van den Broeck C., Kawai R., Lindenberg K. Universality of Synchrony: Critical behavior in a discrete of stochastic phase-coupled oscillators model // Physical Review Letters. 2006. Vol. 96, no. 14. P. 145701.
15. Efimov A., Shabunin A., Provata A. Synchronization of stochastic oscillations due to long-range diffusion// Physical Review E. 2008. Vol. 78, no. 5. P. 056201.
16. Verhulst P. Notice sur la loi que la population suit dans son accroissement // Corr. Math. et Phys. 1838. Vol. 10. P. 113-121.
17. Allee W. C. Cooperation Among Animals with Human Implications. New York: Henry Schuman, Publisher, 1951.
18. Неверова Г.П., Хлебопрос Р.Г., Фрисман Е.Я. Влияние эффекта Олли на динамику популяций с сезонным характером размножения // Биофизика. 2017. Т. 62, № 6. С. 1174-1184.
19. Glendinning P., Perry L.P. Melnikov analysis of chaos in a simple epidemiological model // J. Math. Biol. 1997. Vol. 35. P. 359-373.
References
1. Bailey N. The Mathematical Approach to Biology and Medicine. London: John Wiley and Sons, 1967.
2. Hethcote H. W. The mathematics of infectious diseases. SIAM Review, 2000, vol. 42, no. 4, pp. 599-653.
3. Anderson R., May R. Infectios Deseases of Humans - Dynamics and Control. Oxford: Oxford University Press, 1991.
4. Kermack W., McKendrick A. A contribution to the mathematical theory of epidemics. Proc. R. Soc., 1927, vol. A115, pp. 700-721.
5. Kobrinsky N.E., Trahtenberg B.A. Introduction to the Theory of Finite Automata. Moscow: Fizmatgiz, 1962 (in Russian).
6. Tzetlin M.L. Research on the Theory of Automata and Modeling of Biological Systems. Moscow: Nauka, 1969 (in Russian).
7. Toffoli T., Margolus N. Cellular Automata Machines: A New Environment for Modeling. Cambridge: MIT Press Cambridge, 1987.
8. Boccara N., Cheong K. Automata network SIR models for the spread of infectious diseases in populations of moving individuals. Journal of Physics A: Mathematical and General, 1992, vol. 25, no. 9, p. 2447.
9. Sirakoulis G.C., Karafyllidis I., Thanailakis A. A cellular automaton model for the effects of population movement and vaccination on epidemic propagation. Ecological Modelling, 2000, vol. 133, no. 3, pp. 209-223.
10. Vanag V.K. Study of spatially extended dynamical systems using probabilistic cellular automata. Physics-Uspekhi, 1999, vol. 42, no. 5, p. 413.
11. Provata A., Nicolis G., Baras F. Oscillatory dynamics in low-dimensional supports: A lattice Lotka-Volterra model. JofChem. Phys. 1999, vol. 110, pp. 8361-8368.
12. Shabunin A., Baras F., Provata A. Oscillatory reactive dynamics on surfaces: a lattice limit cycle model. Physical Review E, 2002, vol. 66, no. 3. 036219.
13. Tsekouras G., Provata A., Baras F. Waves and their interactions in the lattice Lotka-Volterra model. Izvestiya VUZ. Applied Nonlinear Dynamics, 2003, vol. 11, no. 2, pp. 63-71.
14. Wood K., Van den Broeck C., Kawai R., Lindenberg K. Universality of synchrony: Critical behavior in a discrete of stochastic phase-coupled oscillators model. Physical Review Letters, 2006, vol. 96, no. 14, 145701.
15. Efimov A., Shabunin A., Provata A. Synchronization of stochastic oscillations due to long-range diffusion. Physical Review E, 2008, vol. 78, no. 5, 056201.
16. Verhulst P. Notice sur la loi que la population suit dans son accroissement. Corr. Math. et Phys. 1838, vol. 10, pp. 113-121.
17. Allee W. C. Cooperation Among Animals with Human Implications. New York: Henry Schuman, Publisher, 1951.
18. Neverova G.P., Khlebopros R.G., Frisman E.Y. The allee effect in population dynamics with a seasonal reproduction pattern. Biophysics, 2017, vol. 62, no. 6, pp. 967-976.
19. Glendinning P., Perry L. P. Melnikov analysis of chaos in a simple epidemiological model. J. Math. Biol. 1997, vol. 35, pp. 359-373.
Шабунин Алексей Владимирович - окончил Саратовский государственный университет (1990). Доктор физико-математических наук (2008), профессор кафедры радиофизики и нелинейной динамики СГУ. Научные интересы - нелинейная динамика, теория колебаний, синхронизация и управление хаосом. Автор более 75 научных публикаций.
410012 Саратов, ул. Астраханская, 83 Саратовский государственный университет E-mail: shabuninav@info.sgu.ru
