eng
Выпуск #3-4/2019
О.П.Познанский
Метод ренормализационной группы: улучшение модели Сена-Бассера для белого вещества мозга
Метод ренормализационной группы: улучшение модели Сена-Бассера для белого вещества мозга
Просмотры: 2918
Разработан метод ренормализационной группы для исследования влияния большого разброса геометрических и физических микропараметров на эффективный коэффициент диффузии. Основной элемент метода – перенормировка структурных параметров от микро- до миллиметрового масштаба случайной геометрии многокомпонентного белого вещества головного мозга. Подход улучшает модель Сена-Бассера для белого вещества мозга благодаря учету влияния беспорядка, позволяет количественно исследовать чувствительность эффективного коэффициента диффузии к вариациям доминирующего набора микропараметров.
Теги: effective diffusion coefficien renormalization group method sen-basser model метод ренормализационной группы модель сена-бассера эффективный коэффициент диффузии
DOI: 10.22184/1993-8578.2019.12.3-4.166.185
Получено: 24.12.2018 г.
Мы разработали метод ренормализационной группы, чтобы исследовать влияние большого разброса геометрических и физических микропараметров на эффективный коэффициент диффузии. Основным элементом этого метода является перенормировка структурных параметров от микро- до миллиметрового масштаба случайной геометрии многокомпонентного белого вещества головного мозга. Наш подход улучшает модель Сена-Бассера для белого вещества мозга [Sen P. N., Basser P. J., Biophys. J., 89 (2005) 2927] благодаря учету влияния беспорядка и позволяет количественно исследовать чувствительность эффективного коэффициента диффузии к вариациям доминирующего набора микропараметров.
We have developed a renormalization group method in order to explore the influence of a large range of geometrical and physical microparameters on the effective diffusion coefficient. The main element of this method comprises of the renormalization of structural parameters from the micro- to the millimeter scale of the random multi-component white matter of the brain. Our approach enhances the Sen-Basser model of the brain white matter [Sen P. N., Basser P. J., Biophys. J., 89 (2005) 2927] by taking the influence of disorder into consideration and allowing quantitative investigation of the sensitivity of the effective diffusion coefficient to the variations of the dominant set of microparameters.
ВВЕДЕНИЕ
Белое вещество головного мозга можно разделить на две области: внутриклеточное пространство (ICS) и внеклеточное пространство (ECS). ICS состоит из плотно упакованного композита аксонов, глии и их клеточных образований, а ECS представляет собой желеобразную матрицу, которая разделяет клетки биологической ткани. ECS поддерживает основные функции мозга, включая межклеточный обмен, транспортировку питательных веществ и доставку лекарств. ECS представляет только ~ 20% объема здорового головного мозга, но может резко сокращаться и динамически расширяться при изменении объема клеток головного мозга при стимуляции нейронов или при введении осмотических агентов. Изменения в размере и геометрии ECS сопровождаются многими патологическими отклонениями, указывающими на ишемию, воспаление и прогрессирующую опухоль.
Принимая во внимание, что мозг представляет собой среду, состоящую из случайно расположенных доменов разных фаз, ее эффективные макроскопические свойства могут быть описаны методами, разработанными для случайных гетерогенных материалов. Используя современные методы магнитно-резонансной томографии (МРТ), можно получить уникальную информацию об организации белого вещества, нейронных путей и их визуализации. Макроскопические свойства можно оценить, используя усреднение микроструктуры. Например, объемная доля ECS (p) и эффективный коэффициент диффузии (Deff) являются двумя основными макроскопическими параметрами, обычно используемыми для описания микроструктуры мозга. Таким образом, проблему можно сформулировать следующим образом: используя методы МРТ, можно оценить макроскопический коэффициент диффузии биологической ткани. Кроме того, можно оценить чувствительность коэффициента диффузии к объемной доле включенных фаз, отношения микроскопических коэффициентов диффузии, неупорядоченной геометрии доменов и концентрации протонов в различных областях.
Понимание влияния всех этих параметров на диффузию важно для клинических приложений, поскольку их существенное отличие от средних значений является признаком патологических состояний мозга, что дает возможность обнаруживать заболевания на самых ранних стадиях.
В экспериментах с МРТ диффузия измеряется как часть ослабления сигнала из-за декогеренции в фазовых сдвигах спинов молекул воды, вызванных поступательным движением в присутствии внешнего линейного магнитного поля. Диффузия задает фундаментальное ограничение на разрешение изображений в МРТ, но, с другой стороны, может обеспечить контрастный механизм, который позволяет визуализировать спектр молекулярного смещения. При этом предполагается идеальная система градиентных катушек и компенсированные локальные магнитные неоднородности для несмещенной оценки параметров диффузии. В процессе затухания сигнала релаксация намагниченности характеризуется продольным (T1) и поперечным (T2) временами. Релаксация может быть разной для разных доменов из-за изменения локальной намагниченности. Тем не менее при диффузии за длительный период время затухания сигнала может быть описано одной экспоненциальной кривой с заданным коэффициентом диффузии:
S/S0 = exp(–b∙Deff), (1)
где b = γ2G2δ2 (Δ – δ/3) – коэффициент, определенный по формуле Стежскал-Таннера [1, 21] (γ является гиромагнитной постоянной, S0 константа, связанная с временами релаксаций, Δ является временем диффузии и G является градиентом магнитного поля, характеризуемым длительностью δ).
В предпринятых попытках моделирования микроструктуры мозга было показано, что анизотропные свойства коэффициента диффузии хорошо согласуются с моделями ограниченной диффузии, особенно в режиме длительного временного периода. Такая модель диффузии хорошо описывает анизотропию и отслеживает волокна в белом веществе головного мозга. Обычный подход, применяемый для моделирования гетерогенной природы ткани головного мозга, заключается в изучении волокон, представленных в виде цилиндров, расположенных в различных типах симметричных решеток, таких как квадрат или шестиугольник, и использовании аппроксимации эффективной среды Максвелла – Гарнетта [2] или последовательно-параллельного приближения [3]. В таких моделях трудно вводить случайность в структуре, которая очень естественна в биологических материалах, но совершенно противоположна кристаллам твердого тела. Используя предложенный в этой статье подход в диффузии, можно исследовать неупорядоченные биологические материалы и учитывать влияние сложной геометрии и топологии на эффективные транспортные (диффузионные) свойства ткани.
Типичная диффузионно-взвешенная МРТ-последовательность ЭПП (эхо-планарная последовательность) позволяет получать информацию о дефазировке спинов молекулы воды в течение Td = 10 ÷ 60 мс. Коэффициент самодиффузии чистой воды при температуре 36,6 °C равен 2,94 ∙ 10–9 м2/с. Мы считаем, что присутствие глии-клеток во внеклеточном пространстве согласуется с ограниченным коэффициентом диффузии, который ниже на 30%, чем у свободной воды при температуре тела и примерно равен 2 ∙ 10–9 м2/с. Кроме того, внутриклеточная диффузия сильно зависит от нейрофиламентов и микротрубочек и, следовательно, меньше на 75% диффузии воды и оценивается как 0,75 ∙ 10–9 м2/с. Используя уравнение диффузии Эйнштейна (<x2> = 2DTd), можно вычислить среднее смещение молекулы воды (<x2>), которое равно 5 ÷ 15 мкм. Диаметр волокон примерно равен 3 мкм. Среднее расстояние между волокнами случайно и находится в интервале 1 ÷ 18 мкм. Легко понять, что молекула воды зондирует все составляющие компоненты биологической ткани, а МРТ-сигнал дает усредненный или гомогенизированный эффект затухания. Численное моделирование процесса дефазировки является довольно сложным и требует большой вычислительной мощности. Действительно, если плотность воды 0,995 г/см3 при температуре 36,6 °С и молярная масса 18,01524 г/моль, то в вокселе (2 мм)3 находится 270 ∙ 1018 молекул или 4,441 ∙ 10–4 моль. Молекулы воды размером 3,2 Å испытывают 106 столкновений во время пробега 15 мкм. Огромное количество сталкивающихся молекул вместе с сильно неоднородной структурой среды серьезно препятствует прямому моделированию дискретного случайного блуждания. В этих условиях среднюю диффузию можно оценить либо с помощью обычных схем эффективной среды [4, 5], предполагая принципы суперпозиции (аппроксимация Максвелла – Гарнетта), либо используя специальную технику пространственного усреднения [6–9].
В этой работе мы исследуем связь между локальными и эффективными (глобальными) диффузионными транспортными свойствами и оцениваем чувствительность коэффициента диффузии к доминирующему набору микропараметров: внеклеточной объемной фракции, внеклеточной диффузии, диффузии миелина, аксонов, среднего размера аксона, среднего размера миелиновой оболочки, плотности внеклеточных протонов, плотности протонов миелина и плотности протонов аксона. Цель работы состоит в том, чтобы предоставить инструмент, охватывающий описание всего диапазона переменных и неоднородностей, флуктуирующих на иерархии пространственных масштабов. В следующем разделе объясняется структурная модель случайного гетерогенного материала, которая проиллюстрирована с помощью разбиения пространства на квадраты. Затем следуют результаты классической теории диффузии и ее применения к локальным оценкам диффузии. Наконец, реализация метода ренормализационной группы (РГ) вместе с его численными результатами обсуждаются в последнем разделе.
СТРУКТУРНАЯ МОДЕЛЬ МОЗГОВОГО БЕЛОГО ВЕЩЕСТВА
В этом разделе мы в общих чертах объясним основные понятия метода ренормализационной группы и применим их к задаче оценки эффективного коэффициента диффузии мозгового белого вещества, моделируемого разбиением на квадраты двумерного пространства.
В изображениях гистологического сечения биологической ткани (белого вещества головного мозга), полученных на электронном микроскопе можно отчетливо распознать два набора в структуре: миелиновые аксоны и внеклеточную гелеобразную матрицу [10–12]. Определение структурных и эффективных физических (диффузионных) параметров материала со случайно распределенными локальными свойствами – довольно сложная обратная задача. Поэтому обычно вводятся некоторые упрощающие исходную проблему предположения. Миелиновые аксоны можно рассматривать как цилиндры, которые статистически однородно и изотропно распределены в матрице флюида. Симметрия задачи позволяет разбить систему на ортогональные 1d продольные и 2d поперечные подпространства. Затем для математической оценки ортогональных признаков ткань может быть отображена на квадратные разбиения (рис.1а). Предположим, что черные квадратные блоки описывают свойства миелинованного аксона (ICS), а белые квадратные блоки относятся к матрице флюида (ECS). Таким образом, белое вещество представляется в виде пучка цилиндров (миелиновых аксонов), случайно внедренных в матрицу, и может рассматриваться как двухфазная среда, где одна из фаз (аксоны) обладает составными свойствами.
Анатомические срезы белого вещества показывают, что миелиновая оболочка нейронов имеет конечную толщину, и продольная структура аксона построена из следующих компонентов: миелиновая оболочка, которая характеризуется внутренними (Ra), и внешними (Rm) диаметрами, аксонной мембраной и узлами Ранвье (рис.1b). Диффузия через миелин затруднена двойными слоями липидов и отдельной оболочкой мембраны. При моделировании свойств белого вещества можно определить проницаемость нейрона, поскольку комбинированный эффект при определении усредненной диффузии включает влияние всех подструктур.
На рис.1c поперечная структура аксона изображена с заданными параметрами диффузии (De) и равновесной концентрацией молекул воды, ce для ECS; в то же время ICS характеризуется аксональным (Dа) и миелиновым (Dm) коэффициентами диффузии и концентрациями ca и cm соответственно.
Рассмотрим двухфазную структуру, заполняющую пространство регулярным разбиением, чьи отдельные белые и черные квадратные блоки представляют собой либо ECS, либо ICS (рис.2a). Такое разбиение может рассматриваться как набор примитивной ячейки Вигнера-Зейтца. Пусть общая объемная доля ECS (белые блоки на рис.2a) в композите будет равна p0. Затем часть белых блоков обладает верхней границей (U) диффузии. ICS занимает долю (1 – p0) и обладает нижней границей (L) диффузии. Пусть D^U и D^L – локальные диффузии чистых фаз U и L чьи свойства определяются линейными уравнениями состояния.
Локальное распределение вероятностей коэффициентов диффузии для этого разбиения:
, (2)
где δ(D^) – функция Дирака, и предполагается, что D^0U = D^U и D^0L = D^L. Эта разбиение однородно в масштабе отдельного блока белого квадрата и гетерогенна для блока черного квадрата, представляющего аксон, покрытый миелином. Вместо отдельных блоков репрезентативная группа смежных блоков может считаться основным строительным блоком разбиения. Такая группа, называемая перенормировочным кластером, может содержать любую возможную конфигурацию фаз и, следовательно, может быть однородной или гетерогенной в зависимости от того, содержит ли она одну фазу или обе фазы. На рис.2b приведены возможные конфигурации ренормируемого кластера для квадратного разбиения. Перенормировочные кластеры могут варьироваться по размеру, но обычно выбираются симметричными, чтобы удовлетворить изотропности разбиения на каждом этапе процесса перенормировки. Кластер, показанный на рис.2b, является самым маленьким, подходящим для квадратного разбиения. Число блоков в перенормировочном кластере Ld, где d – это размерность разбиения и L-масштаб длины кластера. Так как масштаб длины единичного блока разбиения произвольно принимается за единицу, то L – это длина масштабирования. Например, L для кластера, показанного на рис.2b, равно 2.
Целью процесса перенормировки является замена исходного гетерогенного разбиения эквивалентным однородным разбиением, эффективный коэффициент диффузии которого должен быть определен (рис.2b). Это достигается посредством операции последовательного усреднения в масштабе ренормгруппового кластера, так что на каждом шаге получается новое распределение вероятности. Предполагается, что морфология сохраняет свой беспорядок на каждом шаге, а распределение вероятности всегда состоит из ряда дельта-функций. После n перенормировочных шагов распределение получает вид Pn(D^(r→)).
В общем случае любое разбиение перенормируется путем замены каждого перенормировочного кластера блоком, коэффициент диффузии которого эквивалентен эффективному коэффициенту диффузии этого перенормировочного кластера. В результате этого преобразования средний размер кластера ξn+1 связан с ξn уравнением:
ξn+1 = ξn / L. (3)
То есть каждый шаг ренормализации уменьшает объем разбиения на 1/Ld поскольку Ld блок заменяется одним эквивалентным блоком. Если после n → ∞ шагов, ξn→∞ сходится к нулю, процесс завершается. На этом этапе разбиение является однородным и принимает одно из возможных D^Un→∞ или D ^Ln→∞ значений и Pn→∞(D^(r→)) приблизительно задается:
Pn→∞(D^(r→)) = δ(D^(r→) – D^eff), (4)
где D^eff = D ^U|n L →∞ – эффективный коэффициент диффузии композита, который определен теоретической моделью. Эта процедура перенормировки, часто называемая прореживанием, была впервые предложена для перколяции по узлам Рейнольдсом и др. [13].
Для осуществления процесса перенормировки необходимо определить все возможные расположения черных и белых фаз в перенормировочном кластере. Количество возможных конфигураций x фаз в кластере перенормировки Ld-блока равно xLd. Поскольку задание фаз для блоков является случайным, вероятность возникновения конкретной конфигурации фаз в перенормировочном кластере определяется произведением статистически независимых вероятностей возникновения фаз в блоках, составляющих кластер.
Предположим, что геометрия ECS позволяет перемещаться с одной стороны перенормировочного кластера на противоположную, не проникая в другую фазу. Тогда, применяя такое условие к разным реализациям двухкомпонентной системы, мы можем описать связывающий кластер условным конечномерным распределением R(p). Вероятность 1 – R(p) действительна для дополнительного случая.
Функция связности R(p) для кластера "2 × 2" РГ может быть выражена уравнением полинома
R(p) = p4 + 4p3 (1 – p) + 2p2 (1 – p)2. (5)
Уравнение 5 было построено по следующим правилам:
вероятность того, что четыре, три и два белых квадратных блока присутствуют в связывающем кластере равны p4, p3 (1 – p) и, p2 (1 – p)2, соответственно (рис.2b);
для описания связности противоположных сторон доступны только эти конфигурации;
вероятность связности вычисляется путем (табл.1): (a) умножения вероятности того, что rkn-ый блок существует на количество gk конфигураций блоков; (b) суммирования этих произведений по всем m конфигурациям (классам) блоков.
Вся процедура построения вероятности изображена на рис.2b. После упрощения полиномиального соотношения в уравнении (5) получим:
R(p) = 2p2 – p4. (6)
Если учесть свойства перенормируемости, то приходим к уравнению
pn+1 = 2p2n – p4n, (7)
где n – уровень кластера перенормируемости.
Корни (решения) или неподвижные точки алгебраического уравнения (7), которые принадлежат к физически значимому диапазону (0≤p*≤1), являются p* = {0,0.618,1}. Первый и третий корни представляют собой занятые и свободные решетки (полностью занятые или пустые пределы). Эти пределы называются "устойчивыми" или "однородными", поскольку рекурсивная процедура для решения уравнения (7) приводит к одному из этих пределов путем последовательного уменьшения фазы. В этих точках система проявляет трансляционную симметрию, а "устойчивые" неподвижные точки представляют собой "стоки". Вторая неподвижная точка p* = {0,618} обозначается как "неустойчивая" или "источник" потока параметра перенормировки pn.
Если концентрация фаз композита задана точно неустойчивой неподвижной точкой, то средний размер кластера бесконечно большой. Следовательно, для перенормирования такой гетерогенной системы требуется бесконечно большое количество шагов итераций. Таким образом, количество шагов, требуемых в процедуре перенормировки, увеличивается с увеличением размера кластера по мере приближения p0 к неустойчивой точке.
Функция связности R(p) и ее производная R'(p) изображены на рис.2d. Максимальное значение R'(p) находится в окрестности фиксированной точки p* = 0,618. Тонкие пунктирные линии, соединяющие биссектрису и функцию связности, – это перенормировка объемной доли pn согласно уравнению (7). В непосредственной близости p* = 0,618 система обладает свойствами само-подобия и может рассматриваться как стохастический фрактал. Вблизи неподвижной точки типичный размер линейного кластера (или корреляционная длина) ξ = Ln расходится как:
ξ ~ |p – p* |–ν, (8)
и показатели корреляционной функции могут быть оценены с помощью уравнения (7)
, (9)
если для процесса перенормирования принято правило
, (10)
сохраняющееся для каждого уровня масштаба. Критическая оценка экспоненты из уравнения (9) составляет приблизительно 1,64. Таким образом, диапазон масштабирования корреляционной длины гетерогенного кластера находится в пределах
b ≤ ξ(p) ≤ l, (11)
где l – макроскопический размер образца. В области p < p* корреляционная длина представляет собой типичный линейный размер "островов" ECS, окруженных пучками волокон, а когда p > p*, корреляционная длина представляет собой типичный линейный размер сгруппированных волоконных пучков, встроенных в матрицу флюида.
ЛОКАЛЬНЫЕ ДИФФУЗИОННЫЕ СВОЙСТВА
Поскольку мы рассматриваем диффузию D^(r→) в стационарном режиме без дополнительных истоков молекул воды, то локально в каждой точке гетерогенной среды закон Лапласа (консервативный и потенциальный)
∇→ D^(r→) ∇→c(r→) = 0 (12)
и закон Фика (линейное отношение состояния)
J→ (r→) = –D^ (r→) ∇→c(r→) (13)
выполняются вместе с условиями непрерывности Коши на границе двух фаз ∂Г1 и ∂Г2 для полноты. Для нормальной составляющей потока массы J→n→ (r→) на границе поверхности раздела фаз выполняется
J→n→ (r→) | ∂Г1= J→n→ (r→) | ∂Г2, (14)
и аналогично для концентрации справедливо
c(r→) | ∂Г1= c(r→) | ∂Г2. (15)
Для геометрии, описанной на рис.1b,c, граничные условия (уравнения (13, 14, 15)) имеют вид:
. (16)
Здесь мы предположили, что температура и химическая активность однородны по всему объекту, а движущей силой массопереноса является только функция концентрации. В реальном случае аксоны покрыты миелиновой оболочкой, прерывающейся узлами Ранвье. Обе эти компоненты играют роль в электрических свойствах аксона и передают нейронный сигнал. Миелин имеет очень плотную структуру, что приводит к незначительной диффузии (0,3 · 10–9 м2/с) по сравнению с таковой для экстра- или интрааксонного флюида. Диффузия в ECS затруднена из-за присутствия клеток глии и других структур и ниже, чем коэффициент диффузии в чистой воде. Мы будем считать, что присутствие глии в ECS не нарушает изотропность среды. В этом случае ECS можно охарактеризовать сферическим тензором:
. (17)
Цилиндрические волокна ICS обладают поперечной анизотропией, которая выражается в тензорной форме как:
. (18)
После локальной процедуры гомогенизации закон Фика определяется как
<J→(r→)> = –D^eff<∇→c(r→)>, (19)
где объемный поток и концентрация, соответственно, равны:
(20)
Если имеется симметрия относительно вращения вокруг одной из пространственных осей, например оси z, то можно показать, что эффективный тензор диффузии имеет две независимые компоненты:
. (21)
Такой тензор характеризуется поперечной изотропной симметрией и описывает усредненную геометрию цилиндрических волокон, внедренных в ECS.
Продольную компоненту тензора диффузии можно рассчитать как среднее по объему [4]:
. (22)
В формуле (22) f – это объемная доля аксона внутри примитивной ячейки Вигнера – Зейтца и эффективная концентрация ceff молекул воды рассчитывается соотношением [4]
. (23)
Поперечные свойства определяются как [4]:
, (24)
где
, (25)
а также
ε = сeDe, ε = сaDa, ε = сmDm.
МЕТОД РЕНОРМАЛИЗАЦИОННОЙ ГРУППЫ
В биологических тканях можно распознавать уединенные волокна, небольшие группы изолированных волокон и крупные пучки волокон, ограничивающих внеклеточный флюид. Такую систему трудно описать, ограничившись двумя различными масштабами, где диффузионные свойства изменяются быстро или медленно. Вместо этого система демонстрирует все диапазоны характерных размеров ξ неоднородностей, и они находятся в определенном интервале (уравнение (11)). Поэтому, используя модель "2 × 2" – РГ кластера, можно определить разнообразие диффузии во всех масштабах и оценить эффективный коэффициент поперечной диффузии.
Пусть каждый квадратный блок разбиения (рис.2а) представляет собой эффективный коэффициент микроскопической диффузии в локальной области. Структурированный набор состоит из двух типов квадратов с коэффициентами диффузии D111 (белый квадрат) и коэффициентов диффузии D112 (черный квадрат) соответственно. Вероятность появления белого квадрата равна p0, вероятность возникновения черного (1 – p0). Распределение фаз с двумя различными свойствами описывается обобщенной функцией плотности вероятности уравнения (2).
Предположим, что случайная биологическая ткань может быть смоделирована по схеме, изображенной на рис.2с. Если n – это итеративный шаг в построении иерархической решетки, то для каждого итерационного шага функция распределения плотности вероятности (уравнение (2)) имеет вид:
. (26)
Значение Def(U f, j) представляет собой усредненный эффективный коэффициент диффузии, рассчитанный для случая ECS. Такой эффективный коэффициент диффузии рассматривается как верхняя граница коэффициента диффузии. Для начального шага j = 1 эффективный коэффициент диффузии определяется Def(Uf, 1) = D111 , где D111 является внеклеточным коэффициентом диффузии флюида. Коэффициент Def(L f, j) представляет собой обратную ситуацию и, следовательно Def(Lf, 1) = D112 , где D112 – диффузия в аксонах.
Локально для гетерогенной ткани усредненные верхние и нижние границы коэффициента диффузии могут быть оценены с помощью модели Сена-Бассера (уравнения (24, 25)). Верхняя и нижняя границы могут быть представлены белым и черным квадратами на "2 × 2" – РГ кластере. Тогда эффективный коэффициент диффузии можно оценить с помощью рекурсивного алгоритма. Этот алгоритм ранее использовался для решения проблемы деградации в твердом и полимерном материалах с помощью методов Монте-Карло [14]. Верхняя граница эффективного коэффициента диффузии оценивается как [15]
, (27)
а нижняя граница имеет вид:
. (28)
В целом алгоритм ренормализационной группы состоит из следующих шагов:
Инициализировать верхнюю и нижнюю границы коэффициентов диффузии, а именно, Def(U f, j) и Def(L f, j).
Инициализировать вероятность pj появления белого квадрата на РГ кластере "2 × 2".
Оценить коэффициенты диффузии всего "2 × 2" РГ кластера в соответствии с моделью Сена-Бассера (СБ) для верхнего (Def(Uf, j + 1)) (уравнения 24, 25) и нижнего (Def(Lf, j + 1)) (уравнения 24, 25) пределов (табл.1).
Оценить вероятность R(pj) связности белых блоков (уравнение 7) на кластере "2 × 2".
Вернуться к шагу 1, если нормы |Def(Uf, j + 1) – Def(Uf, j)| и |Def(Lf, j + 1) – Def(Lf, j )| между эффективными коэффициентами диффузии текущей и предыдущей итераций больше заданной ошибки.
Остановить процедуру, если нормы |Def(Uf, j + 1) – Def(Uf, j)| и |Def(Lf, j + 1) – Def(Lf, j )| между эффективными коэффициентами диффузии текущей и предыдущей итерации будут меньше заданной ошибки. Полученный коэффициент Def(Uf(L), j + 1) = Deff(p) соответствует эффективной диффузии с p = p1. Таким образом, рассчитанная поперечная диффузия D11eff = Deff(p).
При построении структуры, описанной на рис.2b, c, перенормируются следующие три параметра (табл.1):
. (29)
Система (29) представляет собой систему нелинейных уравнений, образующих поток параметров в шестимерном пространстве.
ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ
Чтобы получить реалистичные коэффициенты диффузии в биологической ткани, мы используем параметры, экспериментально измеренные в [16, 17]. Эти параметры заданы в табл.2. Начальное значение верхней границы нормировалось Def(Uf, 1) / D111 = 1
(где D111 – ECS диффузия), а нижняя граница задавалась отношением Def(Lf, 1) / D111 на начальном итерационном этапе (j = 1). Значение Def(Lf, 1) было рассчитано в соответствии с моделью Сена-Бассера (уравнение (4)), где коэффициент ICS диффузии обозначен как D2 = D112 ). На рис.3а изображена эволюция перенормировки поперечного эффективного коэффициента диффузии D11eff во время j-шагов рекурсивной процедуры. На рисунке представлены верхний и нижний пределы диффузии вместе с их средним значением. Наблюдаемый коэффициент диффузии флуктуирует между этими границами. На небольших масштабах ξ > b флуктуации являются гигантскими, но, когда ξ >> b и ξ → L, система демонстрирует статистически однородные свойства. Это явление явно проявляется на рис.3а, b, где изображена сходимость верхней и нижней границы слева направо с растущим масштабом ξ. Видно, что при эволюции параметров можно выделить три отдельные зоны: I – линейная; II – нелинейная и III – стабильная. Режим I –
область неустойчивости (демонстрирующая масштабную симметрию), во время итерационного процесса рекурсивное значение перемещается к области притяжения III (область с трансляционной симметрией). Режим II можно рассматривать как зону сепаратрисы.
На рис.3а начальные доли ECS были p1 = 0,25 (сплошная линия) и p1 = 0,35 (пунктирная линия). Очевидно, что верхняя и нижняя границы коэффициента диффузии сходятся к некоторому значению (т.е. аттрактору), которое зависит от начальной доли ECS. Как можно видеть на рис.3а, б, большее значение эффективного коэффициента диффузии соответствует большему значению начальной доли ECS фазы. Для случая, когда p1 = 0,61 число итераций j = 11, это строит решетку, состоящую из ~811 связей. Такого количества связей достаточно, чтобы получить описание физических свойств сильно неупорядоченной и гетерогенной биологической ткани. Процедура итерации параметров останавливается в соответствии с заранее заданным условием , где была выбрана приемлемая норма ε = 10–3.
На рис.4а представлена зависимость нормированного коэффициента поперечной эффективной диффузии Deff, Dt, eff, Dl, eff от концентрации p ECS. Доля клеток при этом была (1 – p). Каждой точке кривой соответствует точка аттрактора после j-шагов процедуры рекурсии. Правая и левая асимптоты (p →{0,1}) сходятся к верхней и нижней границам приближения эффективной среды соответственно.
Промежуточная асимптота расположена между верхней и нижней границами эффективного коэффициента диффузии и демонстрирует характерную нелинейность. Эта нелинейность исчезает в пределе очень высокой проницаемости, где преобладает проницаемость над эффектами локализации (не показано на рисунке). Такая нелинейность обычна для многих диффузионных (транспортных) процессов и часто изображается в логарифмическом масштабе [19, 20]. Нелинейность максимальна вблизи пороговой концентрации p* = 0,618 для двумерного случая. В реалистичной ситуации в экспериментах объемная доля внеклеточной воды изменяется от 0,1 до 0,4. Как можно понять из моделирования в таком диапазоне дисперсии объемной доли, эффективный поперечный коэффициент диффузии более чувствителен к изменению проницаемости, а не к объемной доле ECS.
Для прикладных задач полезно знать чувствительность
S = |∂log(Deff) / ∂log(X)| (34)
эффективного коэффициента диффузии, вычисленного как след эффективного тензора диффузии в уравнении (21):
Deff = (2D11eff + D33eff ) / 3. (35)
Исходные значения (X в уравнении (34)) взяты из литературы, где они были экспериментально получены (табл.2, [2]). Результаты представлены на рис.4b. Анализ чувствительности эффективного коэффициента диффузии дает представление о том, как связаны различные микропараметры с развитием мозга или его патологическим состояниям. Аналогичная проблема рассмотрена в [18], где использовались численные эксперименты прямого броуновского случайного блуждания Монте-Карло. В обоих случаях подтверждается нелинейное поведение эффективной диффузии по сравнению с внеклеточной объемной концентрацией. Причиной такого явно выраженного нелинейного поведения является увеличение вероятности блокирования молекул воды с уменьшением внеклеточной объемной концентрации. С другой стороны, нейронная клетка может проявлять эффекты деградации, приводящие к повышению проницаемости; такой эффект уменьшает блокировку молекул воды. Это продемонстрировано на рис.4b.
Дизайн экспериментов по диффузии, который может подтвердить предложенную теорию, приведен в [19, 20]. Существенным моментом таких экспериментов является изменение шкалы времени в качестве аргумента диффузии, связанной с зондированием пространства. На коротких временах ожидается сильная флуктуация из-за неоднородности белого вещества головного мозга. Такие флуктуации уменьшаются на больших временах, что позволяет оценить эффективный коэффициент диффузии.
ВЫВОДЫ
В данной работе мы предложили итерационный алгоритм масштабного усреднения для оценки эффективной диффузии случайной гетерогенной биологической ткани. Диффузионные (транспортные) свойства изучались в широком диапазоне микроскопических и композиционных параметров биологической многокомпонентной гетерогенной ткани, состоящей из клеточных и внеклеточных доменов. В этих доменах молекулы характеризовались разными коэффициентами диффузии и плотностями (концентрациями), а именно внутри клеток и во внешней среде. Результаты предложенной итерационной схемы усреднения по масштабам образца были использованы для изучения влияния большого диапазона микроструктурных и композиционных параметров на эффективный коэффициент диффузии на итерационной решетке, состоящей из ~811 связей. Кроме того, была представлена самосогласованная схема моделирования для прогнозирования эффективного коэффициента диффузии на длительных временных масштабах. Эта схема расчетов позволила выявить зависимость эффективного коэффициента диффузии от отношения микроскопических коэффициентов диффузии фаз композита, их концентраций и проницаемости клеток.
ЛИТЕРАТУРА/REFERENCES
Callaghan P.T. Principles of Nuclear Magnetic Resonance Microscopy (Clarendon Press, Oxford) 1991.
Sen P.N., Basser P.J. Biophys. J., 89 (2005) 2927.
Szafer A., Zhong J., Gore J. Magn. Reson. Med., 33 (1995) 697.
Crank J. The mathematics of diffusion (Clarendon Press, Oxford) 1955.
Carslow H.S., Jager J.C. Conduction of Heat in Solids (Clarendon Press, Oxford) 1959.
Posnansky O., Guo J., Hirsch S., Papazoglou S., Braun J., Sack I. Fractal network dimension and viscoelastic powerlaw behavior: I. A modeling approach based on a coarse-graining procedure combined with shear oscillatory rheometry, Phys. Med. Biol., 57 (12) (2012) 4023.
Posnansky O. Viscoelastic properties of a hierarchical model of soft biological tissue: two-dimensional and three-dimensional cases, Journ. Stat. Phys., 164 (5) (2016) 1043.
Posnansky O. Dynamic characteristics of the effective susceptibility function in random three-component system, Physica A: Stat. Mech. and its Applications 473 (2017) 18.
Posnansky O. Tensor of effective susceptibility in random magnetic composites: application to two-dimensional and three-dimensional cases, Journal of molecular structure, 1160 (2018), 293.
Rollenhagen A., Lübke J.H.R. Cell Tissue Res 326 (2006) 221.
Harsan L. A., Poulet P. et al. J. Neuroscience Res., 85 (2007) 935.
Rushton W. J. Phisiol., 115 (1951) 101.
Reynolds P. J., Stanley H. E., Klein W. Phys. Rev. B, 21 (1980) 1223.
Poznansky O.P., Novikov V.V. J. Polym. Eng., 19 (1999) 223.
Hashin Z. Theory of composite materials (Pergamon Press, New York).
Latour L.L., Svoboda K., Mitra P.P., Sotak C.H. PNAS, 91 (1994) 1229.
Neeb H., Zilles K., Shah N. Neuroimage 31 (2006) 1156–1168.
Hall M., Alexander D. ISMRM 2008, 1779.
Posnansky O., Kupriyanova Y., Shah J. ESMRB 2008, 633.
Posnansky O. Effective diffusion in random composites measured by NMR, Journal of magnetism and magnetic materials 466 (2018) 92.
Posnansky O.P. Dynamic diffusion measured in binary composites measured by NMR, NANOINDUSTRY, 11, 7–8 (2018) 22.
Получено: 24.12.2018 г.
Мы разработали метод ренормализационной группы, чтобы исследовать влияние большого разброса геометрических и физических микропараметров на эффективный коэффициент диффузии. Основным элементом этого метода является перенормировка структурных параметров от микро- до миллиметрового масштаба случайной геометрии многокомпонентного белого вещества головного мозга. Наш подход улучшает модель Сена-Бассера для белого вещества мозга [Sen P. N., Basser P. J., Biophys. J., 89 (2005) 2927] благодаря учету влияния беспорядка и позволяет количественно исследовать чувствительность эффективного коэффициента диффузии к вариациям доминирующего набора микропараметров.
We have developed a renormalization group method in order to explore the influence of a large range of geometrical and physical microparameters on the effective diffusion coefficient. The main element of this method comprises of the renormalization of structural parameters from the micro- to the millimeter scale of the random multi-component white matter of the brain. Our approach enhances the Sen-Basser model of the brain white matter [Sen P. N., Basser P. J., Biophys. J., 89 (2005) 2927] by taking the influence of disorder into consideration and allowing quantitative investigation of the sensitivity of the effective diffusion coefficient to the variations of the dominant set of microparameters.
ВВЕДЕНИЕ
Белое вещество головного мозга можно разделить на две области: внутриклеточное пространство (ICS) и внеклеточное пространство (ECS). ICS состоит из плотно упакованного композита аксонов, глии и их клеточных образований, а ECS представляет собой желеобразную матрицу, которая разделяет клетки биологической ткани. ECS поддерживает основные функции мозга, включая межклеточный обмен, транспортировку питательных веществ и доставку лекарств. ECS представляет только ~ 20% объема здорового головного мозга, но может резко сокращаться и динамически расширяться при изменении объема клеток головного мозга при стимуляции нейронов или при введении осмотических агентов. Изменения в размере и геометрии ECS сопровождаются многими патологическими отклонениями, указывающими на ишемию, воспаление и прогрессирующую опухоль.
Принимая во внимание, что мозг представляет собой среду, состоящую из случайно расположенных доменов разных фаз, ее эффективные макроскопические свойства могут быть описаны методами, разработанными для случайных гетерогенных материалов. Используя современные методы магнитно-резонансной томографии (МРТ), можно получить уникальную информацию об организации белого вещества, нейронных путей и их визуализации. Макроскопические свойства можно оценить, используя усреднение микроструктуры. Например, объемная доля ECS (p) и эффективный коэффициент диффузии (Deff) являются двумя основными макроскопическими параметрами, обычно используемыми для описания микроструктуры мозга. Таким образом, проблему можно сформулировать следующим образом: используя методы МРТ, можно оценить макроскопический коэффициент диффузии биологической ткани. Кроме того, можно оценить чувствительность коэффициента диффузии к объемной доле включенных фаз, отношения микроскопических коэффициентов диффузии, неупорядоченной геометрии доменов и концентрации протонов в различных областях.
Понимание влияния всех этих параметров на диффузию важно для клинических приложений, поскольку их существенное отличие от средних значений является признаком патологических состояний мозга, что дает возможность обнаруживать заболевания на самых ранних стадиях.
В экспериментах с МРТ диффузия измеряется как часть ослабления сигнала из-за декогеренции в фазовых сдвигах спинов молекул воды, вызванных поступательным движением в присутствии внешнего линейного магнитного поля. Диффузия задает фундаментальное ограничение на разрешение изображений в МРТ, но, с другой стороны, может обеспечить контрастный механизм, который позволяет визуализировать спектр молекулярного смещения. При этом предполагается идеальная система градиентных катушек и компенсированные локальные магнитные неоднородности для несмещенной оценки параметров диффузии. В процессе затухания сигнала релаксация намагниченности характеризуется продольным (T1) и поперечным (T2) временами. Релаксация может быть разной для разных доменов из-за изменения локальной намагниченности. Тем не менее при диффузии за длительный период время затухания сигнала может быть описано одной экспоненциальной кривой с заданным коэффициентом диффузии:
S/S0 = exp(–b∙Deff), (1)
где b = γ2G2δ2 (Δ – δ/3) – коэффициент, определенный по формуле Стежскал-Таннера [1, 21] (γ является гиромагнитной постоянной, S0 константа, связанная с временами релаксаций, Δ является временем диффузии и G является градиентом магнитного поля, характеризуемым длительностью δ).
В предпринятых попытках моделирования микроструктуры мозга было показано, что анизотропные свойства коэффициента диффузии хорошо согласуются с моделями ограниченной диффузии, особенно в режиме длительного временного периода. Такая модель диффузии хорошо описывает анизотропию и отслеживает волокна в белом веществе головного мозга. Обычный подход, применяемый для моделирования гетерогенной природы ткани головного мозга, заключается в изучении волокон, представленных в виде цилиндров, расположенных в различных типах симметричных решеток, таких как квадрат или шестиугольник, и использовании аппроксимации эффективной среды Максвелла – Гарнетта [2] или последовательно-параллельного приближения [3]. В таких моделях трудно вводить случайность в структуре, которая очень естественна в биологических материалах, но совершенно противоположна кристаллам твердого тела. Используя предложенный в этой статье подход в диффузии, можно исследовать неупорядоченные биологические материалы и учитывать влияние сложной геометрии и топологии на эффективные транспортные (диффузионные) свойства ткани.
Типичная диффузионно-взвешенная МРТ-последовательность ЭПП (эхо-планарная последовательность) позволяет получать информацию о дефазировке спинов молекулы воды в течение Td = 10 ÷ 60 мс. Коэффициент самодиффузии чистой воды при температуре 36,6 °C равен 2,94 ∙ 10–9 м2/с. Мы считаем, что присутствие глии-клеток во внеклеточном пространстве согласуется с ограниченным коэффициентом диффузии, который ниже на 30%, чем у свободной воды при температуре тела и примерно равен 2 ∙ 10–9 м2/с. Кроме того, внутриклеточная диффузия сильно зависит от нейрофиламентов и микротрубочек и, следовательно, меньше на 75% диффузии воды и оценивается как 0,75 ∙ 10–9 м2/с. Используя уравнение диффузии Эйнштейна (<x2> = 2DTd), можно вычислить среднее смещение молекулы воды (<x2>), которое равно 5 ÷ 15 мкм. Диаметр волокон примерно равен 3 мкм. Среднее расстояние между волокнами случайно и находится в интервале 1 ÷ 18 мкм. Легко понять, что молекула воды зондирует все составляющие компоненты биологической ткани, а МРТ-сигнал дает усредненный или гомогенизированный эффект затухания. Численное моделирование процесса дефазировки является довольно сложным и требует большой вычислительной мощности. Действительно, если плотность воды 0,995 г/см3 при температуре 36,6 °С и молярная масса 18,01524 г/моль, то в вокселе (2 мм)3 находится 270 ∙ 1018 молекул или 4,441 ∙ 10–4 моль. Молекулы воды размером 3,2 Å испытывают 106 столкновений во время пробега 15 мкм. Огромное количество сталкивающихся молекул вместе с сильно неоднородной структурой среды серьезно препятствует прямому моделированию дискретного случайного блуждания. В этих условиях среднюю диффузию можно оценить либо с помощью обычных схем эффективной среды [4, 5], предполагая принципы суперпозиции (аппроксимация Максвелла – Гарнетта), либо используя специальную технику пространственного усреднения [6–9].
В этой работе мы исследуем связь между локальными и эффективными (глобальными) диффузионными транспортными свойствами и оцениваем чувствительность коэффициента диффузии к доминирующему набору микропараметров: внеклеточной объемной фракции, внеклеточной диффузии, диффузии миелина, аксонов, среднего размера аксона, среднего размера миелиновой оболочки, плотности внеклеточных протонов, плотности протонов миелина и плотности протонов аксона. Цель работы состоит в том, чтобы предоставить инструмент, охватывающий описание всего диапазона переменных и неоднородностей, флуктуирующих на иерархии пространственных масштабов. В следующем разделе объясняется структурная модель случайного гетерогенного материала, которая проиллюстрирована с помощью разбиения пространства на квадраты. Затем следуют результаты классической теории диффузии и ее применения к локальным оценкам диффузии. Наконец, реализация метода ренормализационной группы (РГ) вместе с его численными результатами обсуждаются в последнем разделе.
СТРУКТУРНАЯ МОДЕЛЬ МОЗГОВОГО БЕЛОГО ВЕЩЕСТВА
В этом разделе мы в общих чертах объясним основные понятия метода ренормализационной группы и применим их к задаче оценки эффективного коэффициента диффузии мозгового белого вещества, моделируемого разбиением на квадраты двумерного пространства.
В изображениях гистологического сечения биологической ткани (белого вещества головного мозга), полученных на электронном микроскопе можно отчетливо распознать два набора в структуре: миелиновые аксоны и внеклеточную гелеобразную матрицу [10–12]. Определение структурных и эффективных физических (диффузионных) параметров материала со случайно распределенными локальными свойствами – довольно сложная обратная задача. Поэтому обычно вводятся некоторые упрощающие исходную проблему предположения. Миелиновые аксоны можно рассматривать как цилиндры, которые статистически однородно и изотропно распределены в матрице флюида. Симметрия задачи позволяет разбить систему на ортогональные 1d продольные и 2d поперечные подпространства. Затем для математической оценки ортогональных признаков ткань может быть отображена на квадратные разбиения (рис.1а). Предположим, что черные квадратные блоки описывают свойства миелинованного аксона (ICS), а белые квадратные блоки относятся к матрице флюида (ECS). Таким образом, белое вещество представляется в виде пучка цилиндров (миелиновых аксонов), случайно внедренных в матрицу, и может рассматриваться как двухфазная среда, где одна из фаз (аксоны) обладает составными свойствами.
Анатомические срезы белого вещества показывают, что миелиновая оболочка нейронов имеет конечную толщину, и продольная структура аксона построена из следующих компонентов: миелиновая оболочка, которая характеризуется внутренними (Ra), и внешними (Rm) диаметрами, аксонной мембраной и узлами Ранвье (рис.1b). Диффузия через миелин затруднена двойными слоями липидов и отдельной оболочкой мембраны. При моделировании свойств белого вещества можно определить проницаемость нейрона, поскольку комбинированный эффект при определении усредненной диффузии включает влияние всех подструктур.
На рис.1c поперечная структура аксона изображена с заданными параметрами диффузии (De) и равновесной концентрацией молекул воды, ce для ECS; в то же время ICS характеризуется аксональным (Dа) и миелиновым (Dm) коэффициентами диффузии и концентрациями ca и cm соответственно.
Рассмотрим двухфазную структуру, заполняющую пространство регулярным разбиением, чьи отдельные белые и черные квадратные блоки представляют собой либо ECS, либо ICS (рис.2a). Такое разбиение может рассматриваться как набор примитивной ячейки Вигнера-Зейтца. Пусть общая объемная доля ECS (белые блоки на рис.2a) в композите будет равна p0. Затем часть белых блоков обладает верхней границей (U) диффузии. ICS занимает долю (1 – p0) и обладает нижней границей (L) диффузии. Пусть D^U и D^L – локальные диффузии чистых фаз U и L чьи свойства определяются линейными уравнениями состояния.
Локальное распределение вероятностей коэффициентов диффузии для этого разбиения:
, (2)
где δ(D^) – функция Дирака, и предполагается, что D^0U = D^U и D^0L = D^L. Эта разбиение однородно в масштабе отдельного блока белого квадрата и гетерогенна для блока черного квадрата, представляющего аксон, покрытый миелином. Вместо отдельных блоков репрезентативная группа смежных блоков может считаться основным строительным блоком разбиения. Такая группа, называемая перенормировочным кластером, может содержать любую возможную конфигурацию фаз и, следовательно, может быть однородной или гетерогенной в зависимости от того, содержит ли она одну фазу или обе фазы. На рис.2b приведены возможные конфигурации ренормируемого кластера для квадратного разбиения. Перенормировочные кластеры могут варьироваться по размеру, но обычно выбираются симметричными, чтобы удовлетворить изотропности разбиения на каждом этапе процесса перенормировки. Кластер, показанный на рис.2b, является самым маленьким, подходящим для квадратного разбиения. Число блоков в перенормировочном кластере Ld, где d – это размерность разбиения и L-масштаб длины кластера. Так как масштаб длины единичного блока разбиения произвольно принимается за единицу, то L – это длина масштабирования. Например, L для кластера, показанного на рис.2b, равно 2.
Целью процесса перенормировки является замена исходного гетерогенного разбиения эквивалентным однородным разбиением, эффективный коэффициент диффузии которого должен быть определен (рис.2b). Это достигается посредством операции последовательного усреднения в масштабе ренормгруппового кластера, так что на каждом шаге получается новое распределение вероятности. Предполагается, что морфология сохраняет свой беспорядок на каждом шаге, а распределение вероятности всегда состоит из ряда дельта-функций. После n перенормировочных шагов распределение получает вид Pn(D^(r→)).
В общем случае любое разбиение перенормируется путем замены каждого перенормировочного кластера блоком, коэффициент диффузии которого эквивалентен эффективному коэффициенту диффузии этого перенормировочного кластера. В результате этого преобразования средний размер кластера ξn+1 связан с ξn уравнением:
ξn+1 = ξn / L. (3)
То есть каждый шаг ренормализации уменьшает объем разбиения на 1/Ld поскольку Ld блок заменяется одним эквивалентным блоком. Если после n → ∞ шагов, ξn→∞ сходится к нулю, процесс завершается. На этом этапе разбиение является однородным и принимает одно из возможных D^Un→∞ или D ^Ln→∞ значений и Pn→∞(D^(r→)) приблизительно задается:
Pn→∞(D^(r→)) = δ(D^(r→) – D^eff), (4)
где D^eff = D ^U|n L →∞ – эффективный коэффициент диффузии композита, который определен теоретической моделью. Эта процедура перенормировки, часто называемая прореживанием, была впервые предложена для перколяции по узлам Рейнольдсом и др. [13].
Для осуществления процесса перенормировки необходимо определить все возможные расположения черных и белых фаз в перенормировочном кластере. Количество возможных конфигураций x фаз в кластере перенормировки Ld-блока равно xLd. Поскольку задание фаз для блоков является случайным, вероятность возникновения конкретной конфигурации фаз в перенормировочном кластере определяется произведением статистически независимых вероятностей возникновения фаз в блоках, составляющих кластер.
Предположим, что геометрия ECS позволяет перемещаться с одной стороны перенормировочного кластера на противоположную, не проникая в другую фазу. Тогда, применяя такое условие к разным реализациям двухкомпонентной системы, мы можем описать связывающий кластер условным конечномерным распределением R(p). Вероятность 1 – R(p) действительна для дополнительного случая.
Функция связности R(p) для кластера "2 × 2" РГ может быть выражена уравнением полинома
R(p) = p4 + 4p3 (1 – p) + 2p2 (1 – p)2. (5)
Уравнение 5 было построено по следующим правилам:
вероятность того, что четыре, три и два белых квадратных блока присутствуют в связывающем кластере равны p4, p3 (1 – p) и, p2 (1 – p)2, соответственно (рис.2b);
для описания связности противоположных сторон доступны только эти конфигурации;
вероятность связности вычисляется путем (табл.1): (a) умножения вероятности того, что rkn-ый блок существует на количество gk конфигураций блоков; (b) суммирования этих произведений по всем m конфигурациям (классам) блоков.
Вся процедура построения вероятности изображена на рис.2b. После упрощения полиномиального соотношения в уравнении (5) получим:
R(p) = 2p2 – p4. (6)
Если учесть свойства перенормируемости, то приходим к уравнению
pn+1 = 2p2n – p4n, (7)
где n – уровень кластера перенормируемости.
Корни (решения) или неподвижные точки алгебраического уравнения (7), которые принадлежат к физически значимому диапазону (0≤p*≤1), являются p* = {0,0.618,1}. Первый и третий корни представляют собой занятые и свободные решетки (полностью занятые или пустые пределы). Эти пределы называются "устойчивыми" или "однородными", поскольку рекурсивная процедура для решения уравнения (7) приводит к одному из этих пределов путем последовательного уменьшения фазы. В этих точках система проявляет трансляционную симметрию, а "устойчивые" неподвижные точки представляют собой "стоки". Вторая неподвижная точка p* = {0,618} обозначается как "неустойчивая" или "источник" потока параметра перенормировки pn.
Если концентрация фаз композита задана точно неустойчивой неподвижной точкой, то средний размер кластера бесконечно большой. Следовательно, для перенормирования такой гетерогенной системы требуется бесконечно большое количество шагов итераций. Таким образом, количество шагов, требуемых в процедуре перенормировки, увеличивается с увеличением размера кластера по мере приближения p0 к неустойчивой точке.
Функция связности R(p) и ее производная R'(p) изображены на рис.2d. Максимальное значение R'(p) находится в окрестности фиксированной точки p* = 0,618. Тонкие пунктирные линии, соединяющие биссектрису и функцию связности, – это перенормировка объемной доли pn согласно уравнению (7). В непосредственной близости p* = 0,618 система обладает свойствами само-подобия и может рассматриваться как стохастический фрактал. Вблизи неподвижной точки типичный размер линейного кластера (или корреляционная длина) ξ = Ln расходится как:
ξ ~ |p – p* |–ν, (8)
и показатели корреляционной функции могут быть оценены с помощью уравнения (7)
, (9)
если для процесса перенормирования принято правило
, (10)
сохраняющееся для каждого уровня масштаба. Критическая оценка экспоненты из уравнения (9) составляет приблизительно 1,64. Таким образом, диапазон масштабирования корреляционной длины гетерогенного кластера находится в пределах
b ≤ ξ(p) ≤ l, (11)
где l – макроскопический размер образца. В области p < p* корреляционная длина представляет собой типичный линейный размер "островов" ECS, окруженных пучками волокон, а когда p > p*, корреляционная длина представляет собой типичный линейный размер сгруппированных волоконных пучков, встроенных в матрицу флюида.
ЛОКАЛЬНЫЕ ДИФФУЗИОННЫЕ СВОЙСТВА
Поскольку мы рассматриваем диффузию D^(r→) в стационарном режиме без дополнительных истоков молекул воды, то локально в каждой точке гетерогенной среды закон Лапласа (консервативный и потенциальный)
∇→ D^(r→) ∇→c(r→) = 0 (12)
и закон Фика (линейное отношение состояния)
J→ (r→) = –D^ (r→) ∇→c(r→) (13)
выполняются вместе с условиями непрерывности Коши на границе двух фаз ∂Г1 и ∂Г2 для полноты. Для нормальной составляющей потока массы J→n→ (r→) на границе поверхности раздела фаз выполняется
J→n→ (r→) | ∂Г1= J→n→ (r→) | ∂Г2, (14)
и аналогично для концентрации справедливо
c(r→) | ∂Г1= c(r→) | ∂Г2. (15)
Для геометрии, описанной на рис.1b,c, граничные условия (уравнения (13, 14, 15)) имеют вид:
. (16)
Здесь мы предположили, что температура и химическая активность однородны по всему объекту, а движущей силой массопереноса является только функция концентрации. В реальном случае аксоны покрыты миелиновой оболочкой, прерывающейся узлами Ранвье. Обе эти компоненты играют роль в электрических свойствах аксона и передают нейронный сигнал. Миелин имеет очень плотную структуру, что приводит к незначительной диффузии (0,3 · 10–9 м2/с) по сравнению с таковой для экстра- или интрааксонного флюида. Диффузия в ECS затруднена из-за присутствия клеток глии и других структур и ниже, чем коэффициент диффузии в чистой воде. Мы будем считать, что присутствие глии в ECS не нарушает изотропность среды. В этом случае ECS можно охарактеризовать сферическим тензором:
. (17)
Цилиндрические волокна ICS обладают поперечной анизотропией, которая выражается в тензорной форме как:
. (18)
После локальной процедуры гомогенизации закон Фика определяется как
<J→(r→)> = –D^eff<∇→c(r→)>, (19)
где объемный поток и концентрация, соответственно, равны:
(20)
Если имеется симметрия относительно вращения вокруг одной из пространственных осей, например оси z, то можно показать, что эффективный тензор диффузии имеет две независимые компоненты:
. (21)
Такой тензор характеризуется поперечной изотропной симметрией и описывает усредненную геометрию цилиндрических волокон, внедренных в ECS.
Продольную компоненту тензора диффузии можно рассчитать как среднее по объему [4]:
. (22)
В формуле (22) f – это объемная доля аксона внутри примитивной ячейки Вигнера – Зейтца и эффективная концентрация ceff молекул воды рассчитывается соотношением [4]
. (23)
Поперечные свойства определяются как [4]:
, (24)
где
, (25)
а также
ε = сeDe, ε = сaDa, ε = сmDm.
МЕТОД РЕНОРМАЛИЗАЦИОННОЙ ГРУППЫ
В биологических тканях можно распознавать уединенные волокна, небольшие группы изолированных волокон и крупные пучки волокон, ограничивающих внеклеточный флюид. Такую систему трудно описать, ограничившись двумя различными масштабами, где диффузионные свойства изменяются быстро или медленно. Вместо этого система демонстрирует все диапазоны характерных размеров ξ неоднородностей, и они находятся в определенном интервале (уравнение (11)). Поэтому, используя модель "2 × 2" – РГ кластера, можно определить разнообразие диффузии во всех масштабах и оценить эффективный коэффициент поперечной диффузии.
Пусть каждый квадратный блок разбиения (рис.2а) представляет собой эффективный коэффициент микроскопической диффузии в локальной области. Структурированный набор состоит из двух типов квадратов с коэффициентами диффузии D111 (белый квадрат) и коэффициентов диффузии D112 (черный квадрат) соответственно. Вероятность появления белого квадрата равна p0, вероятность возникновения черного (1 – p0). Распределение фаз с двумя различными свойствами описывается обобщенной функцией плотности вероятности уравнения (2).
Предположим, что случайная биологическая ткань может быть смоделирована по схеме, изображенной на рис.2с. Если n – это итеративный шаг в построении иерархической решетки, то для каждого итерационного шага функция распределения плотности вероятности (уравнение (2)) имеет вид:
. (26)
Значение Def(U f, j) представляет собой усредненный эффективный коэффициент диффузии, рассчитанный для случая ECS. Такой эффективный коэффициент диффузии рассматривается как верхняя граница коэффициента диффузии. Для начального шага j = 1 эффективный коэффициент диффузии определяется Def(Uf, 1) = D111 , где D111 является внеклеточным коэффициентом диффузии флюида. Коэффициент Def(L f, j) представляет собой обратную ситуацию и, следовательно Def(Lf, 1) = D112 , где D112 – диффузия в аксонах.
Локально для гетерогенной ткани усредненные верхние и нижние границы коэффициента диффузии могут быть оценены с помощью модели Сена-Бассера (уравнения (24, 25)). Верхняя и нижняя границы могут быть представлены белым и черным квадратами на "2 × 2" – РГ кластере. Тогда эффективный коэффициент диффузии можно оценить с помощью рекурсивного алгоритма. Этот алгоритм ранее использовался для решения проблемы деградации в твердом и полимерном материалах с помощью методов Монте-Карло [14]. Верхняя граница эффективного коэффициента диффузии оценивается как [15]
, (27)
а нижняя граница имеет вид:
. (28)
В целом алгоритм ренормализационной группы состоит из следующих шагов:
Инициализировать верхнюю и нижнюю границы коэффициентов диффузии, а именно, Def(U f, j) и Def(L f, j).
Инициализировать вероятность pj появления белого квадрата на РГ кластере "2 × 2".
Оценить коэффициенты диффузии всего "2 × 2" РГ кластера в соответствии с моделью Сена-Бассера (СБ) для верхнего (Def(Uf, j + 1)) (уравнения 24, 25) и нижнего (Def(Lf, j + 1)) (уравнения 24, 25) пределов (табл.1).
Оценить вероятность R(pj) связности белых блоков (уравнение 7) на кластере "2 × 2".
Вернуться к шагу 1, если нормы |Def(Uf, j + 1) – Def(Uf, j)| и |Def(Lf, j + 1) – Def(Lf, j )| между эффективными коэффициентами диффузии текущей и предыдущей итераций больше заданной ошибки.
Остановить процедуру, если нормы |Def(Uf, j + 1) – Def(Uf, j)| и |Def(Lf, j + 1) – Def(Lf, j )| между эффективными коэффициентами диффузии текущей и предыдущей итерации будут меньше заданной ошибки. Полученный коэффициент Def(Uf(L), j + 1) = Deff(p) соответствует эффективной диффузии с p = p1. Таким образом, рассчитанная поперечная диффузия D11eff = Deff(p).
При построении структуры, описанной на рис.2b, c, перенормируются следующие три параметра (табл.1):
. (29)
Система (29) представляет собой систему нелинейных уравнений, образующих поток параметров в шестимерном пространстве.
ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ
Чтобы получить реалистичные коэффициенты диффузии в биологической ткани, мы используем параметры, экспериментально измеренные в [16, 17]. Эти параметры заданы в табл.2. Начальное значение верхней границы нормировалось Def(Uf, 1) / D111 = 1
(где D111 – ECS диффузия), а нижняя граница задавалась отношением Def(Lf, 1) / D111 на начальном итерационном этапе (j = 1). Значение Def(Lf, 1) было рассчитано в соответствии с моделью Сена-Бассера (уравнение (4)), где коэффициент ICS диффузии обозначен как D2 = D112 ). На рис.3а изображена эволюция перенормировки поперечного эффективного коэффициента диффузии D11eff во время j-шагов рекурсивной процедуры. На рисунке представлены верхний и нижний пределы диффузии вместе с их средним значением. Наблюдаемый коэффициент диффузии флуктуирует между этими границами. На небольших масштабах ξ > b флуктуации являются гигантскими, но, когда ξ >> b и ξ → L, система демонстрирует статистически однородные свойства. Это явление явно проявляется на рис.3а, b, где изображена сходимость верхней и нижней границы слева направо с растущим масштабом ξ. Видно, что при эволюции параметров можно выделить три отдельные зоны: I – линейная; II – нелинейная и III – стабильная. Режим I –
область неустойчивости (демонстрирующая масштабную симметрию), во время итерационного процесса рекурсивное значение перемещается к области притяжения III (область с трансляционной симметрией). Режим II можно рассматривать как зону сепаратрисы.
На рис.3а начальные доли ECS были p1 = 0,25 (сплошная линия) и p1 = 0,35 (пунктирная линия). Очевидно, что верхняя и нижняя границы коэффициента диффузии сходятся к некоторому значению (т.е. аттрактору), которое зависит от начальной доли ECS. Как можно видеть на рис.3а, б, большее значение эффективного коэффициента диффузии соответствует большему значению начальной доли ECS фазы. Для случая, когда p1 = 0,61 число итераций j = 11, это строит решетку, состоящую из ~811 связей. Такого количества связей достаточно, чтобы получить описание физических свойств сильно неупорядоченной и гетерогенной биологической ткани. Процедура итерации параметров останавливается в соответствии с заранее заданным условием , где была выбрана приемлемая норма ε = 10–3.
На рис.4а представлена зависимость нормированного коэффициента поперечной эффективной диффузии Deff, Dt, eff, Dl, eff от концентрации p ECS. Доля клеток при этом была (1 – p). Каждой точке кривой соответствует точка аттрактора после j-шагов процедуры рекурсии. Правая и левая асимптоты (p →{0,1}) сходятся к верхней и нижней границам приближения эффективной среды соответственно.
Промежуточная асимптота расположена между верхней и нижней границами эффективного коэффициента диффузии и демонстрирует характерную нелинейность. Эта нелинейность исчезает в пределе очень высокой проницаемости, где преобладает проницаемость над эффектами локализации (не показано на рисунке). Такая нелинейность обычна для многих диффузионных (транспортных) процессов и часто изображается в логарифмическом масштабе [19, 20]. Нелинейность максимальна вблизи пороговой концентрации p* = 0,618 для двумерного случая. В реалистичной ситуации в экспериментах объемная доля внеклеточной воды изменяется от 0,1 до 0,4. Как можно понять из моделирования в таком диапазоне дисперсии объемной доли, эффективный поперечный коэффициент диффузии более чувствителен к изменению проницаемости, а не к объемной доле ECS.
Для прикладных задач полезно знать чувствительность
S = |∂log(Deff) / ∂log(X)| (34)
эффективного коэффициента диффузии, вычисленного как след эффективного тензора диффузии в уравнении (21):
Deff = (2D11eff + D33eff ) / 3. (35)
Исходные значения (X в уравнении (34)) взяты из литературы, где они были экспериментально получены (табл.2, [2]). Результаты представлены на рис.4b. Анализ чувствительности эффективного коэффициента диффузии дает представление о том, как связаны различные микропараметры с развитием мозга или его патологическим состояниям. Аналогичная проблема рассмотрена в [18], где использовались численные эксперименты прямого броуновского случайного блуждания Монте-Карло. В обоих случаях подтверждается нелинейное поведение эффективной диффузии по сравнению с внеклеточной объемной концентрацией. Причиной такого явно выраженного нелинейного поведения является увеличение вероятности блокирования молекул воды с уменьшением внеклеточной объемной концентрации. С другой стороны, нейронная клетка может проявлять эффекты деградации, приводящие к повышению проницаемости; такой эффект уменьшает блокировку молекул воды. Это продемонстрировано на рис.4b.
Дизайн экспериментов по диффузии, который может подтвердить предложенную теорию, приведен в [19, 20]. Существенным моментом таких экспериментов является изменение шкалы времени в качестве аргумента диффузии, связанной с зондированием пространства. На коротких временах ожидается сильная флуктуация из-за неоднородности белого вещества головного мозга. Такие флуктуации уменьшаются на больших временах, что позволяет оценить эффективный коэффициент диффузии.
ВЫВОДЫ
В данной работе мы предложили итерационный алгоритм масштабного усреднения для оценки эффективной диффузии случайной гетерогенной биологической ткани. Диффузионные (транспортные) свойства изучались в широком диапазоне микроскопических и композиционных параметров биологической многокомпонентной гетерогенной ткани, состоящей из клеточных и внеклеточных доменов. В этих доменах молекулы характеризовались разными коэффициентами диффузии и плотностями (концентрациями), а именно внутри клеток и во внешней среде. Результаты предложенной итерационной схемы усреднения по масштабам образца были использованы для изучения влияния большого диапазона микроструктурных и композиционных параметров на эффективный коэффициент диффузии на итерационной решетке, состоящей из ~811 связей. Кроме того, была представлена самосогласованная схема моделирования для прогнозирования эффективного коэффициента диффузии на длительных временных масштабах. Эта схема расчетов позволила выявить зависимость эффективного коэффициента диффузии от отношения микроскопических коэффициентов диффузии фаз композита, их концентраций и проницаемости клеток.
ЛИТЕРАТУРА/REFERENCES
Callaghan P.T. Principles of Nuclear Magnetic Resonance Microscopy (Clarendon Press, Oxford) 1991.
Sen P.N., Basser P.J. Biophys. J., 89 (2005) 2927.
Szafer A., Zhong J., Gore J. Magn. Reson. Med., 33 (1995) 697.
Crank J. The mathematics of diffusion (Clarendon Press, Oxford) 1955.
Carslow H.S., Jager J.C. Conduction of Heat in Solids (Clarendon Press, Oxford) 1959.
Posnansky O., Guo J., Hirsch S., Papazoglou S., Braun J., Sack I. Fractal network dimension and viscoelastic powerlaw behavior: I. A modeling approach based on a coarse-graining procedure combined with shear oscillatory rheometry, Phys. Med. Biol., 57 (12) (2012) 4023.
Posnansky O. Viscoelastic properties of a hierarchical model of soft biological tissue: two-dimensional and three-dimensional cases, Journ. Stat. Phys., 164 (5) (2016) 1043.
Posnansky O. Dynamic characteristics of the effective susceptibility function in random three-component system, Physica A: Stat. Mech. and its Applications 473 (2017) 18.
Posnansky O. Tensor of effective susceptibility in random magnetic composites: application to two-dimensional and three-dimensional cases, Journal of molecular structure, 1160 (2018), 293.
Rollenhagen A., Lübke J.H.R. Cell Tissue Res 326 (2006) 221.
Harsan L. A., Poulet P. et al. J. Neuroscience Res., 85 (2007) 935.
Rushton W. J. Phisiol., 115 (1951) 101.
Reynolds P. J., Stanley H. E., Klein W. Phys. Rev. B, 21 (1980) 1223.
Poznansky O.P., Novikov V.V. J. Polym. Eng., 19 (1999) 223.
Hashin Z. Theory of composite materials (Pergamon Press, New York).
Latour L.L., Svoboda K., Mitra P.P., Sotak C.H. PNAS, 91 (1994) 1229.
Neeb H., Zilles K., Shah N. Neuroimage 31 (2006) 1156–1168.
Hall M., Alexander D. ISMRM 2008, 1779.
Posnansky O., Kupriyanova Y., Shah J. ESMRB 2008, 633.
Posnansky O. Effective diffusion in random composites measured by NMR, Journal of magnetism and magnetic materials 466 (2018) 92.
Posnansky O.P. Dynamic diffusion measured in binary composites measured by NMR, NANOINDUSTRY, 11, 7–8 (2018) 22.
Отзывы читателей
