Журнал вычислительной математики и математической физики, 2021, T. 61, № 9, стр. 1571-1584

1. ВВЕДЕНИЕ

Церебральная артериовенозная мальформация (АВМ) – это врожденная патология развития сосудов головного мозга, при которой осуществляется прямой сброс крови из артериального бассейна в венозный, минуя капиллярную сеть сосудов. Наличие АВМ нарушает нормальное кровоснабжение мозга, изменяя скорость и давление крови в близлежащей части сосудистой системы. Из-за отсутствия капиллярной сети уменьшается сопротивление соответствующей части кровеносной системы, возникает ишемия близлежащего мозгового вещества. Чаще всего АВМ становятся симптоматическими в возрасте от 20 до 50 лет и выявляются в этот период жизни, причем внутричерепное кровоизлияние является наиболее распространенным клиническим проявлением АВМ, по данным [1] его частота варьируется от 30 до 82%. Высокий риск смертности и инвалидизации из-за внутричерепного кровоизлияния определяет опасность этой патологии [2].

На данный момент эмболизация – это наиболее прогрессивный метод лечения АВМ, представляющий собой избирательное выключение кровеносных сосудов из кровотока путем их заполнения специальной эмболизирующей композицией (эмболизатом) в процессе эндоваскулярной нейрохирургической операции. Эндоваскулярное лечение артериовенозных мальформаций достигается микрокатетерной доставкой агентов, таких как N-бутил-2-цианоакрилат или неадгезивный сополимер этилена и винилового спирта [3], [4]. По данным исследования 408 пациентов с АВМ, которые лечились эндоваскулярно, Baharvahdat и соавт. отмечали 11%-ную частоту геморрагических осложнений, связанных с лечением [5]. По данным исследования группы из 192 пациентов с тотально эмболизированной артериовенозной мальформацией летальность за период госпитализации составила 9.3% [6]. Таким образом, несмотря на хорошо развитую технику операций по эмболизации, риск интраоперационного разрыва сосудов по-прежнему вызывает озабоченность. Гемодинамика, связанная с внутричерепными АВМ, сложна и изменяется по мере изменения ее морфологии и ангиоархитектуры. В связи с вышесказанным, математическое моделирование процесса эмболизации АВМ является актуальной задачей.

Для адекватного моделирования АВМ необходимо знать ее геометрическую структуру. Современные методы реконструкции геометрии АВМ по данным нейровизуализации (КТ, МРТ, церебральная ангиография) позволяют выделить сосуды со средним диаметром не менее $0.5$ мм, чего недостаточно для определения in vivo детальной геометрической структуры АВМ, состоящей из очень большого количества сросшихся переплетенных тонких сосудов, диаметр которых может достигать 0.1 мм. Также измерение in vivo гемодинамических параметров представляется возможным только в достаточно крупных примыкающих к АВМ сосудах, а измерения непосредственно внутри АВМ серьезно затруднены и зачастую невозможны. Из-за отсутствия необходимых данных существующие математические модели имеют упрощенный характер, хотя и позволяют получить качественно правильное описание гемодинамики рассматриваемой патологии.

Взаимодействие АВМ и кровотока в окружающих сосудах часто изучается на основе аналогии с электрическими или гидравлическими сетями, например, в работах [7]–[9]. Такие модели позволяют оценить влияние различных сценариев эмболизации на перестройку кровотока и согласуются с общим медицинским взглядом на гемодинамику АВМ. Для описания совместного течения крови и эмболизата используются различные подходы, например в работе [10] использовали двухфазную модель потока для моделирования движения вязкой капли эмболизата через точку бифуркации или мультифуркации на основе безразмерных уравнений Навье–Стокса для несжимаемых жидкостей. В работе [11] описывается модель эмболизации АВМ на основе концепции двухжидкостного моделирования и скалярного транспорта, а взаимодействие эмболизата с кровью и его затвердевание здесь имитировались увеличением вязкости одной из фаз.

В настоящей работе для моделирования совместного течения крови и эмболизирующей композиции внутри АВМ используется одномерная модель двухфазной фильтрации, построенная на основе клинических данных реальных пациентов, полученных во время нейрохирургических операций в НМИЦ им. акад. Е.Н. Мешалкина. Такой подход, по нашему мнению, оправдан для описания эмболизации мелкососудистых компартментов АВМ. Поскольку АВМ является частью кровеносной системы, то расход крови на входе в АВМ меняется с течением операции, за счет перераспределения крови в здоровые сосуды. Этот эффект учитывается сопряжением фильтрационной модели с моделью здорового сосуда постоянного сопротивления. Математическое описание данного процесса сводится к решению начально-краевой задачи для интегродифференциального уравнения с невыпуклым потоком. Для численных расчетов разработана монотонная модификация схемы CABARET, которая с высокой точностью локализует сильные и слабые разрывы, возникающие в решении данной задачи. Основная цель работы заключается в отыскании оптимального с точки зрения безопасности и эффективности сценария эмболизации АВМ. Целевой функционал и ограничения, возникающие в такой задаче оптимального управления, выбираются в соответствии с медицинскими показаниями. Управлением является зависящая от времени функция, определяющая объемный расход эмболизата на входе в АВМ. Проведено сравнение оптимальных сценариев эмболизации АВМ при различных режимах перетока крови в здоровые сосуды для конкретных пациентов.

2. ПОСТАНОВКА ЗАДАЧИ

Движение крови и эмболизата через мелкососудистую рацемозную АВМ моделируется как процесс двухфазной фильтрации несмешивающихся несжимаемых жидкостей, где вытесняемой фазой является кровь, а вытесняющей – эмболизат. Мы будем предполагать однородность структуры АВМ и рассматривать пространственно-одномерный случай, для которого соответствующее математическое описание, основанное на законе Дарси и законах сохранения масс для каждой фазы, впервые было предложено в работе С. Баклея, М. Леверетта [12]. Пусть $S(t,x) \in [0,\;1]$ – локальная насыщенность АВМ кровью (концентрация крови) во время эмболизации. С учетом этого концентрация эмболизата будет равна $1 - S(t,x)$. Уравнение для $S(t,x)$ внутри АВМ можно записать по аналогии с уравнением Баклея–Леверетта [13]:

где $x \in [0,L]$, $L$ – длина АВМ (которую будет определять по клиническим данным); $m$ – пористость АВМ (которую будем считать постоянной и равной единице); $A$ – площадь поперечного сечения АВМ (которую будем предполагать постоянной и определять по клиническим данным); $f(S)$ – функция Баклея–Леверетта, которая задает закон распределения потоков фаз внутри АВМ и удовлетворяет условиям: (конкретный вид этой функции, получаемый по данным интраоперационного мониторинга, будет приведен ниже); есть расход смеси двух фаз внутри АВМ, где ${{Q}_{b}} = Qf(S)$ – расход крови, а ${{Q}_{e}} = Q(1 - f(S))$ – расход эмболизата, граничное значение которого ${{Q}_{e}}(t,0)$ определяется сценарием проведения операции.

Расход смеси $Q(t)$ не зависит от пространственной координаты $x$, поскольку кровь и эмболизат рассматриваются как несжимаемые жидкости. Здесь и везде далее под расходом жидкости понимается ее объемный расход, представляющий собой массовый расход жидкости, деленный на ее плотность. С учетом этого из соотношения (2.3) следует, что плотность потока смеси крови и эмболизата внутри АВМ задается формулой

где ${{\rho }_{b}}$ – плотность крови и ${{\rho }_{e}}$ – плотность эмболизата.

Рассмотрим систему кровеносных сосудов (фиг. 1), которая в упрощенном виде воспроизводит перераспределение артериального потока крови вблизи АВМ, изображенной на фиг. 1 серой фигурой. Кривыми ${{l}_{{1j}}}$, соединяющими на фиг. 1 узлы ${{B}_{1}}$ и ${{B}_{j}}$, где $j = 2,\;3$, моделируются здоровые мозговые сосуды, через которые в течение всей операции протекает только кровь. ${{B}_{1}} - {{B}_{2}}$ – здоровый кровеносный сосуд постоянного сопротивления ${{R}_{{12}}}$; ${{B}_{1}} - {{B}_{3}}$ – сосуд, подводящий кровь к АВМ (афферент) постоянного сопротивления ${{R}_{{13}}}$; ${{B}_{3}} - {{B}_{4}}$ – участок кровеносной системы, включающий в себя АВМ (серая фигура) с расходом крови ${{Q}_{b}}$ и расходом эмболизата ${{Q}_{e}}$, расход эмболизата на входе в АВМ обозначен ${{q}_{e}}$; ${{Q}_{1}}$ – расход артериальной крови, приходящей в узел ${{B}_{1}}$, ${{Q}_{{12}}}$ и ${{Q}_{{13}}}$ – потоки крови в соответствующих сосудах. Будем предполагать, что геометрия этих сосудов во время операции не изменяется, в силу чего расходы крови ${{Q}_{{1j}}}(t)$ в этих сосудах не зависят от пространственных координат и изменяются в процессе эмболизации только как функции времени. Поскольку за счет механизма ауторегуляции церебральный поток крови изменяется достаточно мало, будем считать, что входной расход артериальной крови ${{Q}_{1}}$, приходящей в узел ${{B}_{1}}$, является постоянным. Обозначим через ${{p}_{i}}$ давление крови в узле ${{B}_{i}}$ на фиг. 1, где $i = 1,\;2,\;3$. Давления ${{p}_{1}}(t)$ и ${{p}_{3}}(t)$ изменяются в течение эмболизации. Начальные значения давлений ${{p}_{1}}(0)$, ${{p}_{3}}(0)$ задаются в соответствии с клиническими данными. Давление ${{p}_{2}}$ на выходе из здорового сосуда ${{l}_{{12}}}$ будем считать постоянным, в то время как давление ${{p}_{4}}$ на выходе из АВМ с течением времени изменяется. Конкретные значения давлений ${{p}_{2}}$, ${{p}_{4}}$ считаются известными и получаются из клинических данных, причем ${{p}_{4}}$ зависит от распределения эмболизата в теле АВМ и задается как функция от средней концентрации крови в АВМ ${{p}_{4}} = {{p}_{4}}(\bar {S})$, где $\bar {S}(t) = 1{\text{/}}L\int_0^L {S(x,t)} dx$.

Из закона гидравлического сопротивления [14] для здоровых кровеносных сосудов получаем

где ${{r}_{{1j}}}$ – локальное (удельное) гидравлическое сопротивление сосуда, соединяющего узлы ${{B}_{1}}$ и ${{B}_{j}}$; $d{{l}_{{1j}}}$ – дифференциал длины дуги кривой ${{l}_{{1j}}}$, которой моделируется этот сосуд. Из закона Дарси для АВМ [13] в общем случае находим: где есть локальное сопротивление прохождению крови внутри АВМ, в котором ${{\eta }_{b}}$ – вязкость крови, $K(x)$ – абсолютная проницаемость АВМ, а ${{k}_{b}}(x,S)$ – относительная проницаемость для крови (способ их определения по клиническим данным будет описан ниже).

Исключая из уравнений (2.4) давление ${{p}_{1}}$ с учетом того, что

получаем Исключая из уравнений (2.5) и (2.6) давление ${{p}_{3}}$, приходим к формуле которую с учетом того, что ${{Q}_{{13}}} = {{Q}_{b}}(t,0)$ и ${{Q}_{b}}(t,x) = Q(t)f(S(t,x))$, можно переписать в виде где $a(\bar {S}(t)) = {{p}_{2}} - {{p}_{4}}(\bar {S}(t)) + {{R}_{{12}}}{{Q}_{1}}$. При известных значениях расхода ${{Q}_{1}}$, давлений ${{p}_{2}}$, ${{p}_{4}}$ и сопротивлений ${{R}_{{12}}}$, ${{R}_{{13}}}$, а также при заданных функциях $f(S)$ и ${{r}_{b}}(x,S)$, уравнения (2.1) и (2.7) образуют замкнутую систему интегродифференциальных уравнений для определения расхода смеси $Q(t)$ и концентрации крови $S(t,x)$ внутри АВМ.

Поскольку в начальный момент времени $t = 0$ внутри АВМ находится только кровь, то начальное значение ее концентрации равно единице, т.е.

С учетом условий (2.2) и (2.8) из уравнения (2.7) получаем начальное значение расхода крови (в смеси отсутствует эмболизат)

Из физической постановки задачи следует, что в течение всего времени $T$ проведения хирургической операции выполнено неравенство

в силу чего характеристики дифференциального уравнения (2.1) распространяются со скоростями $Q(t)f{\kern 1pt} {\text{'}}(S(t,x)) > 0$. Это означает, что для корректной постановки начально-краевой задачи для уравнения (2.1) на отрезке $[0,L]$ граничное условие необходимо задавать на левой границе этого отрезка. Поскольку процесс лечения происходит путем заполнения АВМ эмболизатом через подводящую артерию (точка ${{B}_{3}}$ на фиг. 1), то расход вводимого эмболизата с учетом формулы (2.7) задает интегральное граничное условие для $S(t,0)$: в котором граничное значение концентрации крови $S(t,0)$ зависит от значений этой концентрации на всем отрезке $[0,L]$.

Следует заметить, что в рассматриваемой задаче $a(\bar {S}) > 0$ в силу того, что из физиологии перепад давления ${{p}_{2}} - {{p}_{4}}(\bar {S})$ между венозным концом здорового сосуда ${{B}_{2}}$ и венозным концом мальформации ${{B}_{4}}$ по абсолютной величине не превосходит перепада давления ${{R}_{{12}}}{{Q}_{1}}$ между артериальным ${{B}_{1}}$ и венозным концом здорового сосуда ${{B}_{2}}$, который реализуется при полном прекращении потока крови через АВМ. Поэтому знаменатель в формуле (2.10) не может быть отрицательным, и поскольку по определению $f(S) \geqslant 0$, то из формулы (2.10) следует ограничение на возможный расход эмболизата ${{q}_{e}}(t)$:

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

3. РАЗНОСТНАЯ СХЕМА

В работе [15] для численного решения гиперболических уравнений была предложена трехслойная по времени и двухточечная по пространству схема Upwind Leapfrog, которая имеет второй порядок аппроксимации на гладких решениях, является явной и условно устойчивой при числах Куранта $r \in \left( {0,\;1} \right]$. Детальный анализ этой схемы был проведен в работах [16], [17], в которых с учетом кососимметричности своего пространственного шаблона она была названа схемой CABARET. Основные достоинства этой схемы связаны с тем, что она задана на компактном пространственном шаблоне, является обратимой по времени и точной при двух различных числах Куранта $r = 0.5,\;1$, что наделяет ее уникальными диссипативными и дисперсионными свойствами [17]. Для численного решения уравнений одномерной газовой динамики был разработан балансно-характеристический вариант схемы CABARET [18], который с учетом коррекции потоковых переменных (необходимой для монотонизации разностного решения на ударных волнах) показал высокую точность при расчете классического теста Blast Wave [19]. В настоящее время для численного моделирования пространственно многомерных газодинамических [20] и гидравлических [21] течений широко применяется двухслойная по времени форма записи схемы CABARET [22], [23]. Монотонность этой схемы при аппроксимации скалярного закона сохранения с выпуклым потоком исследовалась в [24] и с невыпуклым потоком – в [25]. В настоящей работе численный расчет начально-краевой задачи (2.8), (2.10) для системы интегродифференциальных уравнений (2.1), (2.7) проводится на основе предложенной в [25], [26] монотонной модификации схемы CABARET. Краткое описание численного алгоритма приводится ниже.

Разобьем отрезок $[0,L]$, на котором решается задача, на $N$ равных частей длины $h = L{\text{/}}N$ и зададим равномерную по пространству прямоугольную разностную сетку

в которой $h$ – постоянный шаг сетки по пространству, а ${{\tau }_{n}}$ — шаг сетки по времени, определяемый из условия устойчивости Куранта. В схеме CABARET используются потоковые $u_{j}^{n} = U({{x}_{j}},{{t}_{n}})$ и консервативные $U_{{j + 1/2}}^{n} = U({{x}_{{j + 1/2}}},{{t}_{n}})$ переменные, заданные соответственно в целых ${{x}_{j}}$ и полуцелых ${{x}_{{j + 1/2}}} = {{x}_{j}} + h{\text{/}}2$ пространственных узлах разностной сетки (3.1), где $U$ – приближенное численное значение концентрации крови $S$.

Пусть $u_{j}^{n}$, $U_{{j + 1/2}}^{n}$ – известное численное решение задачи (2.1), (2.7), (2.8), (2.10) на $n$-м временном слое ${{t}_{n}}$; при $n = 0$ – точная сеточная аппроксимация начального условия (2.8):

Численное решение $u_{j}^{{n + 1}}$, $U_{{j + 1/2}}^{{n + 1}}$ на $(n + 1)$-м временном слое ${{t}_{{n + 1}}}$ будем находить по схеме CABARET в несколько этапов. На первом этапе по формуле где $r = 0.5$, вычисляется шаг схемы по времени на $n$-м временном слое.

На втором этапе путем аппроксимации со вторым порядком интегрального уравнения (2.7) находится численное значение потока смеси на $n$-м временном слое

после чего на третьем этапе по разностным уравнениям вычисляются значения консервативных переменных $U_{{j + 1/2}}^{{n + 1/2}} = U({{x}_{{j + 1/2}}},{{t}_{{n + 1/2}}})$ на полуцелом временном слое ${{t}_{{n + 1/2}}} = {{t}_{n}} + {{\tau }_{n}}{\text{/}}2$. При записи формул (3.2) и (3.3) использованы сокращенные обозначения $f_{j}^{n} = f(u_{j}^{n})$ и $f_{{j + 1/2}}^{n} = f(U_{{j + 1/2}}^{n})$.

На четвертом этапе находятся значения численных потоков $f_{j}^{{n + 1/2}} = f(u_{j}^{{n + 1/2}})$, в которых $u_{j}^{{n + 1/2}} = U({{x}_{j}},{{t}_{{n + 1/2}}})$. В результате аппроксимации со вторым порядком интегрального граничного условия (2.10) находится значение численного потока на левой границе расчетной области

где $f_{{j + 1/2}}^{{n + 1/2}} = f(U_{{j + 1/2}}^{{n + 1/2}})$. Численные потоки $f_{j}^{{n + 1/2}}$ при $j = \overline {1,N} $ вычисляются также, как в работе [25].

На пятом этапе по формуле (3.2), в которой индекс $n$ заменен на $n + 1{\text{/}}2$, определяется численное значение потока смеси ${{Q}_{{n + 1/2}}}$ на полуцелом временном слое, после чего на шестом этапе по разностным уравнениям

вычисляются значения консервативных переменных $U_{{j + 1/2}}^{{n + 1}}$ на $(n + 1)$-м временном слое. На седьмом этапе на $(n + 1)$-м временном слое находятся значения потоковых переменных $u_{j}^{{n + 1}}$, которые вычисляются также, как в работе [25].

4. ЗАДАЧА ОПТИМАЛЬНОЙ ЭМБОЛИЗАЦИИ

Для постановки задачи оптимальной эмболизации необходимо определить критерии эффективности операции и необходимые медицинские ограничения. Эффективность эмболизации мы будем определять тем, насколько полно заполняется АВМ с помощью эмболизата. Математически это условие требует нахождения оптимального расхода эмболизата ${{q}_{e}}(t)$, при котором достигается минимум функционала

где ${{S}_{{{{q}_{e}}}}}(t,x)$ – решение начально-краевой задачи (2.1), (2.7), (2.8), (2.10) с заданным расходом эмболизата ${{q}_{e}}(t)$. Медицинские ограничения требуют выполнения следующих двух условий: Первое из этих условий означает, что давление ${{p}_{3}}$ на артериальном конце АВМ не превышает заданного критического давления ${{p}_{{{\text{max}}}}}$, а второе условие позволяет избежать проникновения эмболизата в венозное русло. Нарушение этих условий увеличивает риск мозгового кровоизлияния.

Рассмотрим решение начально-краевой задачи (2.1), (2.7), (2.8), (2.10), удовлетворяющее ограничениям (2.11), (4.2), (4.3). Одна из возможных методик проведения операции заключается в задании некоторого начального расхода эмболизата, а затем снижение этого расхода до нуля при завершении операции. В соответствии с этим рассмотрим сценарий подачи эмболизата специального вида ${{q}_{e}}(t) = \gamma Q(0)E(t)$, где $\gamma $ – безразмерный параметр, $E(t)$ – безразмерная функция, которая имеет следующий вид:

где ${{\tau }_{1}} > 0$, ${{\tau }_{2}} > 0$, ${{\tau }_{1}} + {{\tau }_{2}} = T$. Из формулы (4.4), в частности, следует, что если $\gamma = 1$, то в начальный момент времени расход эмболизата совпадает с расходом крови через АВМ до эмболизации, т.е. ${{q}_{e}}(0) = Q(0)$, а затем в течение операции расход эмболизата уменьшается по линейному закону и становится равным нулю к концу операции. В сценарии эмболизации, который задается формулой (4.4), управление подачей эмболизата на вход АВМ происходит при помощи трех параметров: $\gamma $, ${{\tau }_{1}}$ и ${{\tau }_{2}}$. Мы расширим класс допустимых управлений, включая в рассмотрение разрывную функцию получаемую из непрерывной функции (4.4) при ${{\tau }_{2}} \to 0$.

Далее решение задачи оптимальной эмболизации ищется в классе функций (4.4), (4.5) при различных значениях входных параметров $\gamma $, ${{\tau }_{1}}$ и ${{\tau }_{2}}$.

5. ПОСТРОЕНИЕ ФУНКЦИИ БАКЛЕЯ–ЛЕВЕРЕТТА ПО КЛИНИЧЕСКИМ ДАННЫМ

Для решения задачи эмболизации используются функции Баклея–Леверетта, построенные по клиническим данным мониторинга гемодинамических параметров во время нейрохирургических операций в НМИЦ им. Мешалкина [27]. Для этого при помощи аппарата Philips ComboMap с датчиком Philips ComboWire (диаметр датчика – 0.36 мм, длина датчика – 1.85 м) в процессе нейрохирургических операций измеряются скорость и давление кровотока внутри различных сосудов головного мозга вблизи патологии. Таким способом получены значения давления и скорости на входе АВМ (в артерии) перед, во время и после проведения операции по эмболизации, а также на выходе из АВМ (в вене) перед и после проведения операции.

Функция Баклея–Леверетта $f(S)$ полностью определяется относительными фазовыми проницаемостями и вязкостями крови и эмболизата:

здесь ${{k}_{b}}(S)$, ${{\eta }_{b}}$ – относительная фазовая проницаемость и коэффициент динамической вязкости крови, ${{k}_{e}}(S)$, ${{\eta }_{e}}$ – эмболизирующего вещества. Далее используются известные вязкости крови [28] и распространенного эмболизата ONYX18 [4], которые составляют: ${{\eta }_{b}} = 4$ сП, ${{\eta }_{e}} = 18$ сП.

Будем считать, что зависимость между падением давления и объемным расходом крови в теле АВМ приближенно описывается законом Дарси:

где $K$ – постоянная абсолютная проницаемость пористой среды, которая вычисляется по данным мониторинга до начала эмболизации при $t = 0$, когда через АВМ течет только кровь, в силу чего $\bar {S}(0) = 1$, ${{k}_{b}}(\bar {S}(0)) = 1$. С учетом этого из формулы (5.2) получаем Из закона Дарси (5.2) с учетом (5.3) можно получить следующую формулу для относительной проницаемости крови: Из формулы (5.4) видно, что полученная таким образом относительная проницаемость крови не зависит от вязкостей крови и эмболизата, а также от геометрических свойств АВМ.

Для аналитической аппроксимации относительной фазовой проницаемости крови используется модель Кори [29], в которой ${{k}_{b}}(S) = {{S}^{\alpha }}$, где параметр $\alpha > 1$ при известной функции $\bar {S}(t)$ с учетом (5.4) определяется из формулы

Относительная проницаемость эмболизата ${{k}_{e}}(S)$ предполагается симметричной относительно проницаемости крови в следующем смысле: Таким образом, функция Баклея–Леверетта принимает следующий вид:

Для нахождения из формулы (5.5) параметра $\alpha $ методом наименьших квадратов используются данные операций, в которых достигалась тотальная эмболизация АВМ, означающая, что поток крови через АВМ в конце операции прекратился. По данным $i$-го измерения в момент времени ${{t}_{i}}$ восстанавливаются объемный расхода крови ${{Q}_{b}}({{t}_{i}})$, перепад давления $\Delta p({{t}_{i}})$ между артериальным и венозным концами патологии, а также средняя по объему концентрация крови $\bar {S}({{t}_{i}})$ в теле АВМ. Давление ${{p}_{4}}$ на венозном конце АВМ во время операции определяется путем линейной интерполяции значений этого давления до и после операции. На фиг. 2 кружками и треугольниками показаны значения функции Баклея–Леверетта для пациентов $1$ и $2$ соответственно, полученные из клинических данных по формуле (5.1), в которой относительная проницаемость крови определяется по формуле (5.4), а относительная проницаемость эмболизата – по формуле (5.6). Линиями показаны графики приближенных функций Баклея–Леверетта, полученных по формуле (5.7), которые используются во всех дальнейших расчетах.

6. ОПТИМАЛЬНЫЙ СЦЕНАРИЙ ЭМБОЛИЗАЦИИ

Оптимальные сценарии эмболизации найдены численно путем расчета значений целевого функционала (4.1) и ограничений (2.11), (4.2), (4.3). При расчетах предполагалось, что структура АВМ достаточно однородна и в силу этого величины $K$, ${{k}_{b}}$ не зависит от $x$, кроме того предполагалось, что ${{R}_{{13}}} \ll {{R}_{{12}}}$ и сопротивлением подводящего сосуда можно пренебречь, т.е. ${{R}_{{13}}} = 0$. Значение критического давления ${{p}_{{{\text{max}}}}}$ выбрано равным $80$ мм рт.ст. Параметры модели АВМ, включенной в систему церебрального кровообращения для двух пациентов, приведены в табл. 1. Чтобы изучить поведение целевого функционала и ограничений для двух пациентов в зависимости от значений параметров ${{\tau }_{1}}$, ${{\tau }_{2}}$ и $\gamma $, рассматривается сетка в области параметров $[\tau _{1}^{l},\tau _{1}^{r}] \times [\tau _{2}^{l},\tau _{2}^{r}] \times [{{\gamma }^{l}},{{\gamma }^{r}}]$. Также изучается вопрос о влиянии на процесс эмболизации исходного распределения кровотока между здоровым сосудом и АВМ, для этого все расчеты проведены для трех различных значений потока через здоровый сосуд до эмболизации, когда ${{Q}_{{12}}}(0) = {{k}_{{12}}}Q(0)$, где ${{k}_{{12}}}$ принимает значения 0.5, 1, 1.5. Характеристики выбранной сетки параметров приведены в табл. 2. Для нахождения экстремума функционала и проверки ограничений использовались сеточные данные, а для визуализации результатов использовалась интерполяция сеточных данных.

На фиг. 3 представлен расчет для пациента $1$ при ${{k}_{{12}}} = 1$. Данная картина является типичной и для остальных случаев. Область допустимых значений параметров лежит выше поверхности ограничений $OPA{{V}_{1}}{{V}_{2}}$, составленной из поверхности критического давления $OPA$ и поверхности венозного ограничения $OA{{V}_{1}}{{V}_{2}}$, две эти поверхности пересекаются по кривой $AO$. Также на фиг. 3 изображена одна из поверхностей уровня ${{J}_{1}}{{J}_{2}}{{J}_{3}}$ целевого функционала $J$, соответствующее ей значение функционала больше оптимального, вследствие чего в область допустимых параметров попадают точки этой поверхности, лежащие в криволинейном треугольнике ${{O}_{1}}{{O}_{2}}{{O}_{3}}$. При уменьшении значения функционала до оптимального в данном случае весь криволинейный треугольник переходит в точку $O$.

Результаты всех расчетов приведены в табл. 3, в которой указаны оптимальные значения параметров и соответствующие им значения целевого функционала. Наблюдаются три качественно различных случая оптимальной эмболизации, которые переходят из одного в другой при изменении расхода крови через здоровый сосуд. Эти случаи проиллюстрированы на фиг. 4, где показаны результаты расчетов на плоскости (${{\tau }_{1}},T$), которая представляет собой сечение пространства параметров (${{\tau }_{1}},T,\gamma $) плоскостью оптимального значения $\gamma $. Область $OCDE$ – область венозного ограничения, граница этой области указана пунктирной линией, а область $OBC$ – область превышения критического давления, граница этой области указана штриховой линией. Значения линий уровня функционала отмечены шкалой справа от фигуры, а белая точка отмечает оптимальное значение параметров ${{\tau }_{1}}$, ${{\tau }_{2}} = T - {{\tau }_{1}}$. В первых двух случаях (фиг. 4а, б) активными являются венозное ограничение (4.3), а также ограничение, соответствующее достижению критического давления (4.2). В этих случаях оптимальное значение функционала достигается на кривой $OA$ (фиг. 3) – кривой пересечения границы области венозного ограничения и поверхности критического давления.

В первом случае (фиг. 4а) кривая $OA$ лежит на весьма малом расстоянии от поверхности уровня целевого функционала, соответствующей режиму оптимальной эмболизации – значения целевого функционала в точках кривой $OA$ различаются на величину порядка 1%. В этом случае с практической точки зрения в качестве оптимального режима эмболизации могут использоваться все точки кривой $OA$, соответствующие как разрывному режиму эмболизации – для граничной точки $O$ (4.5), так и непрерывным кусочно-линейным режимам – для остальных точек кривой (4.4). Данный случай соответствует значениям ${{k}_{{12}}} = 0.5$ для пациента $1$ и пациента $2$. На фиг. 4а изображен один из оптимальных режимов эмболизации с параметрами${{\tau }_{1}} = 4.3$, ${{\tau }_{2}} = 1.55$, $\gamma = 0.175$.

Во втором случае (фиг. 4б) кривая $OA$ пересекается с поверхностью уровня целевого функционала, соответствующей оптимальному режиму эмболизации в единственной точке $O$, в результате чего оптимальный режим эмболизации является разрывным и задается формулой (4.5). Данный случай реализуется для пациента $1$ при ${{k}_{{12}}} = 1,\;1.5$ и для пациента $2$ при ${{k}_{{12}}} = 1$.

В третьем случае (фиг. 4в) для пациента $2$ при ${{k}_{{12}}} = 1.5$ не достигается значение критического давления для всех рассмотренных значений параметров, а оптимальный режим реализуется при максимально допустимом расходе эмболизата, определяемый ограничением (2.11). Точка, соответствующая оптимальному режиму эмболизации, лежит на границе области допустимых параметров, и режим эмболизации также является разрывным (4.5). В этом случае достигается самая маленькая степень содержания крови в АВМ в конце операции.

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

7. ЗАКЛЮЧЕНИЕ

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

Для описания движения эмболизата в теле АВМ используется уравнение двухфазной фильтрации Баклея–Леверетта, в котором функции Баклея–Леверетта находятся для двух пациентов по данным интраоперационных измерений. В модели учитывается перераспределение крови в здоровые сосуды в процессе эмболизации.

Для численного исследования поставленной задачи разработана модификация схемы CABARET.

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

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

Исследование проводилось в соответствии с Хельсинкской декларацией и было одобрено Инспекционной комиссией Клиники им. Мешалкина (заседание № 549). Все включенные пациенты дали информированное согласие.