Введение

Согласно статистике, острое нарушение мозгового кровообращения по причине развития гемодинамически значимого стеноза внутренней сонной артерии (ВСА) — одна из основных причин летальных исходов и длительной инвалидности [1, 2]. Классическая каротидная эндартерэктомия (КЭА) является одним из наиболее распространенных вариантов лечения этих больных [1, 2]. Однако результат применения заплаты во время этого вмешательства — локальная деформация геометрической формы сосуда, приводящая к изменению гемодинамики потока. Изучение гемодинамических характеристик кровотока на стенке сосуда крайне затруднено in vitro, а тем более in vivo. Поэтому вычислительные методы гидродинамики широко используются для оценки гемодинамических эффектов в сосудах персонально-специфической формы [3, 4].

Цель исследования — выявление зон наибольшего риска рестеноза в построенных компьютерных моделях при классической КЭА.

Материал и методы

Исходные данные для математического моделирования

Исходная геометрическая модель сосуда реконструирована на основе предоперационной компьютерной томографии пораженной левой каротидной бифуркации конкретного пациента. На рис. 1, а представлена пораженная часть сосуда (сонная артерия находится на переднем плане), сегмент кривой показывает предполагаемое положение стенки сосуда под атеросклеротической бляшкой. Пунктирной линией отмечена внутренняя поверхность бляшки в изображенной проекции.

Рис. 1. КТ-ангиограммы предоперационной области течения в сонных артериях.

Сосуд на переднем плане (а), контур сегментации в здоровой части ОСА (б), контур сегментации в пораженной части ВСА вблизи бифуркации (в), сегментации ОСА—ВСА и НСА (г), геометрическая форма базовой модели m0 с линией разреза (д).

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

Построение геометрических моделей

Реконструированная трехмерная модель здорового сосуда построена с использованием программного обеспечения SimVascular Updegrove (2016) и Salome Salome (2007). SimVascular использовали для построения сегментов сосудов (см. рис. 1, б—г). Эти сегменты с помощью скрипта Python нашей собственной разработки были затем импортированы в Salome для построения геометрической модели и сеток. Реконструированная трехмерная модель сосуда изображена на рис. 1, д. Далее эта модель называется базовой и обозначается через m0. Черная линия на стенке сосуда на рис. 1, д указывает линию разреза для последующей имитации имплантации заплаты.

Упомянутый скрипт также использовали для визуального моделирования результата операции КЭА. С его помощью можно визуально построить линию разреза и контуры заплаты на контурах сегментации сосуда (рис. 2). Информация о нарисованных линиях экспортируется скриптом в файл данных в виде списка значений ширины заплаты на ее пересечениях с контурами сегментации. Затем скрипт использует эти данные для изменения геометрической формы базовой модели, имитируя результат имплантации заплаты.

Рис. 2. Процесс визуального построения формы заплаты на контурах сегментации ОСА—ВСА.

Модифицированные модели, являющиеся результатами виртуальной КЭА, строили путем увеличения (или уменьшения) периметров всех сегментов модели m0, пересекающихся с линией разреза. Это позволяет имитировать любую форму заплаты или закрытие разреза без заплаты. Изменения периметров достигали путем масштабирования сегментов относительно их геометрических центров в соответствии со значениями, записанными в файле данных. Коэффициент растяжения по умолчанию полагается равным величине относительного приращения периметра, хранящейся в файле, деленной на 2π. Для сегментов, близких по форме к эллипсам с небольшим эксцентриситетом, такой выбор является достаточно точным. При желании для сегментов неправильной формы коэффициент масштабирования может быть откорректирован вручную. Вследствие масштабирования верхние (см. рис. 2) части проксимальных контуров ВСА смещались вверх и начинали пересекаться с проксимальными контурами НСА. Для корректировки этого выполняли небольшое параллельное смещение всех контуров НСА вдоль большей оси дистального контура ОСА. После этого строили геометрическую модель сосуда и расчетные сетки на ней. Сетки затем экспортировали в OpenFoam для последующих расчетов.

Для сравнительного анализа описанным выше способом были построены геометрические модели m1—m10. Модели m1—m9 имитируют результаты операции КЭА на модели m0 с имплантацией заплат p1—p9 соответственно (рис. 3). Модель m10 (не изображенная на рис. 3) имитирует закрытие разреза без имплантации заплаты.

Рис. 3. Формы заплат и сосудов после виртуальной КЭА.

Данные по формам заплат приведены в табл. 1, которая содержит значения ширины заплат в их поперечных сечениях вдоль линии разреза. Точка 0 колонки «Расстояние...» таблицы соответствует проксимальному концу линии разреза, точка 3,9 — дистальному. Линия разреза, изображенная на рис. 1, д одинакова для всех моделей. Как ее длина 3,9 см, так и расположение соответствуют фактическому разрезу, сделанному во время классической КЭА. Модель m10 построена путем уменьшения окружности просвета сосуда вдоль линии разреза. На это указывают отрицательные значения ширины виртуальной заплаты p10 в табл. 1.

Заметим, что в табл. 1 приведены приращения периметра сосуда после имплантации, ширина же самой заплаты до имплантации должна быть несколько большей.

Таблица 1. Значения ширины заплат p1—p10 вдоль линии разреза

Моделирование течения

Скорость течения U и давление p в построенных геометрических моделях описывали с помощью трехмерных нестационарных уравнений Навье—Стокса для вязкой несжимаемой жидкости:

(1)

с постоянной плотностью ρ=1050 кг/м³ и динамической вязкостью μ=3,675∙10^–3 Па∙с, где τ — тензор напряжений сдвига. Граничное условие непротекания было задано для U на боковой поверхности области течения, а условия параллельного течения — на входе и выходах. Границы области течения считались жесткими. Начальное значение скорости выбирали равным константе 0,15 м/с. В качестве граничных условий для p на входе и выходах устанавливали периодически меняющуюся разность давления, которая создавала периодическое течение, соответствующее постоперационным данным УДВ пациента:

— T=1,06 с — период кардиоцикла;

— Q=6,9 мл/с — объемный расход через ОСА;

— r=1,72 — отношение объемного расхода через ВСА к объемному расходу через НСА.

Метод построения граничных условий заключается в следующем. Сначала делали численный расчет, при котором на обоих выходах устанавливали нулевое давление, а на входе шаг за шагом с интервалом по времени 10 мс строили подходящую кривую давления. В ходе этого расчета приращения давления на каждом временном шаге выбирали вручную таким образом, чтобы результирующая кривая скорости на входе (рис. 4, а) соответствовала огибающей спектра УДВ в ОСА пациента. По результатам этого расчета вычисляли значение r, которое оказалось отличным от целевого значения 1,72. Далее для корректировки значения r на выходах вместо нулевого задавали переменное давление. А именно: на выходе НСА задавали такую же кривую давления, что и для ОСА, только уменьшенную по амплитуде с коэффициентом k=0,1 и с небольшим отставанием по фазе. На выходе ВСА задавали такую же кривую давления, что и на НСА, только инвертированную относительно оси абсцисс. Тем самым расход через НСА уменьшался, а через ВСА увеличивался, при этом кровоток через ОСА практически не менялся. После этого с помощью еще одной серии вспомогательных расчетов значение параметра k (и форму кривой давления, если необходимо) корректировали таким образом, чтобы величины Q и r стали близки к своим целевым значениям.

Рис. 4. Пример визуализации расчетов.

График скорости в ОСА (а). Линии тока при t=0,06 с (б), t=0,14 с (в), t=0,6 с (г). б—г — соответственно изображены линии тока в моменты времени 0,06 с, 0,14 с, 0,6 с, которые отмечены точками (а). Изменения цвета вдоль линии тока отражают скорость движения частиц крови в соответствии с предоставленной шкалой.

Численные расчеты и постпроцессинг

Численные расчеты проводили в OpenFoam методом конечных объемов с использованием алгоритма PISO, который наряду с коммерческим программным обеспечением, таким как Ansys Fluent, OpenFoam, является распространенным инструментом для выполнения гидродинамических расчетов и моделирования течения в кровеносных сосудах [4]. В описанных предварительных численных расчетах, выполненных для получения течения, соответствующего данным УДВ, использовали грубые расчетные сетки. После получения удовлетворительных результатов предварительных расчетов провели окончательные расчеты на мелких сетках. В результате расчетов получены динамические поля давления, скорости и градиента скорости в области течения для нескольких сердечных циклов с дискретизацией времени с шагом 10^–2 с. Информацию о параметрах течения и его производных характеристиках извлекали из результатов расчетов путем постпроцессинга, выполненного в ParaView Ayachit.

Для верификации этих результатов исследована их независимость от шага сетки. Установлено, что результаты существенно не изменяются при использовании сеток с числом узлов более 5·10⁵. Размер ячейки сетки на боковой поверхности установили равным половине размера ячейки внутри вычислительной области. Это сделано для повышения точности расчета гемодинамических показателей, которые выражаются через градиент скорости на стенке сосуда. Изучали также стабилизацию периодических колебаний. Как оказалось, процесс пульсовых колебаний можно считать установившимся, начиная со второго сердечного цикла [5]. В связи с этим далее предполагается, что время t=0 соответствует началу систолического нарастания скорости течения во втором кардиоцикле.

Пристеночное напряжение сдвига и гемодинамические показатели

Показатель WSS рассчитывали как тангенциальную составляющую t_w тензора сдвиговых напряжений τ на стенке сосуда. Показатели TAWSS, OSI и RRT рассчитывали по формулам (2) через усредненное значение t_w за 1 сердечный цикл [2, 4]:

, (2)

где T — длительность сердечного цикла, а |t_w| — евклидова норма вектора t_w. Для количественного сравнения показателей в некоторой зоне σ на стенке сосуда рассчитывали их средние значения по формулам:

(3)

где S — площадь зоны σ. Использовали также безразмерные среднее (I_RRT) и логарифмическое максимальное (M_RRT) значения показателя RRT для зоны σ, которые вычисляли по формулам:

, (4)

где TAWSS_CCA — усредненное значение TAWSS в цилиндрической части ОСА на расстоянии трех радиусов ОСА от точки разветвления сосуда [24].

Для оценки патологических значений показателя RRT использовали неравенства TAWSS<0,4 Па и OSI>0,3, приведенные в исследовании G. Harrison и соавт. (2014) [6]. В этом случае формула (2) дает соответствующее критическое значение RRT=6,25 Па^–1. Далее значения RRT, превышающие критическое, будут считать патологическими.

Результаты

Результаты по течению в базовой модели

Проведен анализ результатов расчетов поля скоростей и гемодинамических показателей для базовой модели m0. На рис. 4 изображены некоторые визуализации для следующих значений параметров: T=1,06 с, Q=6,9 мл/с, r=1,72. На рис. 4, а приведен график модуля вектора скорости в центре проксимального поперечного сечения ОСА.

Изучали также влияние частоты и амплитуды пульса на распределение гемодинамических показателей на стенке сосуда в зоне его атеросклеротического поражения. В табл. 2 приведены значения T, Q и r для 4 вариантов расчета (a—d).

Таблица 2. Параметры течения T, Q, r в расчетах a—d.

На рис. 5 изображены линии уровня RRT на выбранном участке стенки сосуда, соответствующие одноименным наборам параметров из табл. 2. Линии уровня 0 соответствуют критическому значению RRT=6,25 Па^–1. Линии уровня 1—3 соответствуют значениям RRT в 2 раза, 4 раза и 8 раз большим, чем критическое значение. Цветом отображается значение ln(RRT+1) между 0 и 8 в соответствии с прилагаемой шкалой. На рис. 6 показано расположение выбранного участка на стенке сосуда. Рис. 5 демонстрирует, что сколько-нибудь существенного качественного изменения в распределении показателя RRT при вариации параметров течения, приведенных в табл. 2, не происходит.

Рис. 5. Линии уровня RRT (в Па^–1) и распределение ln(RRT+1) для наборов параметров a—d в табл. 2.

Рис. 6. Топологическая структура поля TAWSS в окрестности его точек покоя (точки 1—4).

На рис. 7 изображена в несколько увеличенном виде зона локальных максимумов индекса RRT из рис. 5, а в сочетании с фазовым портретом векторного поля TAWSS.

Рис. 7. Зоны z_1 и z_2 и стационарные точки 1—4 поля TAWSS для базовой модели.

Стрелки на рисунке указывают локальное направление поля TAWSS, цвет отражает значение ln(RRT+1) в диапазоне [0, 8]. Черные точки — это неподвижные точки поля TAWSS, жирные линии — это сепаратрисы двух седловых точек 2 и 4, тонкие линии представляют собой регулярные траектории.

Три зоны σ1—σ3 (см. рис. 7) с равными площадями S были выбраны для сравнительной количественной оценки интегральных значений RRT. Центрами этих зон являются стационарные точки 1—3 поля TAWSS. Средние значения индексов, рассчитанных для каждой зоны по формулам (3), приведены в табл. 3.

Таблица 3. Средние значения индексов в зонах σ1—σ3

Гемодинамические показатели для моделей m0—m10

Результаты численных расчетов с параметрами течения Q и r, отлич ающимися не более чем на 0,7% от результата расчета «а» для базовой модели (см. табл. 2), получены для моделей m1—m10, Изменения геометрической формы сосуда, вызванные имплантацией виртуальной заплаты, приводят к изменению (в пределах 10%) расчетных значений Q и r при неизменных граничных условиях. В силу этого для каждой из моделей m1—m10 проводили вспомогательные расчеты для корректировки коэффициента k, чтобы уложиться в интервал [–0,7%, 0,7%] изменения параметров Q и r.

Для получения сравнительных количественных характеристик гемодинамических показателей в каждой модели выбрали области, где одновременно удовлетворялись оба критерия: TAWSS <0,4 Па и OSI >0,3. Эти области для всех моделей состоят из двух зон — z_1 и z_2, первая из которых содержит стационарные точки 1, 2, вторая — стационарные точки 3, 4 поля TAWSS. На рис. 7 изображены эти зоны для базовой модели. Значения показателей I_RRT и M_RRT в зонах z_1 и z_2, рассчитанные по формулам (4) для каждой модели, приведены в табл. 4. В строках таблицы приведены значения I_RRT и M_RRT для моделей m0—m10 в процентах к соответствующим значениям для базовой модели m0.

Обсуждение

Построение геометрических моделей

Наша задача состояла не только в реконструкции сегментации исходного здорового сосуда по КТА данным пораженного сосуда. На неповрежденных участках сосуда (см. рис. 1, б) границы его просвета совпадают с его внутренней стенкой, и сегментацию в SimVascular можно выполнить в автоматическом режиме. Однако в зонах атеросклеротического поражения (см. рис. 1, в) такого совпадения нет, поэтому каждый контур строили с использованием ручной корректировки с последующим Фурье-сглаживанием. Например, на рис. 1, в правую часть контура строили по внешней границе рентгеноконтрастной кальцинированной бляшки. Кроме того, в силу ограниченной разрешающей способности оборудования граница просвета сосуда размыта, следовательно, при ее распознавании присутствует некоторая неопределенность масштабирования. Поэтому площади некоторых распознанных поперечных сечений сравнивали с площадями тех же поперечных сечений, полученных средствами УДВ, вычисляли поправочный коэффициент, а затем все контуры масштабировали в соответствии с этим коэффициентом. В итоге получены сегментации ОСА—ВСА и НСА (см. рис. 1, г), по которым строили базовую модель m0 (см. рис. 1, д).

Моделирование течения

Моделирование течения выполняли в предположении, что кровь — это ньютоновская жидкость. Давно установлено, что неньютоновская реология проявляется в основном в сосудах малого диаметра и капиллярах. Поэтому для моделирования течения в сонной артерии, как правило, используют ньютоновскую реологию [7].

Предположение о жесткой стенке сосуда принято в данной работе по ряду причин. Во-первых, здоровая сонная артерия располагается в ложе из эластичных и податливых структур из соединительной и жировой ткани, внутри которых и происходит перемещение ее стенки. С появлением кальцинированной атеросклеротической бляшки пульсовые движения стенки сосуда в ее зоне прекращаются из-за отсутствия гибкости. Удаление же бляшки с последующей имплантацией заплаты тем не менее не восстанавливает гибкость стенки сосуда в полной мере, так как податливые окружающие структуры вследствие операционной травмы замещаются ригидной соединительной тканью. Заметим, что отсутствие перемещений сосудистой стенки в ответ на изменение давления внутри сосуда с медицинской точки зрения — самый неблагоприятный гемодинамический вариант, который и был выбран в работе в качестве модели с наихудшим прогнозом. Во-вторых, исследуемые участки сосуда значительно короче длины пульсовой волны давления, поэтому в случае упругой стенки можно предположить, что они расширяются и сужаются почти синхронно. Поэтому если оценивать амплитуду пульсаций диаметра сосуда 5%, то амплитуда колебаний пикового значения WSS не будет превышать 15% [8]. Однако для интегральных показателей (2) это пиковое систолическое значение усредняется по всему периоду кардиоцикла, и его влияние на показатели будет кратно меньше. Наконец, введение в модель дополнительных исходных данных, необходимых для расчетов с деформируемой стенкой, таких как ее механические свойства, толщина, требует их точной оценки. В противном случае предполагаемые погрешности в этих данных могут внести неконтролируемую неопределенность в результаты расчетов [8]. По этим причинам предположение о жесткости стенки довольно распространено для сонной артерии.

Анализ результатов для модели m0

Известно, что течение в каротидной бифуркации имеет сложную структуру, оно всегда содержит изменяющиеся во времени зоны рециркуляции и застоя [9]. В нашем случае качественное поведение линий тока, изображенное на рис. 4, б—г, соответствует опубликованным результатам.

Для определения вероятного положения зон инициации атеросклеротических бляшек и их зависимости от параметров кровотока получены распределения гемодинамических показателей (рис. 5) для четырех вариантов расчета (см. табл. 2). На рис. 5 видно, что расположение зон с высокими значениями RRT практически не зависит от изменения параметров T и Q.

На рис. 6 изображен фазовый портрет векторного поля TAWSS, соответствующий набору параметров a из табл. 2. В рассматриваемой области поле TAWSS имеет четыре стационарные точки, отмеченные на рис. 6 номерами 1—4. Точки 1 и 3 являются устойчивыми узлами, точки 2 и 4 — седловыми точками (неустойчивыми). Стрелки на рис. 6 соответствуют направлению действия напряжения сдвига на клетки эндотелия. Фазовые траектории поля TAWSS соответствуют путям переноса компонентов крови и химических веществ вдоль стенки сосуда. Расчеты, проведенные также для наборов параметров b—d из табл. 2, показали, что топологическая структура векторного поля TAWSS в рассматриваемой области идентична для всех четырех вариантов a—d. Этот факт вместе с результатами, приведенными на рис. 5, позволяет сделать вывод, что изменение параметров течения меньше влияет на конфигурацию зон риска по сравнению с изменениями геометрической формы сосуда.

Из табл. 3 видно, что с точки зрения ухудшения значений гемодинамических показателей зоны σ1—σ3 расположены в порядке возрастания их нумерации. Заметим, что расположение реальной атеросклеротической бляшки на рис. 1, а коррелирует с расположением рисковых зон σ1—σ3 на рис. 6. Возникновение бляшки меняет локальную гемодинамику, что приводит к распространению патологических областей и дальнейшему ее росту. Таким образом, можно предположить, что формирование бляшки в здоровом сосуде m0 начиналось с особых точек поля TAWSS. Наибольшая вероятность возникновения патологических явлений, согласно табл. 3, при прочих равных условиях относится к окрестности σ₃.

Таблица 4. Показатели I_RRT и M_RRT для моделей m0—m10, %

Сравнительный анализ результатов для моделей m0—m10

Сравнительные результаты расчетов показателей I_RRT и M_RRT в зонах z_1 и z_2 для всех моделей приведены в табл. 4. Зоны z_1 и z_2 — это зоны наиболее вероятного возникновения атеросклеротических явлений, большие значения показателей в них означают больший риск рестеноза. В соответствии с табл. 1 заплаты p1—p4 имеют приблизительно постоянную ширину по всей длине, за исключением концов. Заплаты такой формы практикующие хирурги применяют наиболее часто, при этом вопрос о выборе ширины заплаты остается за ними. В мировой литературе не обнаружено каких-либо правил выбора ширины заплаты для какого-либо случая. Однако обсуждается также вопрос о выборе между имплантацией заплаты и закрытием разреза без имплантации заплаты [10]. Данные табл. 4 позволяют сделать вывод, что с точки зрения риска возникновения рестеноза в зонах z_1 и z_2 среди заплат p1—p4 предпочтительнее использовать p3. Предпочтительность использования заплаты перед непосредственным закрытием разреза в нашем случае подтверждается данными табл. 4 для модели m10, имитирующей такой вариант операции.

В настоящей работе в отличие от упомянутых выше исследовали также заплаты неравномерной ширины. Как видно из табл. 1, заплаты p5 и p6 были получены из заплаты p3 путем сужения и расширения ее соответственно в проксимальной части. Оба варианта приводят к ухудшению значений показателей в табл. 4.

Заплаты p7 и p8 получены из p2. Заплата p7 — это p2, суженная в своей проксимальной части. Заплата p8 представляет собой укороченную на 0,7 см p2, имплантированную в разрез, с непосредственным закрытием разреза в проксимальной части с уменьшением окружности просвета сосуда на 0,1 см. В обоих случаях наблюдается небольшое улучшение показателей в зоне z_1 и ухудшение их в зоне z_2 по сравнению с заплатой p2. Заплата p9 получена из p4 сужением ее в дистальной части, и это не приводит к сколько-нибудь существенным изменениям показателей в табл. 4.

Следовательно, лучшие результаты демонстрируют заплаты средней ширины, без сужений и расширений на концах. Наилучшим вариантом является выбор заплаты p3, а наихудшим — закрытие разреза без заплаты.

Заключение

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

Авторы заявляют об отсутствии конфликта интересов.