CC BY
51
Список литературы
1. Hodjkin A.L., Haxley A.F. The dual effect of membrane potential on sodium conductance in the giant axon of Loligo.// J. Physiol. 1952. Vol. 116. P. 497-551.
2. Скотт Э. Волны в активных и нелинейных средах с приложением к электронике / Пер. с англ. под ред. Л. А. Островского и М.И. Рабиновича. М.: Сов. радио, 1977. 368 с.
3. Леннинджер А. Биохимия / Пер. с англ. М.: Мир, 1974. 483 с.
4. Зайко Ю.Н. Второе начало термодинамики и неоднородные системы: приложение к космологии и биофизике. // Биомедицинские технологии и радиоэлектроника. 2003. № 2. C.69-72.
5. Румер Ю.Б., Рывкин М.Ш. Термодинамика, статистическая физика и кинетика. М.: Наука, 1972. 400 с.
6. Роу Дж. Теория нелинейных явлений в приборах сверхвысоких частот / Пер. с англ. под ред. З.С. Чернова. М.: Сов. радио. 1969. 615 с.
7. Шредингер Э. Что такое жизнь. С точки зрения физики / Пер. с англ. под ред. А.А. Малиновского и Г.Г. Порошен-ко. М.: Атомиздат, 1972. 88 с.
8. Николс Дж. Г. и др. От нейрона к мозгу / Пер. с англ. под ред. П.М. Балабана и Р. А. Гиниатуллина. М.: УРСС, 2003. 671 с.
УДК 530.182:577.3
СЛОЖНЫЕ КОЛЕБАНИЯ И СИНХРОНИЗАЦИЯ
В ФУНКЦИОНАЛЬНОЙ МОДЕЛИ ВАСКУЛЯРНОГО ДЕРЕВА НЕФРОНОВ
П.А. Щербаков, О.В. Астахов, Д.Э. Постнов*
Саратовский государственный университет E-mail: scherbakovpa@Gmail.com
Предложена математическая модель, на качественном уровне воспроизводящая колебательные процессы в системе авторегуляции почечного кровотока. Модель содержит ансамбль двумодовых осцилляторов, взаимодействующих посредством двух типов (каналов) связи, различающихся как геометрией (древовидная структура и локальное взаимодействие), так и характером воздействия на индивидуальный осциллятор (модуляция энергонесущего параметра и процесс диффузионного типа). Выявлены основные режимы функционирования такой системы и взаимопереходы между ними.
Ключевые слова: синхронизация, функциональная модель, васкулярное дерево, нефрон, гемодинамика, мультистабильность, кластер.
COMPLEX WAVEFORMS AND SYNCHRONIZATION IN FUNCTIONAL MODEL OF VASCULAR NEPHRON TREE
P.A. Scherbakov, O.V. Astakhov, D.E. Postnov
We suggest functional model that qualitatively describes oscillatory processes in renal autoregulation. Our model consists of ensemble of two-mode oscillators that are coupled by means of two different pathways. The above coupling pathways count both the geometry of ensemble (tree-like structure or local interaction) and the specific action of individual oscillator (energy distribution netrwork or diffusive coupling). We study the typical operating regimes of suggested model as well as transitions between them.
Key words: synchronization, functional model, vascular tree, nefron, hemodynamics, multistability, cluster.
Введение
Исследование поведения взаимодействующих автоколебательных систем - одна из основных задач теории колебаний и нелинейной динамики, которая во многом сводится к выявлению и изучению разнообразных проявлений общего и фундаментального явления - синхронизации [1-3]. Диапазон приложений таких исследований весьма широк. Например, в радиофизике - это ансамбли связанных осцилляторов [4], в нейрофизиологии - коллективная динамика нейронов-пейсмейкеров [5].
Замечательно, что при всем разнообразии проявлений механизмы синхронизации универсальны и могут поэтому изучаться на примере относительно небольшого числа модельных систем. При этом расширение круга изучаемых конкретных задач связано с:
- усложнением индивидуальной динамики взаимодействующих автоколебательных систем;
- изучением различных типов их взаимодействия;
© П.А. Щербаков, О.В. Астахов, Д.З. Постнов, 2009
- увеличением числа взаимодействующих осцилляторов;
- учетом способа их организации в ансамбль (геометрия связей).
По каждому из упомянутых направлений за последние десятилетия имеется значительный прогресс. Так, установлены основные закономерности синхронизации для хаотических [6, 7] и индуцированных шумом колебаний. Исследуются различные типы связи и соответствующие им характеристики синхронных режимов ансамблей осцилляторов [8-10]. Изучаются особенности кооперативной динамики осцилляторов при различных способах их организации в ансамбль [11, 12].
Вместе с тем бесконечное разнообразие окружающего мира предлагает все новые задачи со своей спецификой как свойствами самих взаимодействующих осцилляторов, так и характером их взаимодействия. Одна из таких задач возникла в результате исследований физиологами закономерностей авторегуляции почечного кровотока.
Почки, будучи одним из важнейших органов человека, предназначены для поддержания постоянства объема и состава внеклеточной жидкости, омывающей клетки организма, тем самым обеспечивая оптимальные условия их жизнедеятельности [13]. Они выводят из организма избыток воды и растворенных в ней веществ, а также выполняют различные метаболические функции. Почка человека содержит около миллиона нефро-нов, каждый из которых может считаться мельчайшей функциональной единицей почки [14, 15] и способен самостоятельно выполнять ее главную функцию: регулировать химический и ионный состав плазмы крови посредством ее фильтрации. Одним из важных достижений в области ренальной физиологии было открытие автоколебательной динамики в функционировании нефронов.
Авторегуляция почечного кровотока:
основные механизмы и экспериментальные свидетельства
автоколебательной динамики
Каждый нефрон (рис. 1, a) содержит так называемый клубочек, играющий роль фильтра, с кровеносными сосудами, приносящими и отводящими поток крови (аффе-
рентная и эфферентная артериолы). Фильтрат не удаляется из нефрона немедленно, а проходит процедуру регуляции химического состава в процессе продвижения по проксимальному извитому канальцу и так называемой петле Генле. Перед тем как перейти в дистальный каналец, петля Генле примыкает непосредственно к клубочку (macula densa). Здесь имеет место организованная природой обратная связь: в зависимости от химического состава фильтрата меняется тонус клеток гладкой мускулатуры, окружающих стенку афферентной артериолы, и ее диаметр изменяется: артериола сжимается при чрезмерно интенсивной фильтрации и расслабляется при ее недостаточной интенсивности. Такой механизм, наряду с заметным временем прохождения фильтратом петли Генле (порядка 13 с), создает предпосылки для автоколебательной динамики. И действительно в экспериментах на нефронах крысы in vivo P.P. Le-ysec и N.H. Holstein-Rathlow наблюдали автоколебания проксимального гидростатического давления с периодом 20-40 с [16, 17]. Была экспериментально доказана важная роль упомянутых выше факторов: канальцевогломерулярной обратной связи (КГОС), наличия существенной задержки при продвижении фильтрата по петле Генле, а также миогенного отклика афферентной артериолы [18-20].
Колебания той же частоты, но с фазовым сдвигом, наблюдались и для концентрации ионов натрия и хлора в дистальном канальце нефрона. Согласно результатам измерений, высокая концентрация ионов хлора соответствует минимуму давления в канальце [21]. Спектральный анализ записи артериального давления и его сопоставление с колебательными свойствами отдельных нефро-нов подтверждают предположение о том, что наблюдаемая колебательная динамика является результатом внутренних свойств системы, а не привнесена извне. Была выявлена связь гипертонии с характеристиками колебаний в нефронах. В частности, было экспериментально зафиксировано наличие нерегулярных (хаотических) колебаний проксимального давления для крыс с повышенным артериальным давлением [22].
Дистальная трубка
Афферентная артериола Клубочек
Петля Генле
Рис. 1. Схематичное представление основных частей нефрона и путей их взаимодействия. В сегментах 1 и 2 удаляется вода, в сегментах 3-5 регулируется концентрация ионов натрия и хлора (а). Схематичное представление ветвления кровеносных сосудов почки от единственной почечной артерии (внизу) до васкулярного дерева
из 10-20 нефронов (вверху) (б)
Почка как ансамбль осцилляторов: геометрия и характеристики связей
Задача снабжения кровью каждого из нефронов «решена» природой типичным для живых систем образом: единственная почечная артерия неоднократно ветвится на более мелкие сосуды, распределенные по мозговому веществу и коре почки, как это схематично изображено на рис. 1, б. В результате образуется 50-100 тысяч отростков артерий, каждый из которых (уже не ветвясь) снабжает кровью по 10-20 нефронов, соединенных с ним посредством афферентных артериол. Такую структуру ветвления сосудов называют васкулярным деревом (vascular tree).
Для краткости описания далее под васкулярным деревом нефронов будет пониматься совокупность нефронов и общего кровеносного сосуда, к которому они присоединены. Следует заметить, что в таких деревьях около 50% нефронов расположены парами (имеют общую артериолу) и порядка 5% - триплетами.
Очевидно, что авторегуляция интенсивности кровотока на входе в отдельно взятый нефрон влияет на давление крови в ближайшей к нему точке ветвления. Таким образом, изменение состояния одного нефрона заставляет соседние нефроны реагировать на изменившиеся условия в процессе поддержания нужного режима фильтрации. Результатом
является «конкуренция» нефронов за поток крови: его рост на входе одного из нефронов приводит к уменьшению потока соседнего и наоборот.
Кроме того, электрические потенциалы, управляющие сечением артериолы на входе в нефрон (так регулируется поток втекающей крови), способны распространяться вдоль стенок сосуда и достигать артериол близкорасположенных нефронов, вмешиваясь в их работу [23, 24]. Такое взаимодействие приводит к обратному эффекту: уменьшение просвета афферентной артериолы нефрона индуцирует тот же (хотя и ослабленный) эффект для соседей.
Таким образом, неизбежность взаимодействия нефронов следует уже из анатомических особенностей строения почки. Подтверждается это и результатами экспериментов, а именно удается зафиксировать режимы синхронизации нефронов, принадлежащих к одному дереву [25]. При этом наблюдаются как синфазные, так и противофазные режимы колебаний.
Функциональный подход в моделировании
На данный момент известны количественные математические модели одиночных [26-29], а также парных нефронов [30, 31]. Эти модели неплохо воспроизводят такие характеристики, как частота колебаний, а
а
также уровень максимального и минимального давления в проксимальном канальце. Однако работ, связанных с количественным моделированием процессов в рамках васкулярного дерева нефронов, крайне мало [11, 24, 32, 33].
Во многом это связано со сложностями, возникающими при численном моделировании больших ансамблей количественных моделей нефронов в структуре типа васкулярного дерева. Такая математическая модель имеет сотни динамических переменных и десятки управляющих параметров, далеко не все из которых имеют адекватные экспериментальные оценки. В то же время представляет несомненный интерес исследование хотя бы общих свойств такого ансамбля неф-ронов, как нелинейная автоколебательная система и выявление типичных режимов ее функционирования в зависимости от соотношения управляющих параметров.
В свете всего вышеперечисленного привлекательным и адекватным является использование функциональных моделей, а именно как сами нефроны, так и способ их организации в ансамбль представляются упрощенными модельными уравнениями, воспроизводящими, однако, наиболее существенные свойства прототипа. В нашем случае такими свойствами можно считать:
1) наличие двух источников колебательной динамики в одиночном нефроне; при этом «медленная мода» связана с автоколебаниями вследствие наличия задержки в петле Генле и КГОС, тогда как быстрая мода отображает в 3-5 раз более быстрые колебания диаметра афферентной артериолы при регуляции сосудистого тонуса, которые носят затухающий характер;
2) характер зависимости автоколебаний в нефроне от артериального давления: при слишком низком или слишком высоком давлении в афферентной артериоле автоколебания не наблюдаются;
3) древовидная структура связи, которая в пределе упрощения может быть разложена на минимальные ансамбли 2 и 3 нефронов, подсоединенных к общему источнику (кровеносному сосуду), либо на группу, подсоединенную последовательно к общему «стволу» через некоторые промежутки;
4) дуальный характер взаимодействия, обусловленный упомянутыми выше механизмами.
В нашей работе мы предлагаем имитационную модель нефрона, а также технологию численного моделирования, позволяющую воспроизвести различные варианты геометрии васкулярного дерева. Как показывают результаты вычислительного эксперимента, такая модель воспроизводит на качественном уровне динамику нефрона в части быстрой и медленной колебательных мод.
1. Имитационная модель нефрона в виде двумодоеого осциллятора
Наиболее простой вид функциональной модели нефрона был предложен в работе [34]. Эта модель представляет собой двумодовый осциллятор с соотношением частот медленной и быстрой моды от 1:1 до 1: 6. Модель состоит из двух связанных осцилляторов: одного, демонстрирующего автоколебательную динамику (медленная подсистема), и второго - демпфированного (быстрая подсистема). Каждый из осцилляторов представлен системой из двух дифференциальных уравнений первого порядка.
В работе [34] выбор вида нелинейностей быстрой подсистемы основан на анализе уравнений количественной модели Барфреда одиночного нефрона [26], в то время как уравнения для медленной подсистемы записаны из самых общих соображений (наличие автоколебательной динамики) в виде осциллятора ван дер Поля [35]. Как известно, в этой классической модели теории колебаний отсутствует параметр, явно описывающий приток энергии в систему. Для нефрона роль такого энергонесущего параметра играет артериальное давление крови на входе в неф-рон [32], и учет его в модели представляется принципиально важным.
Для целей нашего исследования была предложена модель в следующем виде:
* = у, (1)
у = Р - у(Я + N'(х)) - NN(х) - х + су, (2) и = V, (3)
У = -Бу - и(1 + веи) + Г(х, и); (4)
где
х) = х2 -10х + 20,
N (х) = х3 /3 - 5 х2 + 20 х,
^(х,и) = А 1апЬ(^и (и + йи ))(1 + т).
В качестве медленной подсистемы (1), (2) был использован осциллятор с нелинейностью К-типа из [36]. В гемодинамической интерпретации параметр Я есть безразмерный аналог сопротивления входной артерио-лы, а Р - аналог давления на входе артерио-лы. Выбранная модель медленной подсистемы относится к широко известному классу осцилляторов с нелинейностью К-типа, но отличается наличием параметра, эквивалентного напряжению питания в радиофизической интерпретации. Данная модель анализировалась в [33, 36]. В частности, было показано, что автоколебания возможны в ограниченной области значений Е. Это хорошо согласуется с результатами по моделям неф-ронов, где для слишком низкого и слишком высокого кровяного давления автоколебаний не наблюдается.
В качестве быстрой подсистемы (3), (4) нашей модели был оставлен нелинейный, параметрически возбуждаемый демпфированный осциллятор, использованный в [36], с масштабными коэффициентами ки и ёи для адаптации к выбранной медленной подсистеме. Связь двух подсистем осуществляется посредством нелинейной функции Р(х,и) в быстрой подсистеме и коэффициента с обратной связи в медленной подсистеме.
Типичный набор значений управляющих параметров включал: Я = 0.15, с =13.2, в=0.0009, А = 0.474, у =11.09, Б=0.1061, ки= = 0.4, ёи = -5.
На рис. 2 приведена однопараметрическая диаграмма, характеризующая наличие автоколебаний в модели (1)-(4) в зависимости от значения параметра Р. Как можно видеть, они отсутствуют при значениях Р менее 4.55 и больших, чем 9.3. Наличие быстрой подсистемы, обеспечивающей, по сути, обратную связь, приводит к эффекту гистерезиса: в области 4.55-6.41 режим автоколебаний и режим устойчивой точки сосуществуют.
40
20
у о
-20
Рис. 2. Зависимость амплитуды колебаний модельной системы от параметра Р. При значениях меньших 4.55 и больших 9.3 автоколебания в осцилляторе отсутствуют. В области значений параметра Р между 6.41 (пунктирная линия) и 9.3 автоколебательный режим представляет единственный аттрактор системы, тогда как в области 4.55-6.41 режим автоколебаний и режим устойчивой точки сосуществуют
Помимо самого факта наличия автоколебательной динамики, главное свойство модели (1)-(4) - это наличие двух различных временных масштабов, один из которых определяется автоколебательной подсистемой (1)-(2), а второй - колебаниями в нелинейном осцилляторе (3)-(4). Несмотря на наличие в этой подсистеме постоянного затухания (параметр Б), сложный характер воздействия со стороны медленной подсистемы, задаваемый Р(х,и), порождает разнообразные колебательные режимы. В [34] упоминается наличие квазипериодических режимов, а также детально обсуждается раздельная синхронизация двух пар базовых частот хаотических автоколебаний. Таким образом, модель (1)-(4) демонстрирует две моды колебаний на различных частотах, которые могут находиться как в рациональном соотношении (резонанс), так и быть рассинхронизованы.
Ранее в [31] был обнаружен эффект рассинхронизации быстрой и медленной мод хаотических колебаний в количественной модели нефрона, где каждое из характерных времен имеет ясную физиологическую интерпретацию (время прохождения фильтратом петли Генле и квазипериод колебаний
диаметра афферентой артериолы). Нами была выявлена аналогичная перестройка хаотического автоколебательного режима модели (1)-(4). На рис. 3, а приведена характерная временная реализация колебаний по переменной у. Низкочастотная мода колебаний представлена последовательностью выбросов положительной и отрицательной полярности. Высокочастотная мода колебаний отражена заполнением интервалов между выбросами, причем форма этого заполнения (амплитуда и количество пиков) не повторяется от квазипериода к квазипериоду. Для количественной характеристики соотношения между частотами мод рассчитывалось число вращения г как среднее по реализации число колебаний
быстрой подсистемы на один период колебаний медленной подсистемы. График изменения г при вариации Р представлен на рис. 3, б. Несмотря на то, что интервал значений Р на графике соответствует хаотическому режиму, величины г выявляют две качественно различные ситуации, а именно при Р <5.293 число вращения изменяется с изменением Р, отражая факт независимого изменения двух временных масштабов хаотических колебаний. При Р > 5.293 число вращения г равно 1/6 в пределах точности измерений: на каждый квазипериод низкочастотной моды приходится ровно 6 квазипериодов высокочастотной моды, в то время как колебательный режим в целом остается хаотическим.
Рис. 3. Пример временной реализации для случая хаоса с рассинхронизованными модами Р = 5.288 (а). Зависимость числа вращения г от параметра Р. При значении параметра Р > 5.293 наблюдается хаотический режим с синхронизованными модами, число вращения г = 1/6. При Р < 5.293 моды хаотического режима рассинхронизованы, число вращения меняется при вариации Р (б)
На рис. 4 такая перестройка хаотического режима охарактеризована более детально. Приведены фазовые проекции и спектры колебаний для двух значений Р: Р = 5.3 (левая колонка) и Р = 5.288 (правая колонка). Различие в характеристиках очевидно: при сохранении общей структуры фазового портрета
рассинхронизация мод отображается наличием дополнительных петель траектории на панели (г), выбивающихся из «общего русла» шестиоборотной структуры на панели (а). Сравнение панелей (б) и (д) позволяет увидеть, что в окрестности утолщения пучка траекторий (левая нижняя четверть рисунка)
а
при рассинхронизации мод могут образовываться дополнительные петли, наличие которых и определяет мгновенное значение г. Диагностика с помощью усредненных спектров мощности по переменной х дает интересные результаты: спектр на панели (в) имеет все признаки двухчастотного колебатель-
ного режима (г), хорошо видны пики, соответствующие низкочастотной и высокочастотной моде колебаний, тогда как спектр на панели (в) более похож на типичный развитый спектр одномодового хаоса. При этом низкочастотный пик оказывается размытым и практически не наблюдается.
10.0
У 2.0 -
-6 0 —
5.0
10.С
д 10.0
У 2.0
-6 0
10.0
(I)
(0
Рис. 4. Сравнение характеристик хаотического режима с синхронизованными (слева, Р = 5.3) и рассинхронизованными (справа, Р = 5.288) модами. Проекции фазового пространства на плоскость х:у (а, г); то же в укрупненном масштабе (б, д); усредненные спектры мощности колебаний по переменной х (в, е)
б
Этот факт согласуется с поведением фазовых траекторий: различное количество
мелкомасштабных петель траектории приводит к большой дисперсии времени в левую часть фазового портрета на панелях (а) и (г), что и размывает проявление соответствующего временного масштаба в спектре.
Таким образом, наблюдаемый эффект рассинхронизации мод хаотического режима не сводится ни к классическому случаю синхронизации двух частот «на фоне» хаотических автоколебаний, что было бы близко к активно изучавшимся в последние десятилетия механизмам синхронизации хаоса, ни к
типичной эволюции характеристик хаотического аттрактора по мере увеличения параметра надкритичности (этот вопрос обсуждался в [34]). Для целей нашего исследования важно, что поведение имитационной модели (1)-(4) воспроизводит характерные черты динамики количественных моделей неф-ронов и тем самым обосновывает ее применение для целей моделирования более крупных структур.
2. Моделирование структуры связей васкулярного дерева
Как уже упоминалось выше, взаимовлияние нефронов в васкулярном дереве обусловлено как модуляцией давления в общем кровеносном сосуде (гемодинамическая
связь), так и распространением электрического потенциала по стенкам артериол (мио-генная связь). Для математической модели двух или трёх нефронов указанное взаимодействие нетрудно описать в явном виде путём соответствующей модификации уравнений. Однако для целей моделирования более крупных ансамблей нефронов и, особенно, с учетом их сложной геометрии более гибким и эффективным представляется описанный ниже метод, основанный на дискретном представлении неоднородной структуры среды [37]. А именно пространственная структура моделируемого ансамбля представлена в виде набора элементов с различными свойствами и различным характером локального взаимодействия с ближайшими соседями.
На практике такая матрица элементов задается в виде графического файла, цвет каждого пикселя которого кодирует свойства элемента. На рис. 5 каждой клетке соответствует своя система уравнений.
(0) - Пустая клетка. Правые части уравнений и сами переменные тождественно равны нулю, взаимодействие отсутствует.
(1) - Имитация большого кровеносного сосуда, давление Р в котором постоянно и равно р:
V- - V
р = Р, р = ° = ^-^—.
Здесь уравнение для потенциала V на качественном уровне описывает пассивное распространение электрического сигнала по стен-
Рис. 5. Пример моделирования двух связанных двумодовых осцилляторов при помощи пиксельной матрицы
кам сосуда, К* - погонное сопротивление, ц -параметр, пропорциональный погонной емкости. В случае соседства с элементами типов 0 или 3, для которых величина V не определена, принимается V;. Индекс - в этом и последующих описаниях элементов принимает значения - = 1, г, ¿, Ь, которые указывают на элемент, расположенный слева, справа, сверху и снизу соответственно.
(2) - Имитация малого кровеносного сосуда с меняющимся давлением:
& = У р-р, ^ = У ^,
^ Кс ^ К*
где Яс - погонное гемодинамическое сопротивление, ц - параметр, пропорциональный эластичности сосуда.
(3) - Медленная подсистема (имитация КГОС-осциллятора нефрона). Значения р; и VI не равны нулю только в случае соседнего элемента типа 1, 2 или типа 4, соответственно.
х = у,
у = р; - у(Я + N (х)) - NN(х) - х + .
(4) - Быстрая подсистема (имитация миогенного отклика артериолы нефрона):
и = V,
V - V
V = -Бу - и(1 + ве ) + Г(х,и) + —----.
В отличие от предыдущих элементов, индекс при Vi указывает на элементы, соседние к смежному элементу типа 3.
В процессе компьютерного моделирования все элементы описанной выше матрицы перебирались при выполнении каждого шага численного интегрирования. При этом количество взаимодействующих двумодовых осцилляторов и геометрия их подключения к общему кровеносному сосуду задавались просто соответствующей закраской пикселей графического файла, что позволяло моделировать самые разнообразные случаи. Ниже описаны результаты, полученные для базовых (простейших) конфигураций, исследование которых необходимо для лучшего понимания нелинейных механизмов, определяющих динамику такой системы.
3. Синхронные колебательные режимы в случае простейших структур связи
В почке млекопитающих порядка 50% нефронов расположены парами (делят общий участок артериолы), а 5% - соединены в триплеты. По этой причине вопрос о динамике таких структур представляет очевидный интерес.
В рамках решаемой нами задачи имеются две особенности, определяющие специфику проявления синхронизации и отличия от хорошо исследованных случаев. Каждая из взаимодействующих автоколебательных систем имеет две моды колебаний, и связь осуществляется по двум каналам, по-разному воздействующим на осцилляторы. В рамках данного раздела проиллюстрируем те особенности, которые возникают при симметричном соединении двух осцилляторов. При их моделировании матрицей элементов (см. рис. 5) кровеносный сосуд от точки с постоянным значением параметра Р до точки ветвления был представлен 8 элементами матрицы, участок от точки ветвления до осциллятора - тремя элементами и один элемент описывал точку ветвления.
При наличии только гемодинамической связи через переменную р (что задается условием 1/Я* = 0) имеет место тенденция к противофазной синхронизации медленной моды колебаний.
Рассмотрим два связанных осциллятора при Яс = 0.001, £ = 40. Как можно видеть на рис. 6, в этом случае устанавливается синхронный режим колебаний с разбегом фазы в половину периода. Подобные результаты были получены в [36] для связанных по питанию одномодовых осцилляторов. В нашем случае к колебаниям медленной моды «привязаны» колебания быстрой моды, хорошо различимые на рисунке. Их наличие и определяет сложную геометрию фазовой проекции.
Следует отметить, что сдвиг фаз колебаний быстрой моды (его легко диагностировать по фазовому портрету, см. рис. 6, б) не постоянен в пределах общего периода. Это подтверждает тот факт, что синхронизация осуществляется в результате взаимодействия осцилляторов по медленной моде.
В условиях только миогенной связи посредством переменной V, а также смешанной связи по обоим каналам наблюдается гораздо более сложная картина. В зависимости от выбора величин Я* и Яс наблюдается различное количество сосуществующих синхронных режимов колебаний, различающихся сдвигом фаз - имеет место фазовая мультистабильность. На рис. 7 приведены временные реализации и фазовые проекции для трех таких режимов. Как можно видеть, сдвиг фаз Дф по общему периоду колебаний (то есть по периоду медленной моды) не соответствует ни синфазной, ни противофазной синхронизации и составляет 0.4л, 0.54л, 0.7п. Однако на графиках временных реализаций переменных хорошо видно, что во всех случаях имеется несколько (точнее, 4, 3 и 2) квазипериодов колебаний быстрой моды, минимумы которых совпадают для обоих переменных. Данное наблюдение позволяет диагностировать механизм образования множественных синхронных режимов, опираясь на работы по механизму фазовой мультистабильности [38, 39, 40].
Как уже упоминалось, миогенная связь носит диффузионный (выравнивающий характер), чем обусловлена тенденция к синфазной синхронизации по переменной V, то есть по колебаниям быстрой моды. Однако одному и тому же сдвигу фаз колебаний
X2
*2
і
і
Рис. 6. Противофазная синхронизация двух связанных двумодовых осцилляторов при Р = 8.6, Яс = 0.001, 1/Ял = 0: а - проекция фазового пространства на плоскость х1 : х2; б - проекция фазового пространства на плоскость у1 : у2; в - временные реализации х1 : х2; г - временные реализации у1 : у2
б
а
в
г
Рис. 7. Временные реализации (слева) и фазовые проекции (справа) трех сосуществующих синхронных режимов при Р = 6.0, Я* = 16.5, ц = 1, Яс = 0.01, £ = 40. Сдвиг фаз колебаний составляет: а - Дф = 0.4п, б - Дф = 0.54п, в - Дф = 0.7п. Непрерывная и пунктирная линии соответствуют переменным х1 и х2. Фазовые проекции справа даны по переменным х1 : х2 в пределах [-5; 15] по обоим направлениям
а
б
в
быстрой моды соответствует несколько вариантов сдвига фаз медленной моды колебаний. При подходящих значениях параметра связи сразу несколько вариантов синхронных режимов, синфазных по быстрой моде колебаний, но отличающихся сдвигом фаз медленной моды могут быть одновременно устойчивы. Данный эффект был продемонстрирован для двумодовых колебаний в [38], его особенности при одновременной связи по нескольким переменным обсуждались в [39]. В [40] результаты были обобщены как «субгармонический механизм фазовой мультистабильности». В отличие от упомянутых выше работ, в нашем случае колебания быстрой моды не непрерывны в пределах общего периода: как можно видеть на рис. 7, амплитуда и частота быстрой моды значительно отличаются для положительной и отрицательной полуволны медленной моды колебаний. С одной стороны, это делает невозможным синхронизацию на определенных сдвигах фаз, а с другой - потенциально способствует формированию дополнительных синхронных режимов. Полная картина фазовой мультистабильности в исследуемом случае весьма сложна. В частности, на рисунке приведены не все обнаруженные при данных параметрах синхронные режимы. Однако механизм их формирования представляется вполне ясным и соответствует описанному выше.
Увеличение силы взаимодействия по каналу гемодинамической связи (рост параметра Яс) усиливает тенденцию к синхронизации в противофазе и, таким образом, ослабляет устойчивость в первую очередь тех режимов, которые характеризуются наименьшим сдвигом фаз (см. рис. 7, а), тогда как режимы, более близкие к противофазному (см. рис. 7, в), могут становиться более устойчивыми.
Итак, динамика двух симметрично связанных двумодовых осцилляторов характеризуется значительным разнообразием синхронных режимов: чисто гемодинамической связи соответствует противофазная синхронизация по медленной моде колебаний, а значит и в целом. Чисто миогенной связи соответствует синхронизация в фазе по бы-
строй моде колебаний, но наличие медленной моды порождает несколько сдвигов фазы колебаний в целом, отвечающих данному условию. В результате одновременное действие двух каналов связи приводит к наличию набора синхронных режимов с различным сдвигом фаз в зависимости от силы взаимодействия по каждому из каналов. Полная картина, включающая все реализующиеся синхронные режимы и области их устойчивости на плоскости управляющих параметров, представляется весьма сложной и требует отдельного исследования.
4. Кластерная генерация и индуцированный связью сдвиг частоты колебаний
Структура разветвленного васкулярного дерева нефронов, на функциональное моделирование которой нацелена наша модель, может быть разложена на подансамбли с более упорядоченной геометрией связи. В предыдущем разделе анализировался простейший случай, когда осцилляторы ансамбля равноудалены от точки с постоянной величиной Р. Другой важной ситуации соответствует распределение нефронов вдоль «ствола» дерева и описывается моделью, в которой ряд нефронов (в нашей модели - двумодовых осцилляторов) последовательно присоединен к общему кровеносному сосуду, давление в котором падает по мере удаления от «основания» (точки с постоянным Р). Такой ансамбль имеет структуру, схематически показанную на рис. 8 и рассматривался в работах [32-34, 36].
Рис. 8. Структура связи типа «цепь». Цифрами отмечены номера осцилляторов. В пиксельном шаблоне расстояние от большого кровеносного сосуда до первой точки ветвления малого кровеносного сосуда составляет 4 пикселя, между точками ветвления малого кровеносного сосуда -один пиксель плюс один на точку ветвления, длина арте-риолы каждого элемента составляет 3 пикселя
В частности, был обнаружен эффект так называемой кластерной генерации, когда в силу особенностей связи в автоколебательном режиме оказываются не все элементы цепочки, а лишь их группа (кластер). Этот эффект является следствием двух факторов. Во-первых, область существования автоколебательного режима для отдельно взятого осциллятора ограничена по величине р: при слишком малых и слишком больших его значениях автоколебаний нет. Во-вторых, средняя величина р будет уменьшаться по мере удаления от начальной точки с постоянным Р. Таким образом, положение рабочей точки осцилляторов будет меняться при движении вдоль цепи, и вполне вероятна ситуация, когда лишь часть их них генерирует колебания,
тогда как остальные находятся в передемп-фированном режиме.
На рис. 9 приведены временные реализации по переменной х для 10-ти осцилляторов цепочки при одном (а) и при двух (б) каналах связи. Как можно видеть, в режиме автоколебаний находятся осцилляторы с четвертого по девятый в колонке (а) и со второго по восьмой - в колонке (б). Для осцилляторов, находящихся ближе к началу цепочки, получаемое ими значение р избыточно, в то время как для осцилляторов, находящихся ближе к концу цепочки, - оно недостаточно для возбуждения автоколебаний. На рисунке можно также различить колебания малой амплитуды в осцилляторах, находящихся вне осцилляторного кластера.
i 1 (УУП^УУУ)
гптггттгп п/уууууу1
ттттпптгп [ууууууу!
КШ1ШГ11 [УУУУУУУ1
ггггттгшж! гуууууун
ггттптггттп!
12 x¡ 4
-4
штшжш
1600 1650 1700 1750 1800 1600 1650 1700 1750 1800
t t
Рис. 9. Временные реализации для группы из 10-ти осцилляторов цепочки: с 3 по 12 элемент: a - случай одного канала связи P = 13.0, Rc = 0.001, £ = 40, 1/Rd = 0, ц = 1; в режиме автоколебаний находятся только шесть осцилляторов, колебания во всех шести осцилляторах не синхронизованы, хорошо выражена быстрая мода колебаний; б - случай двух каналов связи, P = 13.0, Rc = 0.001, £ = 40, Rd = 0.25, ц = 1. В режиме автоколебаний находятся шесть осцилляторов, однако за счет модуляции питания в центральной ветви цепочки наблюдаются колебания небольшой амплитуды в десятом осцилляторе; четыре осциллятора совершают синхронные колебания в фазе; хорошо выражена тенденция к синфазному режиму в остальных осцилляторах кластера
б
а
Они обусловлены модуляцией текущей величины р теми из осцилляторов, которые находятся в режиме автоколебаний. В случае только гемодинамической связи частоты колебаний в осцилляторах, составляющих кластер, заметно разнятся, что хорошо видно на рис. 9, а. Включение второго канала связи приводит к понижению частоты автоколебаний осцилляторов, составляющих кластер, примерно в два раза. Кроме того, на рис. 9, б наблюдается хорошо выраженная тенденция
к синфазной динамике. Эти особенности могут быть объяснены следующим образом: синфазные колебания в соседних осцилляторах кластера приводят к более глубокой модуляции величины р. Чрезмерное уменьшение (равно как и увеличение) р смещает каждый из осцилляторов в режим, близкий к границе области генерации. Согласно [36], это приводит к существенному росту периода колебаний.
8.4
1400
1600
1800
2000
а1
61
Рис. 10. Временные реализации колебаний величины р в главной ветви цепочки в середине осцилляторного кластера: а - случай только гемодинамической связи при Р = 13, Кс = 0.001, £ = 40, 1/К = 0, ц = 1; б - случай двух каналов связи при Р = 13, Кс = 0.001, £ = 40, К = 0.25, ц = 1; а1, 61 - фрагменты временных реализаций, выделенные пунктиром на верхних панелях
а
б
і
Дополнительную информацию о поведении системы дает рис. 10, где приведены графики изменения p в центральной ветви цепочки. Как можно видеть, в случае чисто гемодинамической связи (панели a и а1) колебания p носят нерегулярный характер, что объясняется наложением колебаний с различными частотами. Кроме того, заметны слабые модуляции с частотой колебаний быстрой подсистемы. При добавлении миоген-ной связи (панели б-б1) модуляция в центральной ветви цепочки становится более регулярной и глубокой, что увеличивает амплитуду вынужденных колебаний в осцилляторах, граничащих с кластером. Среднее значение p может быть оценено как 7.85 для случая a и 6.9 для случая б. Это говорит о том, что синфазный режим колебаний в данном случае энергетически менее выгоден, чем противофазный режим или режим со сдвигом фазы. В терминах исходной задачи о динамике ансамблей нефронов это означает, что при наличии заметной миогенной связи кровяное давление быстрее спадает по мере удаления от крупных сосудов, что уменьшает предельно достижимую скорость фильтрации. За счет описанных эффектов наблюдается сдвиг осцилляторного кластера в сторону начала цепи при введении второго канала связи.
Выводы
В данной работе нами была предложена модель двумодового осциллятора, на качественном уровне воспроизводящая особенности колебаний в модели одиночного нефро-на. Как было показано, такая модель демонстрирует режимы с различным соотношением частот быстрых и медленных колебаний, а также позволяет наблюдать перестройку хаотического колебательного режима, связанного с рассинхронизацией мод, которая ранее была обнаружена в динамике количественной модели нефрона.
Для построения функциональной модели васкулярного дерева нефронов был предложен способ задания структуры связей в ансамбле двумодовых осцилляторов путем их дискретного представления посредством пиксельной матрицы.
При исследовании динамики двух симметрично связанных двумодовых осцилляторов было обнаружено, что миогенная связь приводит к формированию набора сосуществующих синхронных режимов колебаний, что обусловлено тенденцией к синхронизации в фазе на быстрой моде колебаний. Увеличение относительного вклада гемодинами-ческой связи делает неустойчивыми синхронные режимы с наименьшей величиной сдвига фаз и, в пределе только гемодинами-ческой связи, приводит к противофазной синхронизации на медленной моде колебаний.
Для ансамбля двумодовых осцилляторов в виде цепочки был продемонстрирован эффект кластерной генерации, что согласуется с известными результатами для ансамблей одномодовых осцилляторов [33, 36]. Новым в данной работе является учет влияния мио-генной связи. Как установлено, ее действие приводит к синхронизации в фазе элементов кластера, что усиливает модуляцию переменной p, которая по смыслу соответствует давлению крови в «стволе» васкулярного дерева нефронов. Эта модуляция, в свою очередь, существенно понижает частоту генерации элементов кластера и приводит к сдвигу его в сторону «основания» дерева с постоянным p = P.
Обнаруженные нами различия в характере изменения p (см. рис. 10) интересным образом коррелируют с некоторыми работами физиологов по моделированию почечного кровотока. При обсуждении механизмов его авторегуляции некоторыми физиологами (D.J. Marsh, США) высказывалось предположение, что хаотизация колебаний в нефро-не вызвана, скорее, их взаимодействием в пределах васкулярного дерева, нежели особенностями индивидуальной динамики. Заметим, что в эксперименте измерения проводятся на нефронах в составе васкулярного дерева (а не изолированных). По этой причине физиологические исследования пока не могут подтвердить либо опровергнуть данную гипотезу. Наши результаты показывают, что при слабой миогенной связи колебания имеют весьма сложную форму, что вызвано различием в режиме элементов осциллятор-ного кластера.
В целом результаты данной работы позволяют продвинуться в двух направлениях: первое - это построение функциональных моделей больших ансамблей нефронов, что является все более актуальным по мере появления экспериментальных данных по пространственным кластерам синхронной активности нефронов, а второе направление -это исследование роли пространственной асимметрии в динамике относительно небольших (10-20 нефронов) ансамблей, что облегчается применением описанной в данной работе методики задания геометрии связей при помощи пиксельной матрицы.
Авторы благодарят коллег O. Sosnivtseva, E. Mosekilde (Denmark Technical University), а также А.Н. Павлова (Саратовский госунивер-ситет) за полезные дискуссии в ходе выполнения работы.
Работа выполнена при финансовой поддержке РФФИ (грант № 090201049).
Список литературы
1. LandaP.S., RosenblumM.G. Synchronization and chaotiza-tion of oscillations in coupled self-oscillating systems // Appl. Mech. Rev. 1993. Vol.46. P.414-426.
2. Pikovsky A., Rosenblum M., Kurths J. Synchronization: A Universal Concept in Nonlinear Sciences. Cambridge University Press, 2001. 411 p.
3. Balanov A., Janson N., Postnov D., Sosnovtseva O. Synchronization From Simple to Complex. Springer, 2009. 406 p.
4. Hadley P., Beasley M.R., Wiesenfeld K. Phase locking of Josephson junction series arrays // Phys. Rev. B. 1988. Vol.38. P.8712.
5. Абарбанель Г.Д.И., Рабинович М.И., Сильверстон А., Баженов М.В., Хуэрта Р., Сущик М.М., Рубчинский Л.Л. Синхронизация в нейронных ансамблях // УФН. 1999. Т.166, №4. С.365-390.
6. Анищенко В. С., Астахов В.В., Вадивасова Т.Е., Нейман А.Б., Стрелкова Г.И., Шиманский-Гайер Л. Нелинейные эффекты в хаотических и стохастических системах. Ижевск, 2003. 530 с.
7. Mosekilde E., Maistrenko Y., Postnov D. Choatic synchronization. Application to living systems. World Scientific, 2002. 440 p.
8. Meda P., Atwater I., Goncalves A., Bangham A., Orci L., Rojas E. The topography of electrical synchrony among -cells in the mouse islet of Langerhans // Q. J. Exp. Psychol. 1984. Vol.69. P.719-735.
9. Kiss I.Z., Zhai Y., Hudson J.L. Emerging coherence in a population of chemical oscillators // Science. 2002. Vol.296. P.1676.
10. Dano S., Hynne F., Monte S. de, Ovidio F. d’, Sorensen P.G., Westerhoff H. Synchronization of glycolytic oscillations in a yeast cell population // Faraday Discuss. 2002. Vol. 120. P.261-275.
11. Postnov D.E., Sosnovtseva O.V., Mosekilde E. Oscillator clustering in a resource consumption chain // CHAOS. 2005. Vol.15. P.1-12.
12. Osipov G.V., Ivanchenko M. V., Kurths J., Hu B. Synchronized chaotic intermittent and spiking behavior in coupled map chains // Phys. Rev. E. 2005. Vol.71. P.056209.
13. West J.B. Best and Taylors Physiological Basis of Medical Practice. 12th ed. Baltimore: Williams and Wilkins, 1991.
14. Netter F.H. The CIBA collection of medical illustration // Summit: CIBA; Kidneys, Ureters and Urinary Bladder, 1973. Vol.6.
15. Despopoulos A., Silbernagel S. Color atlas of physiology. Stuttgard: Georg Thieme Verlag, 1991.
16. Holstein-Rathlou N.-H., Leyssac P.P. TGF-mediated oscillations in the proximal intratubular pressure: Differences between spontaneously hypertensive rats and Wistar-Kyoto rats // Acta Physiol. Scand. 1986. Vol.126. P.333-339.
17. Leyssac P.P., Holstein-Rathlou N.-H. Effects of various transport inhibitors on oscillating TGF pressure response in the rat // Pfflugers Arch. 1986. Vol.407. P.285.
18. Holstein-Rathlou N.-H., Wagner A.W., Marsh D.J. Tubu-loglomerular feedback dynamics and renal blood flow autoregulation in rats // Amer. J. Physiol. 1991. Vol.260. P.F53.
19. Chou K.H., Yu-Ming Ch., Mardarelis V.Z., Marsh D.J., Holstein-Rathlou N. -H. Detection of interaction between myogenic and TGF mechanisms using nonlinear analysis // Amer. J. Physiol. 1994. Vol.267. P.F160.
20. Feldberg R., Colding-Jorgensen M., Holstein-Rathlou N.-H. Analysis of interaction between TGF and the myogenic response in renal blood flow autoregulation // Amer. J. Physiol.
1995. Vol.269. P.F581.
21. Holstein-Rathlou N.-H., Marsh D.J. Oscillations of tubular pressure, flow and distal chloride concentration in rats // Amer. J. Physiol. 1989. Vol.256. P.1007-1014.
22. Yip K.-P., Holstein-Rathlou N.-H., Marsh D.J. Chaos in blood flow control in genetic and renovascular hypertensive rats // Amer. J. Physiol. 1991. Vol.261. P.400-408.
23. Moore L.C., Rich A., Casellas D. Ascending myogenic autoregulation: interactions between tubuloglomerular feedback and myogenic mechanisms // Bull. Math. Biol. 1994. Vol.56. P.391-410.
24. Oien A.H., Aukland K. A multinephron model of renal blood flow autoregulation by tubuloglomerular feedback and myogenic response // Acta. Physiol. Scand. 1991. Vol. 143. P.71-92.
25. Holstein-Rathlou N.-H. Synchronization of proximal intra-tubular pressure oscillations: evidence for interaction between nephrons // Pfflugers. Arch. 1988. Vol.408. P.438-443.
26. BarfredM., Mosekilde E., Holstein-Rathlou N.-H. Bifurcation analysis of nephron pressure and flow regulation // Chaos.
1996. Vol.6. P.280-287.
27. Marsh D.J., Sosnovtseva O.V., Chon K.H., Holstein-Rathlou N.-H. Nonlinear interactions in renal blood flow regulation // Amer. J. Physiol. 2005. Vol.288. P.R1143-R1159.
28. Jensen K.S., Mosekilde E., Holstein-Rathlou N.-H. Selfsustained oscillations and chaotic behavior in kidney pressure regulation // Mondes Develop. 1986. Vol.55. P.91-109.
29. Holstein-Rathlou N.-H., Marsh D.J. A dynamic modelof the tubuloglomerular feedback echanism // Amer. J. Physiol. 1990. Vol.258. P.F1448-F1459.
30. Chen Yu-Ming, Kay-Pong Yip, Marsh D.J., Holstein-Rathlou N.-H. Magnitude of TGF-initiated nephron-nephron interactions is increased in SHR // Amer. J. Physiol. 1995. Vo1.269. P.F198-F204.
31. Postnov D.E., Sosnovtseva O.V., Mosekilde E., Holstein-Rathlou N.-H. Cooperative phase dynamics in coupled nephrons // Intern. J. Mod. Phys. B. 2001. Vol.15. P.3079-3098.
32. Постнов Д., Шишкин A., Щербаков П. Нелинейные эффекты в ансамблях осцилляторов со связью через распределение ресурса. 4.I. Динамические режимы авторегуляции кровотока в васкулярном дереве нефронов // Изв. вузов. ПНД. 2007. Т.15, №5. С.3-22.
33. Постнов Д., Шишкин A., Щербаков П. Нелинейные эффекты в ансамблях осцилляторов со связью через распределение ресурса. 4.II. Колебательные режимы одномерного массива связанных через общий источник питания осцилляторов // Изв. вузов. ПНД. 2007. Т.15, №5. С.23-35.
34. Postnov D.E., Shishkin A. V., Sosnovtseva O. V., Mosekilde E. Two-mode chaos and synchronization properties // PRE. 2005. Vol.72. P.056208.
35. Van-der-Pol B. Forced oscillations in a circuit with nonlinear resistance // Phil. Mag. 1927. Vol.3. P. 64-80.
36. Postnov D.E., Sosnovtseva O.V., Scherbakov P.A., Mose-kilde E. Multimode dynamics in a network with resource mediated coupling // CHAOS. 2008. Vol. 18. P. 015114.
37. Постнов Д.Э., Жирин Р.A. Свидетельство об официальной регистрации программы для ЭВМ №2007614145. 2007.
38. Sosnovtseva O.V., Postnov D.E., Nekrasov A.M., Mosekilde E., Holstein-Rathlou N.-H. Phase multistability of selfmodulated oscillations // Phys. Rev. E. 2002. Vol.66. P.0362.
39. Постнов Д.Э., Некрасов A.M. Механизмы фазовой мультистабильности при синхронизации 3D осцилляторов // Изв. вузов. ПНД. 2005. Т.13, №1-2. С.47-62.
40. Некрасов A.M. Фазовая мультистабильность в диффузионно связанных нелинейных осцилляторах: Дис. ... канд. физ.-мат. наук. Саратов, 2007. 160 с.
