Акустический журнал, 2022, T. 68, № 1, стр. 3-13

Компенсация искажений фокусированных ультразвуковых пучков при транскраниальном облучении головного мозга на различной глубине

Д. Д. Чупова a*, П. Б. Росницкий a, Л. Р. Гаврилов b**, В. А. Хохлова a

a Московский государственный университет им. М.В. Ломоносова
119991 Москва, ГСП-1, Ленинские горы, Россия

b АО “Акустический институт им. акад. Н.Н. Андреева”
117036 Москва, ул. Шверника 4, Россия

* E-mail: daria.chupova@yandex.ru
** E-mail: lrgavrilov.1938@mail.ru

Поступила в редакцию 27.08.2021
После доработки 21.09.2021
Принята к публикации 22.09.2021

Полный текст (PDF)

Аннотация

Проанализированы возможности компенсации аберраций при фокусировке ультразвукового пучка через кости черепа с использованием мозаичных решеток с радиусом кривизны и апертурой F = D = 200 мм, частотой 1 МГц и абсолютно плотным рандомизированным заполнением поверхности элементами. Рассматривается влияние количества элементов (256, 512 и 1024) и глубины фокусировки (25–65 мм от внутренней поверхности черепа) на качество компенсации аберраций, т.е. на остроту фокусировки, местоположение фокуса и максимальное давление в нем. Акустическая модель головы человека строится на основе данных магнитно-резонансной томографии (МРТ). Расчет поля и компенсации аберраций проводятся на основе интеграла Рэлея и волнового уравнения в модели Кельвина–Фойгта. Продемонстрирована возможность острой фокусировки с помощью рассматриваемых решеток в указанном интервале глубин, при этом ширина фокальной перетяжки составляет около 2 мм по уровню – 6 дБ. Проанализирован относительный вклад различных волновых эффектов в искажение ультразвукового пучка при прохождении через череп. Показано, что наиболее сильный вклад в ослабление пучка вносят аберрации (7.4 дБ) и поглощение (6.7 дБ). Вклад отражения (2.1 дБ) и генерации сдвиговых волн в черепе (2 дБ) менее существенен.

Ключевые слова: медицинская акустика, фокусированный ультразвук высокой интенсивности, многоэлементные решетки, интеграл Рэлея

ВВЕДЕНИЕ

Использование мощного фокусированного ультразвука для лечения заболеваний мозга на сегодняшний день является одним из наиболее успешных клинических направлений развития неинвазивной ультразвуковой хирургии [16]. При этом с помощью фокусировки ультразвукового пучка через кости черепа в заданные участки мозга реализуется тепловой механизм воздействия: заранее намеченные патологические участки разрушаются путем быстрого локального перегрева. Проведение таких операций без прямого хирургического вмешательства существенно снижает риски занесения инфекций и повреждения здоровых тканей мозга [7]. Однако присутствие костей черепа вносит искажения в пространственную структуру ультразвукового пучка, поэтому для их компенсации и обеспечения острой фокусировки ультразвука используются многоэлементные фазированные решетки, на каждом элементе которых возможно независимо задавать амплитуду и фазу [7, 8]. Для определения нужной фазы на элементах решетки используются как лучевые методы [9], так и дифракционные, основанные на обращении волнового фронта [10] либо инвертировании фазы [11]. В нейрохирургической практике успешно используются 1024-элементные решетки клинических систем ExAblate (InSightec Ltd., Tirat Carmel, Israel), имеющие форму полусферы с радиусом кривизны 150 мм, апертурой 300 мм и рабочей частотой 650–720 кГц [8]. Такая конфигурация решетки позволяет облучать центральные участки мозга в области таламуса и проводить лечение эссенциального тремора, тремора, вызванного болезнью Паркинсона, локализованных опухолей и других заболеваний мозга [16]. В известной нам литературе отсутствуют количественные данные о размере возможной области сканирования в этом случае.

Недавно был предложен новый класс многоэлементных решеток с более компактной формой в виде сегмента сферы с углом схождения $\alpha = 2\arcsin \left( {{D \mathord{\left/ {\vphantom {D {2F}}} \right. \kern-0em} {2F}}} \right)$ = 60°, радиусом кривизны и апертурой F = D = 200 мм, а также более высокой рабочей частотой f = 1 МГц. Активная площадь решетки и, соответственно, излучаемая мощность приблизительно вдвое меньше, чем у полусферической решетки, однако дефицит мощности удается частично скомпенсировать путем использования абсолютно плотного рандомизированного заполнения поверхности решетки элементами с мозаичной структурой (рис. 1а) [1214]. Отметим, что другими авторами предлагалось несколько вариантов терапевтических решеток с близкими углами схождения и мозаичным расположением элементов, однако в них не обеспечивалась максимально достижимая площадь активной поверхности решетки, которая реализуется при использовании рассматриваемой в данной работе модели [15, 16].

Рис. 1.

Схемы исследуемых решеток: (а) 256 элементов, (б) 512 элементов, (в) 1024 элемента, (г) идеализированная решетка. (д) – Схемы расчета поля при облучении мозга и (е) – компенсации аберраций. Стрелками показан порядок проведения расчетов: с использованием интеграла Рэлея и модели Кельвина–Фойгта (К.–Ф.). Фокусировка проводилась на различной глубине (д): центр кривизны решетки смещен относительно центра мозга на –30 мм (1); –20 мм (2); –10 мм (3); 0 мм (4); +10 мм (5).

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

Кроме того, в предыдущих работах было теоретически показано, что для подобной решетки с частотой 1 МГц и 256 элементами при компенсации аберраций, вносимых черепом, возможно достижение давлений, необходимых для образования высокоамплитудных ударных фронтов и механического разрушения тканей с помощью метода гистотрипсии с кипением [12, 13]. При использовании данного метода практически не проявляются тепловые эффекты, что снижает риск перегрева черепа [1719]. Однако для решеток с иным числом элементов и областей мозга, лежащих за пределами таламуса, исследования по компенсации аберраций не проводились даже в приближении линейной фокусировки пучка через кости черепа.

Целью данной работы было теоретическое исследование влияния количества элементов в предложенных решетках и глубины фокусировки в мозге при механическом перемещении решетки на качество компенсации аберраций в сходящемся линейном пучке. Основными критериями при оценке качества компенсации аберраций являются острота фокусировки (ширина и длина фокального пика), правильность местоположения фокуса и уровень максимального давления в нем. Проведена также оценка относительного вклада различных волновых эффектов в искажение структуры пучка и максимально достижимой амплитуды волны. Были рассмотрены эффекты генерации сдвиговых волн в черепе, поглощения, отражения и аберраций, возникающих при прохождении пучка через неоднородные по толщине кости черепа. В численном эксперименте использовалась акустическая модель головы человека, построенная на основе данных магнитно-резонансной томографии (МРТ). Расчет поля решетки и компенсация аберраций проводились на основе разработанного ранее метода, объединяющего различные волновые модели [12].

ТЕОРЕТИЧЕСКАЯ МОДЕЛЬ

Рассмотрим модель многоэлементной решетки в виде сферического сегмента с фокусным расстоянием и апертурой F = D = 200 мм, углом схождения 60° и рабочей частотой f = 1 МГц [13]. С помощью разработанного ранее алгоритма, основанного на использовании мозаики с ячейками равной площади, поверхность решетки разбивалась на элементы, имеющие форму многоугольников, расположенных случайным образом на поверхности решетки [14]. При этом достигалось абсолютно плотное заполнение поверхности излучателя элементами. Такой способ разбиения одновременно позволяет максимизировать активную площадь решетки заданной геометрии и минимизировать побочные дифракционные эффекты, которые возникают при периодическом расположении элементов [20]. Для технической реализуемости модели излучателя был введен зазор между элементами, равный 0.5 мм. С использованием разработанного алгоритма были построены модели решеток с различным количеством элементов: 256, 512 и 1024 (рис. 1а–1в). Коэффициент заполнения решеток с учетом зазора между элементами составил 92% (площадь элемента 121 мм2) для 256-элементной решетки, 88% (площадь элемента 58 мм2) для 512-элементной решетки и 84% (площадь элемента 28 мм2) для 1024-элементной решетки. Была также рассмотрена идеализированная модель решетки с условно бесконечным количеством элементов (рис. 1г). Для ее построения расчетная прямоугольная сетка численного моделирования проектировалась на сферическую поверхность излучателя. Таким образом, была достигнута возможность квазинепрерывно изменять амплитуду и фазу на поверхности излучателя [21].

Для моделирования фокусировки ультразвукового пучка через интактный череп использовалась трехмерная акустическая модель головы человека, построенная на основе данных МРТ [12]. Модель головы была представлена набором из 91 снимка в аксиальной проекции, взятых на разной высоте, размером 191 × 256 пикселей каждый. Пространственное разрешение снимков составляло 1 × 1 × 1 мм. Каждый снимок разделялся методом пороговой обработки на четыре сегмента: внешнее пространство, заполненное водой, кожу, череп и мозг. В модели были учтены геометрические особенности каждого сегмента, представляющего голову человека, при этом каждая среда считалась однородной [22, 23] и характеризовалась следующими параметрами [12]: плотность ${{{{\rho }}}_{0}}$ [24], скорость распространения ${{c}_{p}}$ [25] и коэффициент поглощения ${{{{\alpha }}}_{p}}$ продольных волн [24], а также скорость распространения ${{c}_{s}}$ и коэффициент поглощения ${{{{\alpha }}}_{s}}$ сдвиговых волн (отличны от нуля только для черепа) [24, 25]. Параметры всех сегментов задавались в программе в виде матрицы, их значения приведены в табл. 1.

Таблица 1.  

Значения скорости распространения продольных ${{c}_{p}}$ и сдвиговых ${{c}_{s}}$ волн, плотности ${{{{\rho }}}_{0}}$, коэффициентов поглощения продольных ${{{{\alpha }}}_{p}}$ и сдвиговых ${{{{\alpha }}}_{s}}$ волн для воды, кожи, черепа и мозга

  ${{c}_{p}}$, м/с ${{c}_{s}}$, м/с ${{{{\rho }}}_{0}}$, кг/м3 ${{{{\alpha }}}_{p}}$, дБ/см ${{{{\alpha }}}_{s}}$, дБ/см
Вода 1500 0 1000 0 0
Кожа 1624 0 1109 1.84 0
Череп 2820 1500 1732 8.83 19.15
Мозг 1550 0 1030 0.21 0

Численный алгоритм моделирования фокусировки ультразвука состоял из комбинации двух методов. На первом этапе акустическое поле излучателя рассчитывалось с помощью интеграла Рэлея [26] в однородной среде (воде) от поверхности решетки до горизонтальной плоскости (xy), расположенной вблизи поверхности черепа (рис. 1д):

(1)
$p({\mathbf{r}}) = - \frac{{i{{\omega }}{{{{\rho }}}_{0}}}}{{2{{\pi }}}}\int\limits_S {\frac{{{{v}_{n}}({\mathbf{r}}{\kern 1pt} ')\exp (ikR)}}{R}} dS{\kern 1pt} ',$
где ${{\omega }} = 2{{\pi }}f$ – циклическая частота, $k = {{{\omega }} \mathord{\left/ {\vphantom {{{\omega }} {{{c}_{0}}}}} \right. \kern-0em} {{{c}_{0}}}}$ – волновое число, ${{c}_{0}}$ – скорость звука в среде, ${{{{\rho }}}_{0}}$ – плотность среды, ${{v}_{n}}({\mathbf{r}}{\kern 1pt} ')$ – комплексная амплитуда колебательной скорости на плоскости xy, ${\mathbf{r}}{\kern 1pt} '$ – радиус вектор элемента поверхности $dS'$, $S$ – площадь плоскости xy, $R = \left| {{\mathbf{r}} - {\mathbf{r}}{\kern 1pt} '} \right|$ – расстояние от элемента площади $dS'$ до точки наблюдения с координатой ${\mathbf{r}}$. Рассматривалась компонента комплексной амплитуды поля с временной зависимостью $\exp ( - i{{\omega }}t)$, где $i$ – мнимая единица. Распределение амплитуды колебательной скорости на поверхности элементов предполагалось равномерным: ${{v}_{n}}({\mathbf{r}}{\kern 1pt} ')$ = ${{v}_{0}}$.

Прямой численный расчет интеграла Рэлея (1) с поверхности многоэлементной решетки является достаточно затратным по времени, особенно для проведения серий вычислений. В условиях данной задачи расстояние между излучателем и заданной плоскостью значительно превышает ближнюю зону элемента решетки, поэтому для ускорения расчетов использовался аналитический метод, согласно которому акустическое поле на плоскости xy вычислялось как сумма аналитических решений для дальнего поля каждого элемента [13, 14, 20]. Для этого предварительно многоугольные элементы излучателя разбивались на прямоугольные треугольники. В свою очередь, поле треугольного излучателя рассчитывалось по формуле [14]:

(2)
$p = \frac{{{{p}_{0}}ab\exp (ik{{r}_{0}})[I(a,x) - I(b,y)]}}{{2{{\pi }}{{r}_{0}}({{ax} \mathord{\left/ {\vphantom {{ax} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}} - {{by} \mathord{\left/ {\vphantom {{by} {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}})}},$
где $I(a,x) = \exp ({{ - ika} \mathord{\left/ {\vphantom {{ - ika} {2({x \mathord{\left/ {\vphantom {x {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}})}}} \right. \kern-0em} {2({x \mathord{\left/ {\vphantom {x {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}})}})sinc({{ka} \mathord{\left/ {\vphantom {{ka} {2({x \mathord{\left/ {\vphantom {x {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}})}}} \right. \kern-0em} {2({x \mathord{\left/ {\vphantom {x {{{r}_{0}}}}} \right. \kern-0em} {{{r}_{0}}}})}})$, $a$ и $b$ катеты прямоугольного треугольника, ${{p}_{0}} = {{{{\rho }}}_{0}}{{{\text{c}}}_{0}}{{v}_{n}}$ характерное давление на поверхности элемента, ${{r}_{0}} = {{({{x}^{2}} + {{y}^{2}} + {{z}^{2}})}^{{1/2}}}$ радиус-вектор точки наблюдения, $\left\{ {x,y,z} \right\}$ координаты точки наблюдения. Таким образом, поле акустического давления на плоскости xy, перпендикулярной оси решетки и расположенной вблизи поверхности черепа, было найдено как сумма полей подэлементов в форме прямоугольных треугольников, на которые были разбиты элементы излучателя.

Полученный результат использовался в качестве граничного условия для дальнейшего расчета поля в области головы человека. Для этого применялся численный псевдоспектральный метод решения волнового уравнения в модели Кельвина–Фойгта, реализованный с помощью программы k-Wave (www.k-wave.org) [22, 23]:

(3)
$\begin{gathered} {{{{\rho }}}_{0}}\frac{{{{\partial }^{2}}{\mathbf{u}}}}{{\partial {{t}^{2}}}} = ({{\lambda }} + {{\mu }})\nabla (\nabla \cdot {\mathbf{u}}) + {{\mu }}{{\nabla }^{2}}{\mathbf{u}} + \\ + \,\,({{\chi }} + {{\eta }})\nabla \left( {\nabla \cdot \frac{{\partial {\mathbf{u}}}}{{\partial t}}} \right) + {{\eta }}{{\nabla }^{2}}\frac{{\partial {\mathbf{u}}}}{{\partial t}}, \\ \end{gathered} $
где ${\mathbf{u}}$ – вектор смещения частиц среды, ${{\lambda }}$ и ${{\mu }}$ – параметры Ламе, ${{\mu }}$ – модуль сдвига, ${{\chi }}$ и ${{\eta }}$ – коэффициенты вязкости при сжатии и сдвиге, соответственно. Параметры Ламе связаны со скоростью продольных ${{c}_{p}}$ и сдвиговых ${{c}_{s}}$ волн следующим образом: $c_{s}^{2} = {{{\mu }} \mathord{\left/ {\vphantom {{{\mu }} {{{{{\rho }}}_{0}}}}} \right. \kern-0em} {{{{{\rho }}}_{0}}}}$, $c_{p}^{2} = {{({{\lambda }} + 2{{\mu }})} \mathord{\left/ {\vphantom {{({{\lambda }} + 2{{\mu }})} {{{{{\rho }}}_{0}}}}} \right. \kern-0em} {{{{{\rho }}}_{0}}}}$; коэффициенты поглощения ${{{{\alpha }}}_{s}}$ и ${{{{\alpha }}}_{p}}$ связаны с коэффициентами вязкости: ${{{{\alpha }}}_{s}} = {{{{\eta }}{{{{\omega }}}^{2}}} \mathord{\left/ {\vphantom {{{{\eta }}{{{{\omega }}}^{2}}} {(2{{{{\rho }}}_{0}}c_{s}^{3})}}} \right. \kern-0em} {(2{{{{\rho }}}_{0}}c_{s}^{3})}}$ и ${{{{\alpha }}}_{p}} = {{({{\chi }} + 2{{\eta }}){{{{\omega }}}^{2}}} \mathord{\left/ {\vphantom {{({{\chi }} + 2{{\eta }}){{{{\omega }}}^{2}}} {(2{{{{\rho }}}_{0}}c_{p}^{3})}}} \right. \kern-0em} {(2{{{{\rho }}}_{0}}c_{p}^{3})}}$.

Данная модель (3) учитывает эффекты дифракции, поглощения, неоднородности и генерации сдвиговых волн в черепе. В результате комбинации двух описанных методов (2) и (3) находилось распределение амплитуды акустического давления pA/p0, нормированной на характерную амплитуду давления p0 на поверхности решетки, в области головы в трех основных анатомических плоскостях: сагиттальной (zx), фронтальной (zy) и аксиальной (xy), а также одномерные распределения pA/p0 вдоль оси решетки. Для анализа уровней давления, достигаемых вблизи поверхности черепа, рассчитывались распределения амплитуд давления pA/p0 в объеме прямоугольной формы, грани которого были удалены от поверхности черепа на 10 мм, причем верхняя грань находилась в воде, а нижняя – в мозге. Для моделирования использовалась расчетная сетка с пространственным шагом $\Delta x = \Delta y = \Delta z = 0.5$ мм. Так как размер вокселя трехмерной акустической модели головы превышал в два раза шаг сетки, каждый воксель делился пополам в каждом из трех направлений. При этом тип среды каждого вокселя оставался прежним. Шаг по времени расчетной сетки выбирался в соответствии с критерием Куранта–Фридрихса–Леви и удовлетворял условию $CFL \equiv {{{{c}_{{\max }}}\Delta t} \mathord{\left/ {\vphantom {{{{c}_{{\max }}}\Delta t} {\Delta x}}} \right. \kern-0em} {\Delta x}}$, где $CFL = 0.1$, а ${{c}_{{\max }}}$ – максимальная скорость звука в рассматриваемой модели.

Для компенсации аберраций, вызванных присутствием костей черепа, использовались аналогичные методы расчетов, но в обратном порядке (рис. 1е). Вначале использовалась модель Кельвина–Фойгта (3) для моделирования распространения сферической волны из точки фокуса до плоскости xy вблизи черепа. Далее использовалось численное решение интеграла Рэлея (1) в однородной среде для расчета поля в плоскости xy и получения комплексных амплитуд и фаз в геометрическом центре каждого элемента решетки.

При расчете поля с компенсацией аберраций полученные фазы инвертировались, а в расчетах без компенсации фаза на каждом элементе полагалась равной нулю. Для достижения максимальной интенсивности в фокусе амплитуды на элементах устанавливались одинаковыми для всех элементов решетки [7]. Для решеток с различным количеством элементов и, соответственно, с различной активной площадью, амплитуда на их элементах выбиралась обратно пропорциональной коэффициенту заполнения решетки так, чтобы все излучатели обеспечивали одинаковую амплитуду давления в фокусе при фокусировке в воде. Такая нормировка позволяет “в чистом виде” выявить влияние числа элементов на качество компенсации аберраций. Далее за нормировочную константу p0 принимается характерное начальное давление p0 = c0ρ0v0 на поверхности 256-элементной решетки.

Расчеты проводились для четырех решеток: 256-элементной, 512-элементной, 1024-элементной и идеализированной решетки, на пяти различных глубинах фокусировки c компенсацией аберраций и без нее. Глубина фокусировки изменялась путем механического перемещения решетки вдоль собственной оси и составляла 25, 35, 45, 55 и 65 мм (рис. 1д, 1–5), считая от внутренней поверхности черепа. Центр мозга соответствует значению 55 мм (рис. 1д, 4).

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

Для оценки вклада генерации сдвиговых волн в черепе скорость ${{c}_{s}}$ и коэффициент поглощения ${{{{\alpha }}}_{s}}$ сдвиговых волн в матрице с акустическими параметрами сред устанавливались равными нулю.

Для оценки вклада поглощения продольных волн каждый элемент 256-элементного излучателя соединялся с точкой фокуса, и для каждой поглощающей среды (кожи, кости, мозга) находилось среднее расстояние $\left\langle {{{l}^{k}}} \right\rangle = {{N}^{{ - 1}}}\sum\nolimits_{i = 1}^N {l_{i}^{k}} $ вдоль получаемых лучей. Здесь N = 256, k = “к”, “ч”, “м”, что обозначает разные среды “кожа”, “кости черепа”, “ткани мозга”. Данные значения расстояний использовались для оценки снижения амплитуды давления A за счет поглощения:

(4)
$A = \exp ( - {{\alpha }}_{p}^{{\text{к}}}\left\langle {{{l}^{{\text{к}}}}} \right\rangle - {{\alpha }}_{p}^{{\text{ч}}}\left\langle {{{l}^{{\text{ч}}}}} \right\rangle - {{\alpha }}_{p}^{{\text{м}}}\left\langle {{{l}^{{\text{м}}}}} \right\rangle ).$

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

Оценка вклада отражения была проведена с использованием приближения плоскопараллельного слоя для случая нормального падения. Каждая среда сохраняла свой импеданс, но граница раздела сред считалась плоской. Таким образом, суммарный коэффициент отражения R по давлению вычислялся как произведение коэффициентов отражения на границах вода–кожа: ${{R}_{{{\text{вк}}}}} = {{\left( {{{{{\rho }}}_{{\text{к}}}}{{c}_{{\text{к}}}} - {{{{\rho }}}_{{\text{в}}}}{{c}_{{\text{в}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{{{\rho }}}_{{\text{к}}}}{{c}_{{\text{к}}}} - {{{{\rho }}}_{{\text{в}}}}{{c}_{{\text{в}}}}} \right)} {\left( {{{{{\rho }}}_{{\text{к}}}}{{c}_{{\text{к}}}} + {{{{\rho }}}_{{\text{в}}}}{{c}_{{\text{в}}}}} \right)}}} \right. \kern-0em} {\left( {{{{{\rho }}}_{{\text{к}}}}{{c}_{{\text{к}}}} + {{{{\rho }}}_{{\text{в}}}}{{c}_{{\text{в}}}}} \right)}}$, кожа–череп ${{R}_{{{\text{кч}}}}}$, череп–мозг ${{R}_{{{\text{чм}}}}}$:

(5)
$R = {{R}_{{{\text{вк}}}}}{{R}_{{{\text{кч}}}}}{{R}_{{{\text{чм}}}}}.$

Оценка влияния аберраций проводилась с учетом полученных результатов для сдвиговых волн, поглощения и отражения, приводящих к ослаблению поля в фокальной области пучка. Для этого величина максимума давления, достигаемого при фокусировке в центр мозга без компенсации аберраций pA/p0, делилась на ослабления амплитуды давления за счет генерации сдвиговых волн, поглощения и отражения. Отношение полученной амплитуды давления и максимума давления в воде принималoсь за вклад аберраций в искажение пучка.

РЕЗУЛЬТАТЫ

На рис. 2 представлены двумерные распределения безразмерной амплитуды давления pA/p0, нормированной на начальное характерное значение p0 на поверхности 256-элементной решетки, при облучении центра мозга без компенсации аберраций (рис. 1д, 4). Значение z = 0 соответствует внутренней поверхности черепа. Результаты представлены в трех основных анатомических плоскостях, проходящих через центр кривизны решетки: сагиттальной, разделяющую голову на правую и левую половины (zx), фронтальной – на переднюю и заднюю (zy) и аксиальной – на выше- и нижележащие (xy). В отсутствие коррекции фазы акустическое поле сильно искажено за счет присутствия костей черепа, максимум давления смещен относительно центра кривизны излучателя и составляет всего pA/p0 = 12.5, тогда как при фокусировке в воде коэффициент усиления по давлению в фокусе равен 100. Наблюдаются также дополнительные максимумы давления вблизи черепа и “размытие” фокуса.

Рис. 2.

Распределение амплитуды давления pA/p0, нормированной на амплитуду давления p0, на элементах решетки, в плоскостях (а) xz, (б) zy и (в) xy для случая облучения центра мозга с помощью 256-элементной решетки без компенсации аберраций.

Влияние аберраций при различной глубине фокусировки проиллюстрировано на рис. 3. Для изменения глубины фокусировки излучатель механически смещался вдоль собственной оси на ‒10 мм вниз, на +10, +20 и +30 мм вверх относительно центра мозга (рис. 1д, 5, 3, 2, 1). Результаты расчетов показаны в виде двумерных распределений pA/p0 в плоскости zy. Для всех глубин положения фокуса наблюдается сильное искажение ультразвукового пучка, “размытие” фокальной области и появление дополнительных максимумов давления вблизи черепа. При уменьшении глубины фокусировки эти эффекты усиливаются. На бóльших глубинах фокусировки, –10 и +10 мм относительно центра мозга (рис. 3в, 3г), амплитуда давления вблизи черепа составляет от 60 до 70% от максимального значения pA/p0. При фокусировке и облучении мозга ближе к черепу, на глубине 35–25 мм от его внутренней поверхности (рис. 3а и 3б), амплитуда давления вблизи черепа становится сравнимой с амплитудой в фокусе. Максимум давления и интенсивности в мозге в случае крайнего положения фокусировки (рис. 3а) снижается на 1.3 дБ (соответственно, на 14 и 26%) по сравнению со случаем облучения центра мозга. Таким образом, уменьшение глубины фокусировки и облучение структур, близких к внутренней поверхности черепа, приводит к возрастанию искажений ультразвукового пучка, уменьшению амплитуды давления в фокусе и ее увеличению вблизи кости черепа.

Рис. 3.

Распределения нормированной амплитуды давления pA/p0 в плоскости zy при фокусировке в мозге на различной глубине без компенсации аберраций для 256-элементной решетки. Центр кривизны решетки смещен относительно центра мозга на (а) –30 мм; (б) –20 мм; (в) –10 мм; (г) +10 мм.

Качество коррекции аберраций при различной глубине фокусировки проиллюстрировано на рис. 4, где представлены двумерные распределения амплитуды безразмерного акустического давления pA/p0 в плоскости zy. Коррекция аберраций была проведена для 256-элементной решетки при облучении на тех же глубинах, что и на рис. 3. Видно, что коррекция позволяет значительно улучшить качество фокусировки, а именно восстановить положение фокуса, обеспечить в фокусе узкую фокальную перетяжку и почти двукратное увеличение амплитуды давления. Так, в случае облучения центра мозга давление в фокусе составляет pF/p0 = 23.8. На одномерных распределениях вдоль оси решетки (рис. 5) можно более детально увидеть результат коррекции на разных глубинах. Ширина и длина фокальной перетяжки на всех глубинах фокусировки практически не отличаются от случая фокусировки в воде и составляют соответственно около 2 и 14 мм по полувысоте фокального пика. Однако с уменьшением глубины фокусировки наблюдается снижение давления и интенсивности в фокусе; при поднятии решетки на 30 мм относительно центра мозга (рис. 4а) – на 1.6 дБ (соответственно, на 17 и 30%). При приближении фокуса к поверхности черепа также увеличивается уровень побочных максимумов давления вблизи черепа от 40% (в случае облучения центра мозга) до 76% (в крайнем положении) от значения pF/p0. В терминах интенсивности это составляет, соответственно, 16 и 58%.

Рис. 4.

Распределения нормированной амплитуды давления pA/p0 в плоскости zy при фокусировке на различной глубине в мозге c компенсацией аберраций для 256-элементной решетки. Центр кривизны решетки смещен относительно центра мозга на (а) –30 мм; (б) –20 мм; (в) –10 мм; (г) +10 мм.

Рис. 5.

Распределения нормированной амплитуды давления pA/p0 (а) вдоль оси решетки и (б) в фокальной плоскости при фокусировке на различной глубине в мозге c компенсацией аберраций для 256-элементной решетки. Центр кривизны решетки смещен относительно центра мозга на –30 мм (1); –20 мм (2); –10 мм (3); 0 мм (4); +10 мм (5).

Влияние количества элементов на эффективность компенсации аберраций проиллюстрировано на рис. 6, где представлен график зависимости величины pA/p0 от глубины фокусировки, считая от внутренней поверхности черепа, для решеток с 256, 512, 1024 элементами и идеализированной решетки (рис. 1а–1г). В верхней части графика представлена нормированная амплитуда давления в фокусе излучателя pF/p0 для идеализированной решетки (сплошная кривая), 1024-элементной решетки (пунктир), 512-элементной решетки (штрих-пунктир) и 256-элементной решетки (штриховая кривая). Кривые в нижней части графика представляют значения максимума давления вблизи поверхности черепа pч/p0. Цифры 15 соответствуют глубинам фокусировки, отмеченным на рис. 1д. Напомним, что значения pA здесь нормированы на характерное значение p0 на поверхности 256-элементной решетки с тем расчетом, чтобы все решетки обеспечивали одинаковую амплитуду давления в фокусе при фокусировке в воде.

Рис. 6.

Значения нормированной амплитуды давления в фокусе pF/p0 и максимума амплитуды давления вблизи поверхности черепа pч/p0 для решеток с 1024 элементами (пунктир), 512 элементами (штрих-пунктир), 256 элементами (штриховая линия) и идеализированной решетки (сплошная линия) при фокусировке на различных глубинах: 25 (1), 35 (2), 45 (3), 55 (область таламуса, 4) и 65 мм (5) от внутренней поверхности черепа.

Из графика видно, что с увеличением количества элементов излучателя растет амплитуда давления в фокусе. Так, в случае идеализированной решетки давление и интенсивность увеличиваются на 1.9 дБ (соответственно, на 24 и 55%) по сравнению со случаем 256-элементной решетки и на 0.62 дБ (соответственно, на 7 и 15%) по сравнению с 1024-элементным излучателем. При переходе от 256 к 512 элементам рост давления и интенсивности в фокусе составляет 0.67 дБ (8 и 17%), от 512 к 1024 – 0.57 дБ (7 и 14%). При уменьшении глубины фокусировки амплитуда давления в фокусе снижается приблизительно одинаково для всех излучателей: амплитуды давления в крайних положениях фокуса по глубине (1 и 5 на рис. 6) отличаются на 17%. Максимальные уровни давления вблизи поверхности черепа с увеличением количества элементов изменяются незначительно. При этом амплитуда давления вблизи черепа увеличивается в 1.68 раза при уменьшении глубины фокусировки в рассматриваемом интервале (положения 1 и 5).

Рис. 7 позволяет сравнить форму фокального пика и достижимые уровни давления для решеток с разным количеством элементов. Представлены одномерные распределения нормированной амплитуды давления pA/p0 вдоль оси решетки z и вдоль оси x в фокальной плоскости при фокусировке в центр мозга для решеток с различным количеством элементов и идеализированной решетки. Значение z = 0 соответствует внутренней поверхности черепа. Из распределений видно, что структура поля, продольные и поперечные размеры фокальной области всех излучателей, определяемые по полувысоте и первым нулям давления относительно фокуса, практически одинаковы для всех излучателей. Однако увеличение числа элементов решетки позволяет повысить абсолютный уровень давления в фокусе излучателя.

Рис. 7.

Распределения нормированной амплитуды давления pA/p0 (а) вдоль оси решетки и (б) в фокальной плоскости для решеток: 256 элементов (штриховая линия), 512 элементов (штрих-пунктир), 1024 элемента (пунктир), идеализированный случай (сплошная линия) в случае облучения центра мозга с компенсацией аберраций.

При поэтапном анализе влияния различных волновых эффектов на уменьшение максимально достижимых уровней давления в фокальной области пучка были получены следующие результаты. На первом шаге, при отключении генерации сдвиговых волн, максимум амплитуды давления поля увеличился на 2 дБ, при этом смещение относительно геометрического фокуса излучателя, побочные максимумы и “размытие фокуса” остались прежними. Дополнительное отключение поглощения привело к увеличению амплитуды давления на 6.7 дБ. При отключении отражения максимум по давлению увеличился еще на 2.1 дБ, а при последующей компенсации аберраций – на 7.4 дБ. Таким образом, проведенные оценки показали, что наибольший вклад в ослабление пучка вносят аберрации; несколько меньший, но сравнимый с ним вклад дает поглощение. Эффекты отражения и генерации сдвиговых волн в черепе влияют на распространение пучка менее существенно.

ЗАКЛЮЧЕНИЕ

В работе проанализированы возможности коррекции аберраций при различной глубине фокусировки в структурах головного мозга для предложенного нового типа решеток с различным количеством элементов. Показано, что при транскраниальном воздействии на структуры мозга с помощью предложенных решеток компенсация аберраций позволяет значительно улучшить качество фокусировки и добиться узкой фокальной области с шириной и длиной около 2 и 14 мм по полувысоте фокального пика в интервале глубин от 25 до 65 мм от внутренней поверхности черепа. Ширина и длина фокальной области практически совпадают с таковыми при фокусировке в воде. Таким образом, продемонстрирована возможность путем механического перемещения решетки в осевом направлении изменять глубину воздействия на структуры мозга в пределах 4 см. С уменьшением глубины фокусировки качество коррекции аберраций несколько ухудшается. Например, в случае 256-элементной решетки, для наименее глубокого положения фокуса (глубина 25 мм) амплитуда давления в фокусе составляет 91% от случая самой глубокой фокусировки (глубина 65 мм). Уровень максимумов давления вблизи черепа возрастает в 1.68 раз, в основном, за счет приближения костей черепа к высокоинтенсивной фокальной области в сходящемся ультразвуковом пучке. Этот эффект не связан с наличием черепа на пути пучка. Вспомогательные расчеты показывают, что если убрать из акустической модели все неоднородности и оставить в качестве среды только воду, то вычисления дадут такие же уровни давления в области предполагаемого расположения черепа.

Фокусировка в указанном выше интервале глубин была проанализирована для различного количества элементов решетки: 256; 512; 1024 и идеализированного случая квазинепрерывного изменения фазы по ее поверхности. Для всех решеток при уменьшении глубины фокусировки наблюдается одинаковый характер снижения уровня давления в фокусе: он уменьшается примерно на одинаковый процент при изменении глубины. При этом уровень максимумов давления вблизи черепа практически не зависит от количества излучающих элементов, что вполне соотносится с высказанным выше предположением о том, что данный эффект действительно определяется лишь геометрией фокусировки. При этом увеличение количества элементов позволяет повысить амплитуду давления в фокусе излучателя. Наибольший рост на 8% (на 0.67 дБ) наблюдается при переходе от 256 элементов к 512. Увеличение количества элементов до величин выше 1024 представляется менее выигрышным, поскольку рост амплитуды давления при переходе от 1024 элементов к идеализированной решетке составляет 7%, что даже меньше, чем при переходе от 256 к 512 элементам.

При оценке влияния различных эффектов на качество фокусировки показано, что наибольший вклад в ослабление пучка вносят аберрации, уменьшая максимально достижимый уровень амплитуды давления на 7.4 дБ, поглощение вносит потери величиной 6.7 дБ, отражение – 2.1 дБ, а генерация сдвиговых волн – 2 дБ. Данный результат наглядно демонстрирует важность дальнейшего развития алгоритмов коррекции аберраций как способа борьбы с самым значительным источником ослабления уровня воздействия.

Несмотря на то, что все расчеты в работе проводились с использованием линейных акустических моделей, практический интерес представляет применение результатов коррекции аберраций не только для реализации теплового метода абляции с использованием гармонических волн, но и метода гистотрипсии с кипением, основанного на механической деструкции тканей в фокусе нелинейными ударными волнами. Такой подход допустим, поскольку ранее было показано, что фазы на элементах решетки, рассчитанные в линейном приближении, позволяют скомпенсировать аберрации и при учете нелинейных эффектов [12]. Кроме того, накапливающиеся с расстоянием нелинейные искажения в профиле волны и образование ударных фронтов, необходимые для создания механических разрушений, вблизи поверхности черепа развиты слабее, чем в фокусе, а значит, не будут обладать деструктивными возможностями даже при сравнительно высоких уровнях давления. Для количественного учета нелинейных эффектов при определении безопасной области воздействия данных решеток расчеты должны проводиться с использованием более сложных численных моделей [27, 28]. Отметим также, что в данной работе проанализированы возможности коррекции аберраций при механическом перемещении излучателя. Отдельными вопросами являются обобщение предложенной модели на случай электронного перемещения фокуса вдоль и поперек оси решетки и анализ влияния неоднородностей внутренней структуры костей черепа на искажение ультразвукового пучка, что представляет предмет дальнейших исследований.

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

Работа выполнена в рамках Междисциплинарной научно-образовательной школы Московского университета “Фотонные и квантовые технологии. Цифровая медицина” при поддержке гранта РФФИ 19-02-00035.

Список литературы

  1. Qiu W., Bouakaz A., Konofagou E., Zheng H. Ultrasound for the brain: A review of physical and engineering principles, and clinical applications // IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 2020. V. 68. № 1. P. 6–20.

  2. McDannold N., Clement G., Black P., Jolesz F., Hynynen K. Transcranial MRI-guided focused ultrasound surgery of brain tumors: Initial findings in three patients // Neurosurgery. 2010. V. 66. № 2. P. 323–332.

  3. Jeanmonod D., Werner B., Morel A., Michels L., Zadicario E., Schiff G., Martin E. Transcranial magnetic resonance imaging-guided focused ultrasound: Noninvasive central lateral thalamotomy for chronic neuropathic pain // Neurosurg. Focus. 2012. V. 32. № 1. P. E1.

  4. Elias W.J., Huss D., Voss T., Loomba J., Khaled M., Zadicario E., Frysinger R.C., Sperling S.A., Wylie S., Monteith S.J., Druzgalm J., Shahm B.B., Harrison M., Wintermark M. A pilot study of focused ultrasound thalamotomy for essential tremor // N. Engl. J. Med. 2013. V. 369. № 7. P. 640–648.

  5. Elias W.J. A randomized trial of focused ultrasound thalamotomy for essential tremor // N. Engl. J. Med. 2016. V. 375 № 8, P. 730–739.

  6. Monteith S., Medel R., Kassell N.F., Wintermark W., Eames M., Snell J., Zadicario E., Grinfeld J., Sheehan J.P., Elias W.J. Transcranial magnetic resonance-guided focused ultrasound surgery for trigeminal neuralgia: A cadaveric and laboratory feasibility study // J. Neurosurg. 2013. V. 118. № 2. P. 319–328.

  7. Гаврилов Л.Р. Фокусированный ультразвук высокой интенсивности в медицине. М.: Фазис, 2013.

  8. Hynynen K., Jones R.M. Image-guided ultrasound phased arrays are a disruptive technology for non-invasive therapy // Phys. Med. Biol. 2016. V. 61. P. 206–248.

  9. Sun J., Hynynen K. Focusing of therapeutic ultrasound through a human skull: A numerical study // J. Acoust. Soc. Am. 1998. V. 104(3). P. 1705–1715.

  10. Pernot M., Aubry J.-F., Tanter M., Thomas J.-L., Fink M. High power transcranial beam steering for ultrasonic brain therapy // Phys. Med. Biol. 2003. V. 48. № 16. P. 2577–2589.

  11. Clement G.T., Hynynen K. A non-invasive method for focusing ultrasound through the human skull // Phys. Med. Biol. 2002. V. 47. № 8. P. 1219–1236.

  12. Rosnitskiy P.B., Yuldashev P.V., Sapozhnikov O.A., Gavrilov L.R., Khokhlova V.A. Simulation of nonlinear trans-skull focusing and formation of shocks in brain using a fully populated ultrasound array with aberration correction // J. Acoust. Soc. Am. 2019. V. 146. № 3. P. 1786–1798.

  13. Росницкий П.Б., Гаврилов Л.Р., Юлдашев П.В., Сапожников О.А., Хохлова В.А. О возможности применения многоэлементных фазированных решеток для ударно-волнового воздействия на глубокие структуры мозга // Акуст. журн. 2017. Т. 63. № 5. С. 489–500.

  14. Rosnitskiy P.B., Vysokanov B.A., Gavrilov L.R., Sapozhnikov O.A., Khokhlova V.A. Method for designing multielement fully populated random phased arrays for ultrasound surgery applications // IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 2018. V. 65. № 4. P. 630–637.

  15. Raju B.I., Hall C.S., Seip R. Ultrasound therapy transducers with space-filling non-periodic arrays // IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 2011. V. 58. № 5. P. 944–954.

  16. Ramaekers P., Ries M., Moonen C.T.W., de Greef M. Improved intercostal HIFU ablation using a phased array transducer based on Fermat’s spiral and Voronoi tessellation: A numerical evaluation // Med. Phys. 2017. V. 44. № 3. P. 1071–1088.

  17. Khokhlova T.D., Canney M.S., Khokhlova V.A., Sapozhnikov O.A., Crum L.A., Bailey M.R. Controlled tissue emulsification produced by high intensity focused ultrasound shock waves and millisecond boiling // J. Acoust. Soc. Am. 2011. V. 130. № 5. P. 3498–3510.

  18. Khokhlova V.A., Fowlkes J.B., Roberts W.W., Schade G.R., Xu Z., Khokhlova T.D., Hall T.L., Maxwell A.D., Wang Y.N., Cain C.A. Histotripsy methods in mechanical disintegration of tissue: Towards clinical applications // Int. J. Hyperthermia. 2015. V. 31. № 2. P. 145–162.

  19. Maxwell A., Sapozhnikov O., Bailey M., Crum L., Xu Z., Fowlkes B., Cain C., Khokhlova V. Disintegration of tissue using High Intensity Focused Ultrasound: two approaches that utilize shock waves // Acoustics Today. 2012. V. 8. P. 24–37.

  20. Ильин С.А., Юлдашев П.В., Хохлова В.А., Гаврилов Л.Р., Росницкий П.Б., Сапожников О.А. Применение аналитического метода для оценки качества акустических полей при электронном перемещении фокуса многоэлементных терапевтических решеток // Акуст. журн. 2015. Т. 61. № 1. С. 57–64.

  21. Bobkova S., Gavrilov L., Khokhlova V., Shaw A., Hand J. Focusing of high intensity ultrasound through the rib cage using therapeutic random phased array // Ultrasound Med. Biol. 2010. V. 36. № 6. P. 888–906.

  22. Treeby B.E., Cox B.T. Modeling power law absorption and dispersion in viscoelastic solids using a split-field and the fractional Laplacian // J. Acoust. Soc. Am. 2014. V. 136. № 4. P. 1499–1510.

  23. Treeby B.E., Jaros J., Rohrbach D., Cox B.T. Modelling elastic wave propagation using the k-Wave Matlab toolbox // IEEE Int. Ultrasonics Symp. 2014. P. 146–149.

  24. Duck F.A. Physical Properties of Tissue: A Comprehensive Reference Book // Academic Press, London, 1990.

  25. White P.J., Clement G.T., Hynynen K. Longitudinal and shear mode ultrasound propagation in human skull bone // Ultrasound Med. Biol. 2006.V. 32. № 7. P. 1085–1096.

  26. O'Neil H.T. Theory of Focusing Radiators // J. Acoust. Soc. Am. 1949. V. 21. № 5. P. 516–526.

  27. Юлдашев П.В., Xoxлова В.А. Моделирование трехмерных нелинейных полей ультразвуковых терапевтических решеток // Акуст. журн. 2011. Т. 57. № 3. С. 337–347.

  28. Rosnitskiy P.B., Yuldashev P.V., Sapozhnikov O.A., Maxwell A.D., Kreider W., Bailey M.R., Khokhlova V.A. Design of HIFU transducers for generating specified nonlinear ultrasound fields // IEEE Trans. Ultrason. Ferroelectr. Freq. Contr. 2017. V. 64. № 2. P. 374–390.

Дополнительные материалы отсутствуют.