Влияние пассивных элементов на синхронизацию осцилляторных ансамблей

Прикладные задачи

^^^^^^^^^^»нелинейной теории колебаний и вслн

Изв. вузов «ПНД», т. 18, № 3, 2010 УДК 621.391.01

ВЛИЯНИЕ ПАССИВНЫХ ЭЛЕМЕНТОВ НА СИНХРОНИЗАЦИЮ ОСЦИЛЛЯТОРНЫХ АНСАМБЛЕЙ

В.С. Петров, Г.В. Осипов

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

Ключевые слова: Синхронизация, автоколебательный ансамбль, пассивный элемент.

Введение

Изучение коллективной динамики больших систем связанных элементов является на сегодняшний день одной из наиболее важных задач, которые возникают в физике, химии, биологии. В таких системах могут наблюдаться всевозможные динамические эффекты, например, синхронизация, распространение плоских, концентрических и спиральных волн, а также формирование различных пространственно-временных структур. Несмотря на то, что большинство естественных систем являются неоднородными, то есть состоят из различных элементов, отличающихся индивидуальной динамикой, на данный момент большая часть опубликованных работ посвящена исследованию динамики однородных сред [1-5]. Тем не менее изучение именно неоднородных систем, представляющих собой смеси пассивных, возбудимых и автоколебательных элементов, является важным и интересным. Примером такой системы является сердце. В сердце клетки рабочего миокарда, так называемые кардиомиоциты, являются возбудимыми; клетки, принадлежащие водителям ритма, являются автоколебательными; кроме этого, существует особый тип клеток -фибробласты, которые являются пассивными. Автоколебательные элементы могут генерировать периодический ритм независимо от внешних условий. Разница между

возбудимыми и пассивными элементами состоит в том, что первые могут производить так называемый потенциал действия в ответ на входящий стимул, в то время как последние лишь релаксируют к состоянию равновесия. С точки зрения качественной теории динамических систем разница сводится к наличию у возбудимых элементов порогового множества (например, кривая медленных движений или входящая сепаратриса седла). Очевидно, что только рассмотрение совокупности всех указанных элементов, а также особенностей их взаимодействия позволяет получить наиболее адекватные результаты при исследовании динамики подобных биологических систем. Заметим, что такие динамические процессы, как распространение концентрических или спиральных волн, имеют прямое отношение к работе сердца [6]. В частности, в отсутствие каких-либо патологий сино-атриальный узел -один из водителей ритма в сердце - генерирует концентрические волны электрической активности. Когда фронт такой волны проходит через клетки рабочего миокарда (кардиомиоциты), он инициирует сокращение мышечной ткани. Однако при развитии патологий может произойти разрушение фронта концентрической волны и рождение спиральной волны или хаоса спиральных волн. В сердце такие эффекты соответствуют различным аритмиям: тахикардии и фибрилляции, борьба с которыми является крайне актуальной задачей. Совмещение подходов биологии и нелинейной динамики позволяет наиболее эффективно исследовать подобные эффекты.

Одним из весьма важных вопросов динамики разнородных ансамблей является проблема синхронизации [7]. Так, в частности, пассивные элементы могут оказывать существенное влияние на свойства синхронизации в осцилляторных ансамблях, представляющих собой совокупность связанных неидентичных автоколебательных элементов [8]. Эта проблема находится в центре рассмотрения данной статьи. Кроме того, приложение подобной задачи может быть найдено при рассмотрении функционирования синусного узла сердца. Синусный узел состоит из автоколебательных и пассивных сердечных клеток - пейсмейкеров и фибробластов, соответственно, причем число последних может достигать 70%. Необходимым условием нормальной работы сердца является синхронизация клеток в пределах синусного узла, так как именно эти клетки формируют ритм всего сердца в целом. Поэтому данная задача является важной и актуальной.

В данной работе представлены результаты численных и аналитических исследований с использованием модели Бонхоффера-ван дер Поля [9]. В контексте моделирования сердечной активности данная модель является простой, но тем не менее достаточно точной аппроксимацией сердечной клетки. Кроме этого, важным достоинством данной модели является возможность моделирования как пассивных, так и автоколебательных элементов одним и тем же набором уравнений, что заметно упрощает исследования.

Данная статья состоит из двух частей. В первой части представлены результаты исследований трех связанных элементов Бонхоффера-ван дер Поля: двух автоколебательных и одного пассивного. Здесь рассматривается влияние пассивного элемента на порог синхронизации автоколебательных элементов. Во второй части статьи полученные результаты перенесены на случай распределенной двумерной системы.

1. Синхронизация двух автоколебательных элементов под влиянием пассивного

Прежде чем перейти к рассмотрению трех связанных элементов, остановимся коротко на одном элементе Бонхоффера-ван дер Поля и различных динамических режимах, наблюдаемых в данной модели.

1.1. Модель Бонхоффера-ван дер Поля. Используемая нами при исследовании модель Бонхоффера-ван дер Поля представляет собой систему обыкновенных дифференциальных уравнений второго порядка:

х = х — х3/3 — у,

у = е(х — у + а).

При моделировании динамики сердечных клеток динамическим переменным х, у можно придать физический смысл потенциала действия мембраны клетки и суммарного ионного тока, протекающего через нее, соответственно. Для введения понятий автоколебательного, возбудимого и пассивного элементов рассмотрим фазовое пространство системы (1). Изоклина вертикальных наклонов у = х — х3/3 (кривая медленных движений) изображена сплошной линией на рис. 1. Изоклины горизонтальных наклонов у = х + а для различных значений параметра а изображены штрихпунктирными линиями. Параметр а является контрольным и позволяет получать различные динамические режимы в данной модели. Если |а| < (1 — е)3/2/3, то система является автоколебательной. В этом случае имеется единственное неустойчивое состояние равновесия, расположенное на неустойчивой ветке кривой медленных движений (отрезок АС) и устойчивый предельный цикл (последовательность АБСБ на рисунке). Если |а| > 8/3, то элемент является пассивным. Для всех остальных значений а элемент является возбудимым. Разница между пассивным и возбудимым режимом состоит в данном случае в расположении состояния равновесия.

Для возбудимого режима оно находится на отрезке АБ (СБ) (см. рис. 1). В этом случае достаточно большой внешний стимул, прикладываемый к первому уравнению системы (1), выводит изображающую точку вправо (влево) от кривой медленных движений, являющейся в данном случае пороговым множеством. В результате этого, прежде чем вернуться в состояние равновесия система пройдет по участку БС (АБ) кривой медленных движений и произведет, таким образом, потенциал действия. В случае же пассивного элемента потенциал действия не производится, так как система возвращается в состоя-

-3.0 -2.0 -1.0 0 1.0 2.0 х

Рис. 1. Фазовое пространство одного элемента Бонхоффера-ван дер Поля. Сплошная линия - кривая медленных движений у = х — х3/3, наклонные штрихпунктирные - изоклины горизонтальных наклонов у = х + а для различных значений а: а1 = 1/3(1 — е)3/2, а2 = —1/3(1 — е)3/2, аз = —8/3, а4 = 8/3

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

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

X1 = Х1 - х\/3 - У1 + ёра(Хр - Х\) + йоо(Х2 - Х\),

У1 = е(х\ - ух + а1),

X2 = Х2 - х2/3 - У2 + dpo(Хр - Х2) + (ооХ - Х2),

У2 = е(Х2 - У2 + а2),

Хр — Хр Хр/3 ур + dop(x1 + Х2 2хр) 1

ур = £(Хр - Ур + ар)-

Здесь (х1,у1) и (х2,у2) - переменные автоколебательных элементов, для которых \аг\ < 1/3(1 - е)3/2; (хр, ур) - переменные пассивного элемента, причем ар < -8/3. Три коэффициента d оо; (ор> (ро описывают силу диффузионной связи: между автоколебательными элементами, со стороны автоколебательных элементов на пассивные и со стороны пассивных элементов на автоколебательные, соответственно. Для простоты начальных исследований ограничимся следующим предельным случаем: (1оо = dop = 0, dpo = d. Это соответствует несвязанным автоколебательным элементам и предельному случаю односторонней связи со стороны пассивного элемента (ор = 0, что упростит дальнейшие аналитические выкладки. Пассивный элемент при этом находится в устойчивом состоянии равновесия хр = Р = - ^3ар. В численных экспериментах Р = 2.5. Параметр (ро = d является контрольным и меняется в ходе численных экспериментов. Более общий случай ненулевых значений (оо и (ор рассмотрен далее.

Для исследования поставленной задачи был выполнен следующий численный эксперимент: при фиксированном значении а1 = 0.3 и различных значениях параметров а2 € [0.3, 0.35], d € [0, 0.2] вычислялись средние частоты автоколебательных элементов Ш1, Ю2. Затем для каждой пары (а2, d) вычислялась эффективная частотная расстройка между этими элементами. Результаты вычислений представлены на рис. 2. Здесь по оси абсцисс отложены значения параметра связи с пассивным элементом d от 0 до 0.2, по оси ординат отложена разность а1 - а2, определяющая различные значения частотной расстройки между автоколебательными элементами в отсутствие связи с пассивным элементом. Рассмотренный диапазон изменения параметра а2 охватывает все возможные частотные расстройки для элемента Бонхоффера-ван дер Поля. Градации серого цвета представляют различные значения эффективной частотной расстройки \ю1 - Ю2\. Из рисунка видно, что введение связи со стороны пассивного элемента может приводить как к уменьшению, так и

¡а^-аъ^Ю"3

ага2 -'-'-:- """ 3.0

0 0.05 0.10 0.15 <3

Рис. 2. Зависимость эффективной частотной расстройки |ю1 — ю2| между двумя автоколебательными элементами (градации серого цвета) от разности параметров а1 — а2 и от связи с пассивным элементом й. Наблюдается увеличение и уменьшение |ю1 — ю2| в зависимости от й. Прямая Ь1 соответствует переходу второго автоколебательного элемента в возбудимый режим через субкритическую бифуркацию Андронова-Хопфа. Прямая Ь2 - аналитическая аппроксимация той же зависимости в кусочно-линейном случае. Белые контурные линии - аналитическая аппроксимация зависимости |ю1 — ю2| от а1 — а2 и й в кусочно-линейном случае. Прямая Ь3 - соответствует соотношению ю1 (а1 — а2,й) = ю2(а1 — а2,й) в кусочно-линейном случае.

к увеличению эффективной частотной расстройки между автоколебательными элементами. Так, например, для а\ — а2 = 0.325 наблюдается уменьшение |ш1 — Ю2| от примерно 0.002 при 1 = 0 практически до нуля при 1 ~ 0.06 и дальнейшее увеличение расстройки до 0.003 при 1 & 0.15. Также очевидно, что при а1 — а2 = 0 значение |ш1 — Ю2| равняется нулю для любого значения 1 (горизонтальная область черного цвета), однако следует заметить, что сами значения частот Ш1, Ю2 меняются при изменении 1. Прямая Ь1 разделяет плоскость параметров (а1 — а2,1) на две области: ниже ¿1, где оба элемента являются автоколебательными, и выше ¿1, где один или оба элемента перешли в возбудимый режим. Так как а2 < а1 = 0.3, то при заданных а1 и а2 второй элемент всегда переходит в возбудимый режим при меньших значениях 1 по сравнению с первым, поэтому ¿1 соответствует переходу второго элемента из автоколебательного режима в возбудимый. Линия ¿1 была получена аналитически. Для ее получения и для исследования бифуркации перехода элемента от автоколебательного режима к возбудимому из-за влияния пассивного элемента рассмотрим систему вида

Хо Х0 Х0 /3 — Уо + 1(хр — Хо),

Уо = е(хо — Уо + ао),

• = 3

X р — Хр Хр/3 ур)

ур = £(хр — Ур + ар). Так как первый и второй элементы в системе (2) не связаны между собой, то переход к системе (3), состоящей только из двух элементов - автоколебательного (хо, уо) и

пассивного (хр,ур), никак не влияет на исследуемую бифуркацию и искомую линию Ь\. Полагая правые части системы (3) равными нулю, получим уравнение для координаты х0 состояния равновесия

х^ + 3йх* + 3(а0 — йР) = 0. Уравнение (4) имеет единственный действительный корень

х* = 1 № — ^

Хо = 2 * ,

где

I = — 12(а* — йР) + 4й3 + 9(а* — йР )2.

Устойчивость состояния равновесия (5) определяется корнями Х12 характеристического уравнения

X2 + Х

Действительная и мнимая части корней Х12 уравнения (7) для фиксированных значений а* = 0.3, Р = 2.5 представлены на рис. 3. Из рисунка видно, что изначально неустойчивый фокус (Б,е(Х1)2) > 0, 1ш(Х1)2) = 0, й £ [0, 0.041]) с ростом й становится сначала неустойчивым узлом (Б,е(Х1)2) > 0, 1ш(А,1>2) = 0, й £ [0.041, 0.22]), а затем снова неустойчивым фокусом. Затем для й ~ 0.29 реальные части двух комплексно-сопряженных корней уравнения (7) обращаются в нуль и становятся отрицательными при дальнейшем росте й. Таким образом, состояние равновесия системы (3) приобретает устойчивость в результате бифуркации Андронова-Хопфа. Детальный анализ показал, что данная бифуркация является субкритической и сопровождается рождением неустойчивого предельного цикла. В итоге, при й ~ 0.298 устойчивый и неустойчивый циклы сливаются и исчезают в результате бифуркации слияния предельных циклов. После этого элемент переходит в возбудимый режим, и колебания в системе (3) вымирают [10].

Для построения ¿1 (см. рис. 2) достаточно построить кривую, определяемую уравнением

Не

что и было сделано для второго автоколебательного элемента системы (2).

Рис. 3. Реальная (а) и мнимая (б) части корней Х1>2 характеристического уравнения (7) для различных значений параметра й при фиксированных а0 = 0.3, Р = 2.5. В результате субкритической бифуркации Андронова-Хопфа состояние равновесия системы (3) приобретает устойчивость при й ~ 0.29

Таким образом, было получено, что введение пассивного элемента может приводить как к увеличению, так и к уменьшению эффективной частотной расстройки между автоколебательными элементами (см. рис. 2). Приведем также аналитическое описание данного эффекта.

1.3. Три связанных элемента. Кусочно-линейная аппроксимация. Рассмотрим систему трех связанных элементов Бонхоффера-ван дер Поля в кусочно-линейной аппроксимации

где

Р (Хг,Уг) = <

х1 = Р (х1,у{) + й(хр — х1),

У1 = е(х1 — у1 + 01),

Х2 = Р(х2,У2) + й(хр — х2),

у2 = е(х2 — У2 + 02),

хр = Р (хр,Ур),

ур = е(хр — Ур + ар),

4 зхг — 2 — Уг, хг < — 1,

2 зхг — Уг, \хг\ < 1,

V 4 зхг + 2 — Уг, хг > +1,

1, 2, р.

Ограничимся случаем односторонней связи со стороны пассивного элемента, тогда последний все время будет находиться в состоянии равновесия хр = Р. Кривая медленных движений Р(х) + й(Р — х) автоколебательных элементов изображена сплошной линией на рис. 4. Заметим, что данная кривая изменяется с ростом й,

в частности, при й = 2/3 она не имеет локальных максимумов и в этом случае не может быть использована в качестве аппроксимации элемента Бонхоффера-ван дер Поля. Однако мы ограничиваемся диапазоном изменения параметра й от 0 до 0.2 и, таким образом, не выходим за рамки кусочно-линейного приближения.

Найдем частоту колебаний первого элемента. Для этого в системе (9) сделаем замену времени т = Ь/е. Тогда, используя соотношение Р(х\,у\) + +й(Р — Х1) = 0 и решая систему (9), найдем траекторию движения первого элемента У1(Ь) по правой ветви кривой

чз

\ (-8/(4+3(1), У \ (0 ¿Р) '-1)+2/3) (1ДР-1)+2/3) Уи/\

(0,0) \

(-1ДР+1)-2/3) (8/(4+ЗЙ Уа Л '), 4^+1)-2/3 )\

х

Рис. 4. Кривая медленных движений Е(х) + +й(Р — х) кусочно-линейного автоколебательного элемента Бонхоффера-ван дер Поля, находящегося под воздействием пассивного элемента с состоянием равновесия (хр = Р,ур = Р + ар) и коэффициентом связи й

медленных движений

У1

где

. = 3й

* = 7 + 3й

и

= 4 + 3й

П = . (13) 7 + 3й

Используя решение (11) и условия

У1 (0)= Ул = (Р + 1)й — 2/3,

У1 (Тг)= Уи = (Р — 1)й + 2/3, найдем время движения по правой ветви кривой медленных движений:

? — У л

Тг = п 1п

I — У и

Аналогично находится время движения первого элемента по левой ветви кривой медленных движений

гт , V — Уи

Т1 = п 1п-,

V — У Л

где

3й(Р + а1) + 4а1 — 6

V = —----.

7 + 3й

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

е

Ы1

Т + Т

Аналогичной формулой дается частота колебаний второго элемента. Тогда, используя (15), можно получить зависимость ю1(а1 — а2,й) — ю2(а1 — а2,й). Линии уровня данной поверхности построены на рис. 2 белыми линиями. Как видно из рисунка, аналитическое приближение с использованием кусочно-линейной модели Бонхоффера-ван дер Поля дает достаточно точное описание результатов, полученных в численном эксперименте.

В данном приближении также имеет место вымирание колебаний при достаточно больших значениях коэффициента связи с пассивным элементом. Условием вымирания автоколебательного элемента в данном случае является касание изоклины горизонтальных наклонов У = 0 с вершиной (1, й(Р—1)+2/3) кривой медленных движений (см. рис. 4):

а =

Линия ¿2 на рис. 2, построенная согласно (16), является аналогом ¿1 на том же рисунке. Разница между ними объясняется тем, что К-образная изоклина горизонтальных наклонов в нелинейном и кусочно-линейных случаях различным образом изменяется при увеличении влияния пассивного элемента.

На плоскости параметров (а\ — а2,() получим также прямую, на которой эффективная частотная расстройка между автоколебательными элементами равна нулю, то есть Ю1(а1 —а2, ()—ю2(а1—а2, () = 0. Используя полученное ранее выражение для частоты (15), получим

— Уй VI — Уи = ^2 — Уй У2 — Уи

— Уи VI — Уй Ъ — Уи V2 — Уй

Решая это квадратное уравнение относительно а2, получим два корня

а2 = а1, а2 = 2(Р — а1.

Первый корень - прямая а1 — а2 = 0, ( € [0, 0.2], представлен горизонтальной областью черного цвета на рис. 2. Второй корень, переписанный в виде а1 — а2 = = 2(а1 — (Р), при а1 = 0.3 представляет собой линию Ь3 на рис. 2. Как видно, данная линия очень хорошо ложится на результаты численных экспериментов (узкая вытянутая черная область на рис. 2, расположенная под углом).

1.4. Синхронизация в ансамбле из трех элементов. Ранее нами рассматривался случай (оо = 0, соответствующий не связанным осцилляторным элементам. В данной ситуации говорить о синхронизации в системе не имеет смысла, поскольку связь между элементами отсутствует. Поэтому до сих пор речь шла лишь об изменении эффективной частотной расстройки между двумя автоколебательными элементами при введении влияния со стороны пассивного элемента. В данном разделе мы рассмотрим систему (2) в случае ( оо = 0. В частности, нас интересует зависимость порогового значения коэффициента связи между автоколебательными элементами (оо, при котором возникает синхронизация, от связи с пассивным элементом (ро. Данное пороговое значение связи будем обозначать .

Ранее было показано, что введение влияния пассивного элемента в системе (2) приводит сначала к уменьшению эффективной частотной расстройки между автоколебательными элементами до нуля и при дальнейшем росте ( ро - к ее возрастанию до исходной величины и даже более. Этот эффект имеет место при любых значений параметров а1, а2, поэтому в данном разделе мы без потери общности ограничимся конкретными значениями а\ =0.3, а2 = 0.25.

На рис. 5 представлены зависимости порога синхронизации от связи с пассивным элементом для трех значений обратной связи со стороны автоколебательных элементов (ор: 0; 0.05; 0.1. Как видно из графика для ( ор = 0, увеличение влияния пассивного элемента приводит к снижению порога синхронизации до нуля ((ро = 0.11). Синхронизация автоколебательных элементов при бесконечно малом значении связи между ними = 0) в данном случае возможна за счет того, что эффективная частотная расстройка между элементами

0.02

-ЛЛ- d0р=0 -о- d0р=0.1 ^ ¿ор=0.0

5

0

0.1

0.2

Рис. 5. Зависимости порога синхронизации d r от связи с пассивным элементом для трех значений обратной связи со стороны автоколебательных элементов dop

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

Ранее, во всех предыдущих разделах нами рассматривался только предельный случай одностороннего влияния со стороны пассивного элемента на автоколебательные. Определим теперь, как влияет обратная связь со стороны автоколебательных элементов (dop) на характер зависимости dthr(dpo). Очевидно, что для исследуемой системы (2) даже при doo = 0 введение ненулевого значения dop означает введение ненулевой связи между автоколебательными элементами посредством пассивного. В этом случае возможно возникновение синхронизации в системе для doo = 0 в некотором диапазоне значений dpo. Это подтверждается соответствующей кривой на рис. 5 (круглый маркер). Здесь синхронизация с нулевым порогом dthr = 0 имеет место для dpo Е [0.09; 0.15]. Заметим также, что увеличение связи dop приводит к увеличению полосы синхронизации с нулевым порогом dthr.

Таким образом, введение обратной связи со стороны автоколебательных элементов на пассивный приводит к появлению области параметра dpo, в которой синхронизация в системе возможна даже при нулевых значениях связи doo. Заметим также, что, несмотря на это, качественная картина зависимостей на рис. 5 сохраняется, а именно, наблюдается уменьшение порога синхронизации dthr до нуля с ростом влияния пассивного элемента dpo, и затем возрастание dthr с увеличением dpo.

2. Синхронизация в двумерном автоколебательном ансамбле

В данном разделе результаты, полученные для трех связанных элементов, перенесены на случай двумерных решеток. Отметим, что последние являются аппроксимацией некоторых реальных автоколебательных сред, например, синусного узла сердца, вытянутого в двух измерениях и имеющего маленькую толщину. Топология рассматриваемой системы представлена на рис. 6. Данная топология представляет собой две регулярные N х N N = 120) решетки, расположенные параллельно друг другу. Белым цветом изображена решетка пассивных элементов с хр = Р = 2.5, серым цветом - решетка автоколебательных элементов со случайным разбросом сц Е [0.05; 0.25], г = 1,..., N * N. Линии соответствуют связям между элементами. Нотация коэффициентов (1оо, (рр, (ор, (ро осталась прежней.

Для исследования синхронизации в данной системе были выполнены следующие расчеты. Для случая односторонней связи со стороны пассивных элементов ((ор = (рр = 0) и связи между автоколебательными элементами (оо = = 0.002 выполнялось численное интегрирование системы на временном интервале £ е [0; 2 000 000] для различных значений параметра (ро. Для каж-

Рис. 6. Топология рассматриваемой системы: две

дого случая выполнялся расчет частот регулярные N х N (М = 120) решетки, располо-колебаний элементов второй решетки. женные параллельно друг другу

Результаты данных вычислений представлены на рис. 7. Здесь для каждого рассмотренного значения (ро изображены значения частот всех автоколебательных элементов, кроме этого для трех значений (ро представлены распределения N (ю) числа автоколебательных элементов, имеющих частоту ю. Из рисунка видно, что изначально ((ро = 0) синхронизации в системе нет, о чем говорит достаточно широкое распределение N(ю) (вставка (а)). Заметим, однако, что это распределение не является равномерным, несмотря на случайность распределения параметра а между автоколебательными элементами. Это обуславливается наличием ненулевого параметра связи (оо = 0.002. Увеличение воздействия со стороны пассивных элементов приводит к сближению частот колебаний автоколебательных элементов, то есть к возникновению все более синхронного поведения в системе. Так, например, для (ро = 0.03 распределение N(ю) становится заметно уже и выше (вставка (б) на рис. 7), что говорит о проявлении синхронизации. В итоге для достаточно больших ( ро в системе устанавливается режим глобальной синхронизации, при котором все элементы осциллируют с одной частотой. Об этом свидетельствует, например, точка на рис. 7 для (ро = 0.06, а также распределение N(ю) (вставка (в)), ставшее очень узким и высоким. Заметим также, что пик данного распределения не достигает максимальной величины N * N = 14400 из-за конечности времени счета и, как следствие, вычисления частот не в полностью установившемся режиме. Тем не менее данная погрешность очень мала.

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

Рис. 7. Зависимость частот ю^ всех автоколебательных элементов от параметра связи с пассивными элементами dpo. На вставках - соответствующие распределения N

Рис. 8. Моментальные снимки переменной х элементов решетки после интегрирования системы на интервале г е [0; 2 000 000] для четырех значений йро\ 0

состоянию системы. Эти моментальные снимки для четырех значений с!ро представлены на рис. 8. Цвет на рисунках соответствует величине переменной х каждого автоколебательного элемента. Как видно, сначала для ёро = 0 (рис. 8, а) поведение в системе беспорядочное. С ростом ёро до 0.02 (рис. 8, б) появляются небольшие области одного цвета, говорящие о том, что некоторые группы элементов начинают колебаться синхронно, то есть имеет место кластерная синхронизация. Для больших значений ёро, например 0.06 (рис. 8, г), в системе наблюдается режим глобальной синхронизации, выраженный распространяющейся концентрической волной, причем частота данной волны - есть частота глобальной синхронизации в системе. Кроме того, для случая ёро = 0.05 также наблюдается режим синхронизации согласно вычисленным частотам (см. рис. 7). Однако на моментальном снимке (рис. 8, в) мы видим конкуренцию нескольких источников волн. Это объясняется конечностью времени счета, то есть на интервале времени Ь £ [0; 2 000 000] мы все еще имеем дело с процессом установления, который из-за близости частот наблюдаемых источников волн может протекать весьма долго.

Заключение

В данной работе нами рассмотрено влияние пассивных элементов на свойства синхронизации в неоднородных автоколебательных ансамблях с использованием модели Бонхоффера-ван дер Поля.

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

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

Во второй части статьи результаты, полученные ранее, были обобщены на случай двумерного неоднородного автоколебательного ансамбля. Была продемонстрирована возможность получения эффекта синхронизации за счет введения пассивных элементов. О возникновении синхронизации в данном случае можно судить как по изменению частотных характеристик системы, так и по виду устанавливающихся в системе пространственно-временных режимов. Основной результат, полученный в последней части данной статьи, является наиболее важным, особенно в контексте исследования динамики сердца. Из результатов многочисленных биологических экспериментов известно, что количество пассивных клеток в синусном узле сердца может существенно увеличиваться с возрастом человека и достигать 70 %. В свою очередь, согласно полученным результатам, влияние пассивных элементов может приводить, в том числе, к разрушению синхронного режима работы синусного узла сердца. Механизмы возникновения и срыва синхронных режимов в подобных системах, предложенные в данной работе, могут лежать в основе эффектов, наблюдаемых в биологических экспериментах.

Работа выполнена при поддержке Федеральной целевой программы «Научные и научно-педагогические кадры инновационной России» на 2009-2013гг. (контракты № П2018, П15, П2308, 02.740.11.5138, П942, 02.740.11.5188), при поддержке РФФИ (гранты 08-02-92004, 08-02-970049, 10-02-00940) и при поддержке корпорации Intel.

Особая благодарность А.К.Крюкову и М.А.Комарову за полезные обсуждения.

Библиографический список

1. Bub G., Shrier A. and Glass L. Spiral wave generation in heterogeneous excitable media // Phys. Rev. Lett. 2002. Vol. 88. 058101.

2. Bub G., Shrier A. and Glass L. Global organization of dynamics in oscillatory heterogeneous excitable media // Phys. Rev. Lett. 2005. Vol. 94. 028105.

3. Blasius B. and Tonjes R. Quasiregular concentric waves in heterogeneous lattices of coupled oscillators // Phys. Rev. Lett. 2005. Vol. 95. 084101.

4. Kheowan O.-U., Mihaliuk E., Blasius B., Sendina-Nadal I. and Showalter K. Wave mediated synchronization of nonuniform oscillatory media // Phys. Rev. Lett. 2007. Vol. 98. 074101.

5. Hwang S., Yea K. and Lee K.J.Complex-periodic spiral waves in confluent cardiac cell cultures induced by localized inhomogeneities // Phys. Rev. Lett. 2004. Vol. 92. 198103.

6. Ten Tusscher K.H.W.J. and Panfilov A.V. // Europace. 2007. Vol. 9. P. 38.

7. Пиковский А.С., Розенблюм М.Г., Куртс Ю. Синхронизация. Фундаментальное нелинейное явление. М: Техносфера, 2003.

8. Pazo D. andMontbrio E. // Phys. Rev. E. 2006. Vol. 73. 055202(R).

9. Bonhoeffer K.F. Modelle der Nervenerregung // Naturwissenschaften. 1953. Vol. 40. P. 301.

10. Kryukov A.K., Petrov VS., Averyanova L.S., Osipov G.V., Chen W., Drugova O. and Chan C.K. // Chaos. 2008. Vol. 18. 037129.

Нижегородский государственный Поступила в редакцию 10.11.2009

университет им. Н.И. Лобачевского После доработки 25.02.2010

INFLUENCE OF PASSIVE ELEMENTS ON THE SYNCHRONIZATION OF OSCILLATORY ENSEMBLES

VS. Petrov, G.V. Osipov

This paper deals with the influence of the passive elements on the synchronization in the ensembles of coupled non-identical Bonhoeffer-van der Pol oscillators. With a help of numerical experiment it was demonstrated that the introduction of passive elements may lead to both increase and decrease of global synchronization threshold in the system. These results were confirmed analytically using piecewise linear approximation of the Bonhoeffer-van der Pol model. The results obtained for the ensemble of three coupled elements are generalized to the case of two dimensional lattices of non-identical oscillators.

Keywords: Synchronization, oscillatory ensemble, passive element.

Петров Валентин Сергеевич - родился в 1986 году, окончил Нижегородский государственный университет им. Н.И. Лобачевского (2008). Аспирант кафедры теории управления и динамики машин факультета ВМК ННГУ. Область научных интересов: структуры синхронизации в нелинейных системах, моделирование динамики биологических сред. Имеет 12 публикаций (в соавторстве) в том числе 4 журнальных. Лауреат стипендии имени академика Г.А. Разуваева (2009).

603950, Россия, г. Нижний Новгород, пр. Гагарина, 23 Нижегородский государственный университет им. Н.И. Лобачевского E-mail: valentin.s.petrov@gmail.com

Осипов Григорий Владимирович - родился в Горьком (1960), окончил Горь-ковский государственный университет им.Н.И.Лобачевского (1982). В 19861988 годах учился в аспирантуре Научного совета по комплексной проблеме «Кибернетика» (НСК) АН СССР. Защитил диссертацию на соискание ученой степени кандидата физико-математических наук в НСК (1988) по математическому моделированию и доктора физико-математических наук в ННГУ (2004) по радиофизике. С 1988 года работает в ННГУ, с 2007 года заведует кафедрой Теории управления и динамики машин. Один из авторов монографий «Устойчивость, структуры и хаос в нелинейных сетях синхронизации» и «Synchronization in oscillatory networks». Опубликовал 120 научных статей (в том числе 2 обзора) по теории колебаний и волн и математическому моделированию.

603950, Россия, г. Нижний Новгород, пр. Гагарина, 23 Нижегородский государственный университет им. Н.И. Лобачевского E-mail: osipov@vmk.unn.ru