CC BY
24
Прикладные задачи
^^^^^^^^^^»нелинейной теории колебаний и вслн
Изв. вузов «ПНД», т. 15, № 5, 2007 УДК 530.182:577.3
НЕЛИНЕЙНЫЕ ЭФФЕКТЫ В АНСАМБЛЯХ ОСЦИЛЛЯТОРОВ СО СВЯЗЬЮ ЧЕРЕЗ РАСПРЕДЕЛЕНИЕ РЕСУРСА
Часть 1: Динамические режимы авторегуляции кровотока в васкулярном дереве нефронов
Д.Э. Постное, А.В. Шишкин, П.А. Щербаков
Исследованы характерные колебательные режимы и нелинейные эффекты, возникающие в условиях особого типа связи, который широко распространен в природе. А именно, во многих случаях взаимодействие в ансамбле осцилляторов осуществляется посредством потребления и распределения некоего энергонесущего ресурса. Динамика таких систем имеет ряд особенностей. В первой части работы показано, как детализация модели авторегуляции почечного кровотока приводит к системе интересующего нас класса и каковы ее типичные динамические режимы.
Введение
Классический подход к изучению явления синхронизации [1-3] широко использует парадигму, согласно которой взаимодействующие автоколебательные системы имеют каждая свой собственный источник энергии, тогда как наличие относительно слабой связи между ними ответственно за подстройку частот и фаз (а при сильной связи - и амплитуд) колебаний.
Разумеется, конкретный тип связи во многом определяет поведение ансамбля осцилляторов. Хорошо изученные типы связи включают, например, локальную связь [4], когда осцилляторы в решетке (или иной структуре) взаимодействуют лишь с ближайшими соседями, а также глобальную связь [5,6], когда каждый осциллятор ансамбля связан со всеми остальными его элементами. Для живых систем примером первого типа связи может служить взаимодействие в кластере бета-клеток поджелудочной железы посредством щелевых контактов [7]. Примерами систем с глобальной связью могут служить системы электрохимических осцилляторов [8] или же колонии дрожжевых клеток [9]. В последние годы популярным обьектом исследования становятся так называемые small-world networks [10,11] как системы, сочетающие в себе свойства локальных и дальних взаимодействий.
Для всех перечисленных выше типов связи математическая формулировка задачи предполагает, что нелинейные свойства каждого элемента ансамбля осцилляторов (собственная частота колебаний, чувствительность к возмущениям) управляются
его собственными параметрами, тогда как взаимодействие между осцилляторами задается отдельным набором параметров, характеризующих структуру и силу связи. Таким образом, всегда возможно разделить индивидуальный колебательный режим осциллятора и его поведение в ансамбле как результат действия связи.
Однако имеется обширный класс задач физики, химии или биологии, которые не могут успешно решаться в рамках этой парадигмы, поскольку связь между осцилляторами возникает в процессе распределения энергии (в виде несущего ее ресурса), без которой автоколебательная динамика попросту невозможна. Примером может служить связь (как правило, паразитная) нескольких электронных схем, питающихся от общего источника, или же конкуренция популяций бактерий в питательном растворе [12].
Поскольку потребление энергии (ресурса) одним из элементов ансамбля так или иначе влияет на остальные, такой тип связи формально может быть отнесен к глобальному, или локальному в зависимости от конкретной реализации. Однако его специфика заключается в том, что такие параметры, как частота и амплитуда колебаний, как правило, напрямую зависят от количества потребляемой энергии и по этой причине также зависят от характеристик связи! Если, например, геометрия системы близка к линейной (цепочка осцилляторов), то количество ресурса (энергии, питательных веществ) будет убывать по мере удаления от его источника. В результате даже изначально идентичные осцилляторы могут оказаться в различных режимах функционирования, с различными частотами и амплитудами.
В первой части работы изучаются особенности подобной связи через распределение ресурса на примере авторегуляции почечного кровотока в васкулярном дереве нефронов.
Почка играет существенную роль в регуляции давления крови и ее состава. Нефрон как мельчайшая функциональная единица почки представляет собой фильтрующее устройство с обратной связью, которая контролирует входящий в него поток крови с целью поддержания оптимального режима фильтрации. С момента экспериментального обнаружения колебательных режимов в динамике отдельных нефронов [13] предпринимаются попытки математического моделирования соответствующих процессов на уровне описания изменений давления в потоке крови и фильтрата [14-19]. Эти модели более или менее успешно воспроизводят сам процесс авторегуляции, лежащий в основе возникновения колебаний давления в отдельно взятом нефроне. Однако в живой почке каждый нефрон составляет часть сложно организованного ансамбля. А именно, группы по 10-20 нефронов образуют структуру васкулярного дерева, так как имеют общий кровеносный сосуд, доставляющий кровь в этот ансамбль и разветвляющийся на пары, тройки и одиночные нефроны подобно дереву [20]. В свою очередь, васкулярные деревья включены в аналогичную структуру более высокого уровня и так далее, вплоть до единственной артерии, входящей в почку.
Очевидно, все экспериментальные результаты следует трактовать как измерения активности нефрона, включенного, по меньшей мере, в ансамбль из 10-20 подобных единиц. При этом некоторые предположения и упрощения, использованные при моделировании одиночного нефрона, становятся неадекватными. В частности, необходимо отказаться от предположения о постоянстве артериального давления на входе в нефрон. Следует также заметить, что современные экспериментальные мето-
дики не позволяют получать информацию о режиме функционирования глубинных нефронов, расположенных у основания васкулярного дерева. Эти глубинные нефро-ны отличаются от корковых прежде всего параметрами петли Генле (см. ниже), что существенно влияет на характеристики (задержка в цепи обратной связи и ее глубина) процесса авторегуляции.
В данной работе модифицирована наиболее известная модель нефрона [16] путем включения в нее описания процессов реабсорбции в петле Генле, что улучшает соответствие с результатами экспериментов и позволяет в перспективе моделировать различия в динамике корковых и глубинных нефронов. Использование артериального давления на входе в нефрон в качестве управляющего параметра позволяет получить дополнительную информацию, важную для интерпретации результатов моделирования васкулярного дерева простейшей конфигурации (одномерная цепочка нефронов). По результатам исследования динамики васкулярного дерева обнаружены пространственно неоднородные состояния (кластерная осцилляция), возникновение которых обусловлено структурой связи нефронов. В заключение приводятся результаты численного эксперимента в условиях, когда входное артериальное давление принималось флуктуирующим относительно среднего уровня. Такой шаг в сторону приближения к реальным условиям представляется необходимым с точки зрения сопоставления с данными экспериментов.
1. Модифицированная модель авторегуляции кровотока в одиночном нефроне
Структура одиночного нефрона с некоторыми упрощениями приведена на рис. 1. Основные его части, участвующие в фильтрации и авторегуляции давления крови, это: гладкие мышечные клетки, расположенные на стенках афферентной арте-риолы и меняющие ее радиус, капсула Боумана, содержащая glomerulus (клубочек), петля Генле и клетки macula densa (плотное пятно).
Входящий поток крови проходит через афферентную артериолу к капиллярной системе клубочка, где происходит первичная фильтрация. Далее фильтрат проходит через проксимальную трубку и попадает в петлю Генле. В петле Генле происходит окончательная фильтрация, играющая важную роль в авторегуляции внутринефрон-ного давления крови. Петлю условно можно разделить на две части: тонкий ниспадающий и толстый восходящий участки. В тонком ниспадающем участке происходит отделение воды, в то время как в толстом восходящем преимущественно ре-абсорбируются ионы Na+ и Cl_, тогда как обьем жидкости меняется мало. Прежде чем покинуть систему трубок, жидкость проходит через специализированные клетки macula densa, которые реагируют на изменения в концентрации NaCl и вырабатывают сигнал обратной связи электрической природы, который управляет диаметром афферентной артериолы путем сокращения гладких мышечных клеток, расположенных на стенках артериолы. Этот механизм получил название тубуло-гломерулярной обратной связи, которая обеспечивает поддержание основных функций нефрона в условиях внешних изменений, например, давления крови. При этом нелинейный характер процессов, протекающих в системе трубок, и наличие существенной задержки в цепи обратной связи (порядка 10-20 с для нефрона крысы) приводят к
возникновению осцилляций давления, потоков и концентраций ионов с периодом 20-40 с. Было обнаружено, что такие колебания могут носить почти регулярный, либо (в других случаях) сильно нерегулярный характер [21-23].
В настоящей работе использована в качестве прототипа наиболее известная модель одиночного нефрона [16]. Она содержит шесть дифференциальных уравнений, полный вид которых приведен в Приложении. Все уравнения можно разделить на три группы: уравнения для системы трубок (проксимальной и петли Генле), уравнения для клубочка и уравнения, описывающие авторегуляцию диаметра афферентной артериолы. Три линейных уравнения учитывают необходимый сдвиг фаз в цепи обратной связи. Однако эта модель дает неудовлетворительную оценку параметров обратной связи, при которых нефрон переходит в хаотический режим. Кроме того, для целей нашего исследования представляется важным более детальный учет процессов фильтрации и реабсорбции. В рамках предложенной нами модификации были заменены уравнения, описывающие процессы в петле Генле.
На выходе из капсулы Боумана давление Pt в проксимальной трубке определяется балансом притока от клубочка Ffiit, оттока в петлю Генле Fhi и потоком реабсорбции Freab сквозь стенки.
dP 1
= ^~[FfHt(Pa, Pt, Г) " Freab - Fhl], (1)
dt Btub
где Ffiit определяется скоростью фильтрации в клубочке и зависит как от артериального давления Pa, так и от радиуса артериолы r; Freab составляет, примерно, постоянную часть от общего количества фильтрата (порядка 70%); Fhi соответствует потоку жидкости, втекающему в первый сегмент петли Генле; Btub - параметр эластичности проксимальной трубки.
Попытка усовершенствования модели [16] в части описания процессов, происходящих в петле Генле, была предпринята в работах [19, 24] с помощью связанных нелинейных уравнений в частных производных. Это позволило дать более детальное описание процессов, однако, создало ряд новых проблем таких, как адекватность применяемых методов численного интегрирования. При построении модифицированной модели нами была поставлена цель выдержать баланс между детальностью модели, с одной стороны, и простотой ее численного исследования и интерпретации результатов, с другой стороны.
В рамках предлагаемой модели петля Генле условно делится на пять сегментов (см. рис. 1). В первых двух - имеет место интенсивная реабсорбция воды, тогда как в последних трех - меняется баланс ионов NaCl. Такое деление соответствует современнным представлениям о физиологии нефрона, а с точки зрения использования модели позволяет, например, отдельно менять степень реабсорбции воды в первых сегментах, что на качественном уровне воспроизводит действие некоторых медикаментов, регулирующих артериальное давление.
Для каждого из сегментов записываются дифференциальные уравнения для давления P, для количества Q и A вещества NaCl, а также алгебраическое уравнение для концентрации C ионов NaCl. Так, для первого сегмента уравнения для давления записываются в виде
Pi = [Fhi - Fh2 - Fabsi], (2)
Bhi
Fh1=^, (3
Fh2 = P—2, (4)
Rh2
Fabsi = SiFhi, (5)
где Fhi и Fh2 описывают втекающий и вытекающий для данного сегмента потоки жидкости, соответственно; Fabsi отражает процесс реабсорбции воды; Bhi - константа, описывающая эластичность стенок сегмента; Rhi и Rh2 - гидродинамическое сопротивление первого и второго сегмента, соответственно; Si - коэфициент реабсорбции воды в первом сегменте; здесь и далее Pi, i = 1...5 - давление в конце i-го сегмента.
Уравнения для количества Q, А и концентрации С вещества №С1 в пределах первого сегмента записываются в виде
= Со^ы - ^ , (6) ТН1
А1 = ^ - СЕ^, (7)
С1 = ^ • ®
где уравнение для Ql в линейном приближении моделирует конечную скорость распространения №С1 по сегменту, а уравнение для А1 определяет количество вещества в конце первого сегмента с учетом реабсорбции. Параметр характеризует время прохождения сегмента фильтратом. Со есть концентрация ионов №С1 в начале сегмента, тогда как С1 - в его конце. Аналогично уравнениям (2), (6)-(8) записываются уравнения для второго сегмента петли Генле
Р2 — "5— [ЕН2 — ЕЬ,3 — ЕаЬв2], ВН2
Q2
0?2 — СХ¥Н2 — ^, А2 — ^ — С2
ТН2
С2 — ^, 2 ВН2Р2'
(9)
где
Р1 — Р2
ЯН2 '
Р2 — Р3 Янз '
Е«Ь«2 —
Рн2 —
ЕЪ —
(10)
(11) (12)
Здесь и описывают входящий в сегмент поток и выходящий, соответственно; РаЬз2 описывает процесс поглощения воды.
В отличие от тонкого ниспадающего участка петли Генле в толстом восходящем сегменте практически отсутствует реабсорбция воды, в то же время происходит поглощение №С1. Соответственно уравнения для давления не содержат член, отвечающий за реабсорбцию воды, но в уравнениях для количества №С1 добавляется слагаемое, описывающее поглощение №С1.
При этом уравнения для третьего сегмента принимают вид
Рз = ^[Fh3 - FM],
Bh3
■ Q3
Q3 = C2Fh3--,
Q %h3 (13)
A3 = ^^ - C3Fh4 - f (C3) -,
Th3 То
C = A3
C3 = BhP'
где то есть общее время прохождения фильтратом петли Генле, а нелинейная функция f (C) = M0C/(Cm + C) описывает интенсивность процесса реабсорбции NaCl. В уравнениях (13) входящий и выходящий потоки определяются как Fh3 = = (Р2 — P3)/Rh3 и Fh4 = (Р3 - P4)/Rh4, соответственно. Аналогично с (13) записываются уравнения и для следующего, четвертого сегмента.
Для последнего, пятого сегмента уравнения претерпевают изменение, связанное с возможностью обратного тока ионов NaCl в фильтрат на этом участке, и принимают вид
Р5 = 7^[Fh5 - Fd],
Bh5
■ q5
Q5 = C4Fh5--,
■ Q5 %h5 T3 (14)
AL5 = — - C5Fd - [f (C5) - L(Ci - C5)] -,
Th5 То
C5 = A5 _
Bh5P5
Здесь Fh5 = (P4 - P5)/Rh5 и Fd = (P5 - Pd)/Rd; Fd описывает поток фильтрата на выходе из нефрона; дистальное давление Pd принято постоянным и оценено из данных эксперимента.
Специализированные клетки macula densa измеряют концентрацию C5 ионов NaCl и вырабатывают сигнал обратной связи Ф (электрической природы), который влияет на сокращение мышц афферентной артериолы,
т ^max - ^min .-,,4 Ф = ^max - ^-г (п ,п-7VT • (15)
1 + exp[a(C5/Ceg - 1)]
В этом соотношении, вид нелинейности которого взят из [16], величины ^max и ^m¡n определяют максимальную и минимальную активацию мышечных клеток, а определяет наклон кривой обратной связи, Ceq - равновесную концентрацию ионов NaCl.
Оставшаяся часть модели одиночного нефрона [16] нами не изменялась. Она содержит уравнения для регуляции радиуса атериолы, которые записываются в виде
Г = Vr, Vr = Ю {Pav (Pt, Г) - Peq (r, C5, a, T) - wdVr } , (16)
где r - радиус активной части артериоллы, vr - скорость его изменения, ю определяет эластичность стенки артериолы, а определяет силу обратной связи. Величины Ffiit, Pav, Peq находятся из алгебраических соотношений, которые должны решаться
вместе с интегрированием всей системы. Детальное описание нелинейных функций Раю, Рещ и рекомендации по моделированию процессов в этой части нефрона можно найти в работах [14-19]. Значения, размерность и смысл параметров математической модели нефрона в рамках предлагаемой модификации приведены в таблице 1.
Таблица 1
Параметры модели одиночного нефрона
Со 150 м-моль/л Начальная концентрация ионов №С1
Ям 1.0 кПа-с/(н-л) Гидродинамическое сопротивление в первом сегменте
Ян2 2.0 кПа-с/(н-л) Гидродинамическое сопротивление во втором сегменте
Пьз 0.8 кПа-с/(н-л) Гидродинамическое сопротивление в третьем сегменте
Як4 0.8 кПа-с/(н-л) Гидродинамическое сопротивлениее в четвертом сегменте
Як5 0.7 кПа-с/(н-л) Гидродинамическое сопротивление в пятом сегменте
Яа 1.0 кПа-с/(н-л) Гидродинамическое сопротивление в дистальной трубке
Вгиь 1.9 н-л/кПа Эластические свойства проксимальной трубки
Мо 90 м-моль/(с-л) Максимальная реабсорбция №С1
Ст 35 м-моль/л Константа половинной скорости для реабсорбции в 3, 4, 5 сегментах
БН1 0.8 н-л/кПа Эластичность стенок первого сегмента
Вь2 0.7 н-л/кПа Эластичность стенок для второго сегмента
Ем 0.2 н-л/кПа Эластичность стенок для третьего сегмента
Вы 0.2 н-л/кПа Эластичность стенока для четвертого сегмента
Еьь 0.2 н-л/кПа Эластичность стенок для пятого сегмента
ТН1 0.6 с Временная задержка в первом сегменте
ТН2 0.6 с Временная задержка во втором сегменте
Тья 2.8 с Временная задержка в третьем сегменте
Тм 2.8 с Временная задержка в четвертом сегменте
Тьь 2.8 с Временная задержка в пятом сегменте
Сец 42.0 м-моль/л Равновесная концентрация №С1
ь 0.38 с"1 Диффузионная константа обратного тока №С1
С 45.0 м-моль/л Концентрация N0 вне петли Генле
81 0.35 Коэфициент поглощения воды в первом сегменте
82 0.27 Коэфициент поглощения воды во втором сегменте
То 12 с Время задержки в петле
2. Колебательные режимы модели одиночного нефрона
Модель одиночного нефрона имеет большое количество параметров. Для сопоставления с данными эксперимента имеет смысл выбрать в качестве основных один-два, интервал значений которых известен и/или оценка которых представляет интерес. Такими параметрами были выбраны величина артериального давления Ра и параметр крутизны кривой обратной связи а.
Необходимым шагом после разработки модифицированной модели нефрона было ее тестирование: выявление характерных режимов функционирования и значений основных переменных и сопоставления с экспериментальными оценками. Следует отметить, что при разработке модели в работе [16] выбор вида и количества (в данном случае трех) дифференциальных уравнений, описывающих динамику петли, осуществлялся чисто эмпирически, путем подгонки получаемых результатов под экспериментальные. В нашей модели, основанной на более поздних экспериментальных оценках большего количества параметров, количественные характеристики колебаний должны бы находиться в разумных пределах без дополнительных подстроек.
На рис. 2 приведены временные реализации и проекции фазового портрета динамических режимов модели (2)-(16) при Ра = 13.3 кПа. Во-первых, отметим сам факт наличия автоколебаний при выбранных значениях параметров петли Ген-ле. Это означает, что данная совокупность значений, как минимум, правдоподобна. Более тонкой деталью является наличие мелкомасштабных осцилляций на вершине каждого максимума, ожидаемых с точки зрения механизма регуляции и экспериментальных данных, но слабо выраженных в динамике прежней модели [16]. Экспе-
Рис. 2. Временная реализация и проекция фазового портрета колебаний в модифицированной модели нефрона при Ра = 13.3 кПа: а - периодические колебания при а = 10.5; б - хаотические колебания при а = 12.5
риментальные измерения дают интервал 1.5-2.0 кПа. Таким образом, диапазон изменения для модели находится в пределах экспериментальной оценки.
Более сложная ситуация с временными характеристиками. Использованная нами аппроксимация распространения фильтрата вдоль петли Генле подразумевает, что в пределах каждого из пяти сегментов распределение давления и концентраций близко к линейному. Адекватность такого представления трудно оценить заранее. В нашем случае, совместное действие всех параметров временных масштабов дает период колебаний порядка 40 с, что укладывается в рамки экспериментальных оценок, хотя и выглядит несколько завышенным.
Особый интерес представляет сопоставление значений параметра а и типа колебательного режима. Как было показано в экспериментальных работах [13], колебания становятся существенно нерегулярными для лабораторных крыс с повышенным кровяным давлением. При этом, оценки дают а ~ 12... 15 против а < 12 для крыс с нормальным кровяным давлением. Интересно, что в модели [16] эти оценки воспроизвести не удается, там переход к хаотическому режиму наблюдается лишь при а ~ 26. Как следует из рис. 2, в модифицированной модели при Ра = 13.3 кПа увеличение а до значения 12.5 приводит к развитию хаотического режима с характерной структурой фазового портрета, в котором мелкие петли отражают относительно быстрые колебания радиуса артериолы.
В целом, можно сказать, что предложенная модификация модели по ряду характеристик лучше соответствует экспериментальным данным и потому может быть использована для целей моделирования динамики более крупной структуры нефронов.
При использовании модели (2)-(16) в качестве элемента группы взаимодействующих нефронов принципиальную роль играет ее реакция на изменение артериального давления Ра на входе в афферентную артериолу. На рис. 3 приведены однопараметрические бифуркационные диаграммы, полученные следующим образом. Для текущего значения Ра при интегрировании уравнений модели фиксируются значения в момент их локального минимума или максимума. Таким образом, единственная линия значений для некоторого значения Ра означает отсутствие колебаний (минимум и максимум совпадают). Наличие двух или более четко
Рис. 3. Однопараметрическая бифуркационная диаграмма при а = 12.5 (а) и при а = 22.0 (б). Серым фоном дан размах колебаний по переменной Рг
выраженных линий означает наличие периодических колебаний. Неупорядоченная структура значений Р1 соответствует хаотическому режиму.
Как можно видеть из рис. 3, имеется ограниченный диапазон значений артериального давления 12.90кПа<Ра< 14.86кПа, в пределах которого модель демонстрирует колебания. При этом области как низких, так и высоких значений артериального давления из этого диапазона соответствуют области хаотической динамики, тогда как при промежуточных значениях Ра наблюдаются периодические колебания.
Одномерная аппроксимация структуры васкулярного дерева нефронов
Кровеносные сосуды почки в целом образуют сложную, ветвящуюся, многоуровневую структуру (см., например, недавние работы [25], где эта структура визуализирована без нарушения целостности обьекта). Расположение нефронов также привязано к этой структуре. Как уже упоминалось во Введении, типичной структурой, которую образуют нефроны, является «дерево» из 12-20 единиц. По мере продвижения от «основания» дерева к его «верхушке» падение давления на входе нефронов может составлять до 5.32 кПа. Очевидно, изменения кровотока на входе в отдельный нефрон в большей или меньшей степени скажутся на режиме работы остальных нефронов дерева.
Несмотря на то, что реальная пространственная структура нефронного дерева нерегулярна (нефроны могут располагаться как поодиночке, так и парами или тройками), на начальном этапе моделирования разумной аппроксимацией представляется простейшая линейная геометрия дерева, где нефроны расположены поодиночке, последовательно и через равные расстояния (рис. 4).
В рамках такой упрощенной модели дерева предполагается, что каждый нефрон связан с соответствующей точкой ветвления интерлобулярной артерии через арте-риолу с гемодинамическим сопротивлением Щ, ] = 1...10. Артериальное давление Ра, входящее в модель одиночного нефрона, теперь становится давлением крови в конкретной точке ветвления Ра. Участки между точками ветвления интерлобулярной артерии характеризуются их гемодинамическими сопротивлениями Щ, ] = 1...10. Тогда, изменение давления крови в каждой из точек ветвления васкулярного дерева описывается уравнением
Рис. 4. Схематическое изображение упрощенной структуры васкулярного дерева. Десять нефронов расположены последовательно и через равные расстояния. Р0 - артериальное давление в основании дерева.
йРа Б0;-3
3 <и
аа Р3-1 Р3
Щ
аа Р3 Р3+1
Щ+1
Ра _ р9
Рз Р3
(17)
где Ба = 1-1,] = 1...10, характеризует эластичность стенок артерии, а Р3 есть аналог давления в клубочке Рд для модели одиночного нефрона. Таким образом, васкулярная модель дерева нефронов состоит из набора уравнений, описывающих одиночные нефроны, а также набора уравнений вида (17), которые определяют ди-
намику давления в каждой из точек ветвления. Общее число дифференциальных уравнений модели составляет 190.
В качестве управляющих параметров были выбраны Ро и набор гемодина-мических сопротивлений Щ, ] = 1...10. Так как наша основная цель состоит в выявлении основных свойств такого типа связи в ансамбле, все 10 нефронов были взяты идентичными. Приведенные ниже результаты численного моделирования получены для случая равных сопротивлений сегментов интерлобулярной артерии, Щ = 0.035 кПа с/(н л), ] = 1...10.
4. Индуцированная связью неоднородность режимов
На рис. 5 приведена амплитуда колебаний по переменной в каждом из десяти нефронов дерева при различных значениях артериального давления Ро. Как можно видеть, при Ро = 15.3 кПа (рис. 5, а) амплитуды колебаний отличны от нуля только у основания дерева (точки ветвления 1, 2, 3, 4). Нефроны, расположенные у этих точек ветвления, находятся в автоколебательном режиме, в то время как расположенные дальше по дереву нефроны поддерживают постоянное (отличное от нуля) значение давления в проксимальной трубке Ръ. При большем Ро = 18.3 кПа (рис. 5, б) ситуация иная: амплитуда колебаний плавно нарастает от первого к третьему нефрону; далее резкий скачок: нефроны с номерами 4, 5 и 6 имеют значительно больший размах колебаний. Далее амплитуда резко снижается, в 7-м нефроне она примерно такая же, как и в 3-м, и продолжает плавно убывать до верхней оконечности дерева. Таким образом, в одномерном массиве нефронов образуется кластер элементов, находящихся в автоколебательном режиме, тогда как остальные находятся в пассивном состоянии.
Дополнительную информацию о таком режиме функционирования исследуемой системы дает рис. 6, где приведены проекции фазового портрета на плоскости переменных и г каждого из 10 нефронов при Ро = 18.3 кПа. Четко видно, что фазовая проекция для 4, 5 и 6 нефронов существенно отличается не только большим размером, но и структурой петель траектории. На трех вставках внизу рисунка
Рис. 5. Амплитуда колебаний по переменной Рг для каждого из нефронов в зависимости от их расположения в дереве ]. Давление в основании дерева Р0 составляло: а - 15.3 кПа, б - 18.3 кПа
0.85
0.30
1 2 3 4 5
д (Цр 6 -1- 7 0 8 9 10
1.215
1.200
1.375
1.360
3.050 3.075
Рис. 6. Фазовые проекции колебаний во всех нефронах дерева. На вставках внизу показана структура аттракторов для 1,7 и 10 нефронов. Значение артериального давления Р0 = 18.3 кПа
приведены в укрупненном масштабе фазовые проекции для тех нефронов, где не наблюдаются колебания большой амплитуды. Вид фазовой проекции в них обьяс-няется тем, что вызванные другими нефронами колебания Ра на входе в нефрон играют роль сигнала воздействия сложной формы, вызывая непрерывную подстройку Рг и г.
Качественное объяснение наблюдаемого эффекта может быть предложено на основе информации о динамике модели одиночного нефрона. Как мы уже видели на рис. 3, область его автоколебательной динамики ограничена в интервале между значениями артериального давления 12.9 кПа и 14.86 кПа. Разумно предположить, что за образование кластеров генерации отвечает величина среднего значения артериального давления в точках ветвления для каждого из нефронов. Другими словами, если среднее значение артериального давления Р а пересекает пределы зоны, где возможна автоколебательная динамика, происходит смена режима, в котором находится нефрон. Таким образом, появление осцилляторного кластера объясняется в терминах падения давления от одной точки ветвления к другой, вызывающего изменения динамического режима, в котором находится нефрон.
Одно из существенных упрощений модели одиночного дерева нефронов заключается в том, что давление Ро в его основании принято постоянным. Очевидно, в живой системе это не так. Помимо неизбежного влияния на его величину соседних деревьев постоянно меняется артериальное давление в организме в целом, отзываясь на любое изменение его активности. Для режима функционирования дерева нефронов существенным может оказаться не только модуляция давления в точках ветвления элементами осцилляторного кластера, но и сходное по действию изменение Ро. Мы провели отдельный численный эксперимент по выяснению влияния флуктуаций артериального давления на режим кластерной осцилляции. При этом,
поскольку флуктуации артериального давления обусловлены не только случайными событиями, но и наложением многих ритмических процессов (сердцебиение, дыхание и др.), они не отвечают популярной у исследователей идеализированной модели 8-коррелированного (белого) шума.
Не претендуя на то, чтобы детально воспроизвести характер природных флук-туаций, мы пользовались следующим методом для изготовления суррогатных данных. По экспериментальным данным оценивалась верхняя граничная частота /тах сигнала артериального давления в сосудах почки крысы. Затем использовался компьютерный генератор гауссова белого шума, последовательность значений которого подвергалась фильтрации, при которой отбрасывались все частотные составляющие, превосходящие /тах = 1Гц. Полученный шумовой сигнал с нулевым средним и интенсивностью О, кПа2 смещался на величину 18.3 кПа и использовался в качестве временной реализации при интегрировании уравнений системы методом Рунге-Кутты четвертого порядка с фиксированным временным шагом.
Рис. 7 показывает, как меняется положение и размер кластера в присутствии шума по отношению к рис. 5, б. Несмотря на то, что среднее значение артериального давления в обоих случаях одинаково, видно, что в присутствии шума кластер расширяется, одновременно меняя профиль амплитуд. Интересно, что нефроны у вершины дерева находятся вне кластера, их амплитуда колебаний отражает действие флуктуаций, но практически совпадает для 8, 9 и 10 нефронов.
По всей видимости, действие шума на исследуемую систему, в первую очередь, определяется тем, что флуктуирующее Ра на какую-то долю времени оказывается в зоне автоколебаний даже для нефронов, далеких от автоколебательного режима. Можно предположить, что наиболее важную роль при этом играют низкочастотные составляющие шумового сигнала, так как нефроны просто не успевают среагировать на скачки артериального давления на интервалах времени много короче, чем характерный период колебаний индивидуального нефрона. Такой вывод подтверждается тем, что выбор различных оценок /тах в диапазоне 1 . . . 3 Гц не оказывает видимого влияния на результаты.
Рис. 7. Показаны максимальная амлитуда Рг для каждого из нефронов в дереве в присутствии шума. В сравнении с рис. 5, б хорошо заметно расширение кластера в сторону основания нефронного дерева
Обсуждение результатов и выводы
Нами предложена модификация модели одиночного нефрона, заключающаяся в значительно более детальном описании процессов распространения, фильтрации и реабсорбции в петле Генле. Подход, по сути, является компромиссом между, с одной стороны, реалистичной моделью в виде системы связанных нелинейных уравнений в частных производных [24], которая чрезвычайно сложна и неудобна для компьютерного счета и бифуркационного анализа, и, с другой стороны, чрезмерно упрощен-
ным феноменологическим описанием в модели [16], которое не вызывает проблем в численном моделировании, однако априори не допускает вариации некоторых существенных параметров, таких, например, как степень реабсорбции, или особенности геометрии петли Генле для глубинных нефронов. По результатам тестирования модифицированная модель показала хорошее соответствие данным экспериментов, в частности, по значениям параметра обратной связи а.
Разработанная модификация модели использована при построении структуры более высокого уровня - нефронного дерева. Заметим, что на данный момент авторам известно крайне небольшое число работ, где предприняты такого рода попытки [20,26]. Наша модель является наиболее детальной, во всяком случае, по числу дифференциальных уравнений.
В данной работе мы не ставили задачу детального исследования динамических режимов столь сложной системы, сосредоточившись на первом этапе на наиболее общих закономерностях, к числу которых, в первую очередь, относится эффект индуцированной связью неоднородности режимов: изначально идентичные элементы нефронного дерева оказываются в различных режимах функционирования в силу особенностей такого рода связи, а именно, падения артериального давления от одной точки ветвления к другой. В результате возникает осцилляторный кластер, то есть группа элементов системы, в которых устанавливается автоколебательный режим.
Несмотря на вполне однозначную трактовку результатов в целом, можно задавать целый ряд более тонких вопросов, например:
• Как более точно определить границы осцилляторного кластера? Как следует из рис. 3, зона хаотической динамики при малых Pin характеризуется относительно малым размахом колебаний, сравнимым с величиной модуляции артериального давления соседними элементами. Возможно ли диагностировать переход элемента системы в активный режим методами теории устойчивости?
• При каких условиях колебания в элементах осцилляторного кластера могут синхронизоваться, а при каких - нет? Ведь изменение степени связи соседних элементов одновременно меняет и характеристики колебаний в них, в том числе, частоту. Таким образом, привычный при решении задач синхронизации анализ плоскости параметров (расстройка частот - сила связи) попросту невозможен!
• Каков механизм активации шумом колебаний в элементах на границах кластера? Как он соотносится с интенсивно изучаемым в последнее время эффектом когерентного резонанса? Какую роль при действии шума играет модуляция режима, вызванная активностью элементов осцилляторного кластера?
Заметим, что сформулированные выше вопросы, по сути, относятся не конкретно к авторегуляции почечного кровотока, но к свойствам более общего класса систем, где автоколебательные элементы объединены в структуру распределения некоего ресурса. К решению таких задач разумно подходить с использованием по возможности более общей и простой модели, сохраняющей наиболее существенные для исследуемых эффектов свойства. Подобный подход и будет развит во второй части данной работы.
Авторы благодарят О.Сосновцеву и Э.Мозекильде (Датский технический университет) за плодотворные обсуждения и консультации при разработке модифицированной модели нефрона.
Статья написана по материалам исследований, проведенных в рамках работ по гос. контракту 02.512.11.2111.
Приложение
Модель одиночного нефрона по работе [16]
Физический (биологический) смысл переменных модели [16] следующий: Pt - давление в проксимальной трубке; r - текущий радиус артериолы;
vr - скорость изменения текущего радиуса артериолы; Xi - вспомогательная переменная 1; X2 - вспомогательная переменная 2; X3 - вспомогательная переменная 3.
В приведенной ниже системе первое уравнение описывает процесс фильтрации в капсуле Боумана. Второе и третье уравнения моделируют процесс регуляции радиуса сосуда под воздействием потенциала активации Ф, который, помимо прочего, управляется переменной X3. Три оставшихся уравнения для переменных Xi, X2, X3 обеспечивают задержку сигнала на время T при распространении по петле Генле. Уравнения имеют вид
Pt ={Ffiit - Freab - (Pt - Pd)/Rnen)},
ctub
Г = vr,
Vr = - {Pav - Peq - -dVr},
- (18) Xi = Pt/RHen - Pd/RHen - (3/T)Xi,
3
X2 = T (X1 - X2)' 3
X3 = T (X2 - X3),
где Ffnt есть скорость фильтрации в клубочке, Ctub определяет эластичность стенок проксимальной трубки. Скорость клубочковой фильтрации определяется из соотношения
C P P
Ffiit = (1 - Яа)(1 - ca)Pa—Pg, (19)
Ce Ra
где Ha - гематокрит крови в афферентной артериоле, Ca и Ce - концентрация протеинов в афферентной и эфферентной плазме, Ra - гидродинамическое сопротивление афферентной артериолы, Pa - артериальное давление и Pg - давление в клубочке. Обратная связь через macula densa описывается сигмоидным соотношением между функцией активации мышечных клеток артериолы и потоком X3 из петли Генле
Ф = W--W ~XmÍ"-• (20)
1 + exp(a(- 1))
1 FHen0
В этом соотношении ^max и определяют максимальную и минимальную активацию мышечных клеток, а определяет наклон кривой обратной связи.
Афферентная артериола разделена на две части. Первая - имеет постоянное гидродинамическое сопротивление, в то время как во второй части диаметр может меняться и, следовательно, гидродинамическое сопротивление зависит от обратной связи
Щ = Яа0[|3 + (1 - |3)Г-4], (21)
где Па0 есть начальная величина сопротивления, а г - радиус активной части арте-риолы. В уравнении для скорости изменения радиуса ьг среднее давление в активной части артериолы определяется как
Р =1
Раь = 2
Ра - (Ра - Рд)вЩ0 + Рд
(22)
Величина давления при котором артериола находится в равновесии с данным радиусом и величиной мышечной активации определяется как
4 7
Ред = 2.4е10'0(г-1'4) + 1.6(г - 1.0) + + ^3.0(0.4^) + 7.2(г + 0.9)]. (23)
Полный набор всех констант, параметров, вычисляемых величин и промежуточных функций по модели (18) вынесен в таблицу 2.
Таблица 2
Параметры модели одиночного нефрона по работе [16]
Ра 13.3 кПа Артериальное давление крови
РV 1.3 кПа Венозное давление крови
Ра 1.3 кПа Давление в дистальной трубке
Ва0 2.3 кПа-с/(н-л) Гидродинамическое сопротивление афферентной артериолы
Ке 1.9 кПа-с/(н-л) Гидродинамическое сопротивление эфферентной артериолы
Внеп 5.3 кПа-с/(н-л) Гидродинамическое сопротивление в петле Генле
СгиЬ 3.0 н-л/кПа Эластические свойства проксимального канальца
ЁНепО 0.2 н-л/с Поток в петле Генле
ЁгеаЬ 0.3 н-л/с Поглощение в проксимальном канальце
На 0.5 Артериальный гематокрит
Са 54 г/л Концентрация протеинов в плазме
а 22-10~3 кПа-л/г Линейный осмотический коэффициент
Ь 0.39-10"3 кПа-(л/г)2 Нелинейный осмотический коэффициент
0.20 Минимальный уровень активации
^шах 0.44 Максимальный уровень активации
0.38 Равновесный уровень активации
а 0.04 с"1 Коэфициент затухания
ю 20 кПа-с2 Параметр контролирующий частоту колебаний радиуса артериолы
в 0.67 Относительная часть пассивного куска артериолы
Т 16 с Временная задержка в петле Генле
а 12 Коэффициент обратной связи
Библиографический список
1. Winfree T. The geometry of biological time. 2nd ed. New York: Springer, 2001.
2. Strogatz S.H. Sync: the emerging science of spontaneous order. New York: Hyperion, 2003.
3. Pikovsky A., Rosenblum M. and Kurths /.Synchronization: A universal concept in nonlinear science. Cambridge: Cambridge University Press, 2001.
4. Afraimovich VS. and Nekorkin V.I. Chaos of traveling waves in a discrete chain of diffusively coupled maps // Int. J. Bifurcation and Chaos. Appl. Sci. Eng. 1994. Vol. 4. P. 631.
5. Kuramoto Y Chemical oscillations, waves and turbulence. Berlin: Springer, 1984.
6. Kaneko K. Globally coupled circle maps // Physica D. 1991. Vol. 54. P 5.
7. Meda P., Atwater I., Goncalves A., Bangham A., Orci L. andRojas E.The topography of electrical synchrony among |3-cells in the mouse islet of Langerhans // Q. J. Exp. Psychol. 1984. Vol. 69. P. 719.
8. Kiss Z., Zhai Y and Hudson J.L. Emerging coherence in a population of chemical oscillators // Science. 2002. Vol. 296. P. 1676.
9. Dano F. Hynne, De Monte S., d'Ovidio F., Sorensen P.G., and WesterhoffH. Synchronization of glycolytic oscillations in a yeast cell population // Faraday Discuss. 2002. Vol. 120. P. 261.
10. Watts /.Small Worlds: The dynamics of networks between order and randomnes. Princeton: Princeton University Press, 1999.
11. Mathias N. and Gopal V. Small worlds: How and why // Phys. Rev. E. 2001. Vol. 63. P. 021117.
12. Mosekilde E., Maistrenko Y and Postnov D. Chaotic synchronization // Application to living systems, World Scientific, 2002.
13. Holstein-Rathlou N.-H. and 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.
14. Jensen K.S., Mosekilde E. and Holstein-Rathlou N.-H. Selfsustained oscillations and chaotic behavior in kidney pressure regulation // Mondes Develop. 1986. Vol. 54/55. P. 91.
15. Holstein-Rathlou N.-H. and Marsh D.J. A dynamic model of the tubuloglumerular feedback mechanism // Am. J. Physiol. 1990. Vol. 258. F1448-F1459.
16. Barfred M., Mosekilde E. and Holstein-Rathlou N.-H. Bifurcation analysis of nephron pressure and flow regulation // Chaos. 1996. Vol. 6. P. 280.
17. Holstein-Rathlou N.-H. and Marsh D.J. A dynamic model of renal blood flow autoregulation // Bull. Math. Biol. 1994. Vol. 56. P. 441.
18. Young D.K. and Marsh D./.Pulse wave propagation in rat renal tubules: Implications for GFR autoregulation // Am. J. Physiol. 1981. Vol. 240. P. F446-F458.
19. Marsh D.J., Sosnovtseva O.V., Chon K.H. and Holstein-Rathlou N.-H. Nonlinear interactions in renal blood flow regulation // Am. J. Physiol. 2005. Vol. 288. P. R1143.
20. Oien A.H. and Aukland K. A multinephron model of renal blood flow autoregulation
by tubuloglomerular feedback and myogenic response // Acta Physiol. Scand. 1991. Vol. 143. P. 71.
21. Holstein-Rathlou N.-H. and Leyssac P.P. TGF-mediated oscillations in the proximal intratubular pressure: Differences between spontaneously hypertensive rats and Witstar-Kyoto rats // Acta Physical. Scand. 1986. Vol. 126. P. 333.
22. Holstein-Rathlou N.-H. and Marsh D.J.Oscillations of tubular pressure, flow and distal chloride concentration in rats // Am. J. Physiol. 1989. Vol. 256. P. F1007.
23. Yip K.-P., Holstein-Rathlou N.-H. and Marsh D.J. Chaos in blood flow control in genetic and renovascular hypertensive rats // Am. J. Physiol. 1991. Vol. 261. P. F400.
24. Holstein-Rathlou N.-H. and Marsh D.J. A dynamic model of renal blood flow autoregulation // Bulletin of Mathematical Biology. 1994. Vol. 563. P. 411.
25. Casellas D., Dupont M., Bouriquet N., Moore L.C., Artuso A. and Mimran A. Anatomic pairing of afferent arterioles and renin cell distribution in rat kidneys // Am. J. Physiol. 1994. Vol. 267. P. F931.
26. Postnov D.E., Sosnovtseva O.V. and Mosekilde E. Oscillator clustering in resource distribution chain // Chaos. 2005. Vol. 15. P. 1.
Саратовский государственный Поступила в редакцию 12.05.2006
университет После доработки 15.05.2007
NONLINEAR EFFECTS IN ENSEMBLES OF OSCILLATORS WITH RESOURCE DISTRIBUTION COUPLING
Part 1: Dynamical regimes of blood flow autoregulation in vascular nephron tree
D.E. Postnov, A.V. Shishkin, P.A. Shcherbakov
We study the typical oscillatory regimes and nonlinear affects related to the specific but widely distributed in nature type of coupling. Namely the interaction in an ensemble of oscillators occurs due to the consumption and distribution of some energy-related resource. Such systems manifest the number of specific features. In the first part of our work we show that the detailed model of kidney blood flow autoregulation delivers such resource sharing system that belongs to such class and investigate its typical regimes.
Постнов Дмитрий Энгелевич - профессор кафедры радиофизики и нелинейной динамики СГУ. Доктор физико-математических наук (2001). Область научных интересов - сложная динамика математических моделей биологических систем, индуцированные шумом эффекты в нелинейных динамических системах. Автор 67 научных статей и книг «Chaotic Synchronization. Application to Living Systems» (World Scientific, 2002) и «Synchronization: from Simple to Complex» (Springer, 2007). E-mail: [email protected]
Шишкин Александр Владиславович - окончил Саратовский государственный университет, кафедру радиофизики и нелинейной динамики. Кандидат физико-математических наук (2007, СГУ). Область научных интересов - сложная динамика математических моделей нефронного дерева, индуцированные шумом эффекты, наблюдаемые в нефронном дереве. Автор 7 научных статей. E-mail: [email protected]
Щербаков Павел Александрович - окончил Саратовский государственный университет, кафедру радиофизики и нелинейной динамики (2004). В настоящее время является аспирантом второго года обучения. Научные интересы -динамика автоколебательных систем со связью посредством распределения ресурсов, радиофизическое моделирование. E-mail: [email protected]
