Акустический журнал, 2021, T. 67, № 3, стр. 250-259

Влияние траектории перемещения фокуса на равномерность температурного поля при импульсном воздействии мощного ультразвукового пучка на биологическую ткань

П. А. Пестова a*, М. М. Карзова a, П. В. Юлдашев a, У. Крайдер b, В. А. Хохлова ab

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

b Университет штата Вашингтон, Лаборатория прикладной физики
1013 NE 40th Street, Сиэтл, США

* E-mail: ppolina-98@yandex.ru

Поступила в редакцию 17.07.2020
После доработки 04.02.2021
Принята к публикации 08.02.2021

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

Аннотация

Представлены результаты численного эксперимента по облучению образца ткани говяжьей печени ex vivo мощным фокусированным ультразвуком с помощью терапевтической решетки клинической системы Sonalleve V2 3.0T (Profound Medical Corp., Canada) в режимах, характерных для метода гистотрипсии с кипением. Облучение проводилось по траекториям, состоящим из дискретных фокусов, расположенных на двух либо четырех концентрических окружностях радиусами от 2 до 8 мм. Сравнивался режим облучения каждого дискретного фокуса последовательностью из пятнадцати импульсов и переходом к следующему фокусу с режимом, при котором происходит поочередное облучение всех дискретных фокусов на окружностях по одному разу за пятнадцать обходов по всему набору облучаемых точек. Проанализировано влияние размера траектории и количества импульсов на степень проявления тепловых эффектов. Фокусировка ультразвукового пучка в ткани описывалась с помощью уравнения Вестервельта, температурное поле моделировалось с помощью уравнения теплопроводности. Показано, что более однородная структура температурного поля обеспечивается при большем временном интервале между повторным облучением каждого из фокусов, а оптимальное значение временного интервала лежит в диапазоне от трех до шести периодов следования импульсов.

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

ВВЕДЕНИЕ

В последнее время технология неинвазивной ультразвуковой хирургии HIFU (High Intensity Focused Ultrasound) активно развивается как с целью совершенствования уже существующих подходов, так и разработки новых режимов облучения [13]. Разрушения биологической ткани, полученные при воздействии мощным фокусированным ультразвуком, могут быть двух видов – механические и тепловые [4]. Одним из методов получения механических разрушений ткани является метод гистотрипсии с кипением, в котором ткань облучается последовательностью импульсов миллисекундной длительности с коэффициентом заполнения около 1%, при этом в профилях давления в фокальной области ультразвукового пучка происходит образование высокоамплитудных ударных фронтов [4, 5]. В режимах гистотрипсии с кипением важно отсутствие тепловой денатурации ткани, однако на практике ее тепловой нагрев неизбежен за счет поглощения энергии ультразвуковой волны на ударных фронтах. Степень проявления тепловых эффектов необходимо уметь предсказывать при выборе протокола облучения, позволяющего минимизировать эффекты перегрева ткани для получения чисто механических разрушений.

В недавних экспериментах с использованием клинической системы MRg HIFU Sonalleve V1 3.0T (Profound Medical Corp., Canada) по получению объемных разрушений в ткани говяжьей печени ex vivo методом гистотрипсии с кипением было показано, что при длительности импульсов менее 10 мс и коэффициенте заполнения 1% возможно добиться чисто механического разрушения ткани без проявления эффектов тепловой денатурации [6]. В то же время было показано, что полученные разрушения сопровождались повышением температуры в объеме образца на 15–20°C. Такой нагрев не приводил к тепловой денатурации содержимого области разрушения, однако выявил важность разработки траекторий облучения, при которых достигается наиболее равномерный нагрев объема образца с наименьшими пиковыми температурами.

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

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

Геометрия численного эксперимента схематично показана на рис. 1а. В качестве излучателя рассматривалась обновленная версия V2 мощной фазированной решетки клинической системы MRg HIFU Sonalleve 3.0T (Profound Medical Corp., Canada). Рабочая поверхность решетки состоит из 256 элементов круглой формы с диаметром d = 6.6 мм и рабочей частотой   f0 = 1.2 МГц. Элементы решетки расположены на сегменте сферической поверхности с диаметром D = 136 мм и фокусным расстоянием F = 140 мм, разделенной на 8 секторов [7]. В работе был выбран ударно-волновой режим облучения, соответствующий акустической мощности решетки 645 Вт, при этом номинальное значение интенсивности на элементах решетки, рассчитанное как отношение мощности решетки к суммарной площади ее элементов, составило I0 = 7.4 Вт/см2.

Рис. 1.

(а) – Схема численного эксперимента. Ультразвуковой пучок создается HIFU-решеткой, состоящей из 256 элементов диаметром 6.6 мм с рабочей частотой 1.2 МГц; диаметр решетки 136 мм; фокусное расстояние F = 140 мм; фокусировка происходит на глубину h = 2 см в образец ткани говяжьей печени; излучатель и образец ткани помещены в воду. (б) – Две траектории дискретного перемещения фокуса ультразвукового пучка: первая траектория “4 мм” состоит из 2-х концентрических окружностей с радиусами 2 и 4 мм (последовательность электронного перемещения фокуса решетки отмечена числами); вторая траектория “8 мм” включает в себя 4 окружности с радиусами 2, 4, 6, 8 мм.

Ультразвуковой пучок, проходящий через согласующую среду (воду), фокусировался на глубине h = 2 см в образце ткани печени (рис. 1а). Фокус решетки перемещался в ее фокальной плоскости по траекториям, состоящим из набора дискретных фокусов, расположенных на двух либо четырех концентрических окружностях радиусами 2, 4, 6 и 8 мм, соответственно (рис. 1б). При этом облучение единичных фокусов на каждой из окружностей происходило в такой последовательности, чтобы поочередно облучаемые фокусы были по возможности максимально удалены друг от друга, а каждый из них не был расположен в соседней точке с фокусом, облучаемым перед предыдущим (показано с помощью нумерации на рис. 1б). Подобные траектории используются в клинической практике для получения объемных тепловых разрушений, а также соответствуют недавним экспериментам по разрушению биологических тканей методом гистотрипсии с кипением [6, 8].

Для траектории, состоящей из двух окружностей с радиусами 2 и 4 мм, рассматривалось два протокола облучения. В первом протоколе каждый дискретный фокус на окружностях траектории последовательно облучался пятнадцатью импульсами подряд, а затем происходило перемещение фокуса в следующую точку траектории. Во втором протоколе последовательно облучалась каждая точка траектории по одному разу, а всего совершалось пятнадцать обходов по двум окружностям. Поскольку в обоих протоколах производится одинаковое число единичных разрушений (360 штук), то время облучения для них будет одинаковым. Для траектории, состоящей из четырех окружностей с радиусами 2, 4, 6 и 8 мм, рассматривался только второй протокол облучения, при этом количество обходов по всем четырем окружностям было увеличено с пятнадцати до тридцати. Для анализа влияния размера траектории на характер температурного распределения аналогичный протокол облучения с тридцатью обходами всех окружностей был рассчитан для траектории меньшего размера, состоящей из двух окружностей с радиусами 2 и 4 мм.

Время облучения единичного фокуса theat в образце печени подбиралось таким образом, чтобы при нагреве от начальной температуры T0 = 23°C температура кипения (100°C) не достигалась. Затем, для нахождения интервала времени между началом последовательных облучений единичных фокусов, найденное значение theat умножалось на сто, что обеспечивало характерный для метода гистотрипсии с кипением коэффициент заполнения (1%).

ЧИСЛЕННАЯ МОДЕЛЬ

Акустическое поле и тепловые источники

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

(1)
$\frac{{{{\partial }^{2}}p}}{{\partial \tau \partial z}} = \frac{{{{c}_{0}}}}{2}\Delta p + \frac{\beta }{{2{{\rho }_{0}}c_{0}^{3}}}\frac{{{{\partial }^{2}}{{p}^{2}}}}{{\partial {{\tau }^{2}}}} + \frac{\delta }{{2c_{0}^{3}}}\frac{{{{\partial }^{3}}p}}{{\partial {{\tau }^{3}}}} + L(p),$
где p = p(x, y, z, τ) – давление, $\Delta = {{{{\partial }^{2}}} \mathord{\left/ {\vphantom {{{{\partial }^{2}}} {\partial {{x}^{2}}}}} \right. \kern-0em} {\partial {{x}^{2}}}}$ + $ + \,\,{{{{\partial }^{2}}} \mathord{\left/ {\vphantom {{{{\partial }^{2}}} {\partial {{y}^{2}} + {{{{\partial }^{2}}} \mathord{\left/ {\vphantom {{{{\partial }^{2}}} {\partial {{z}^{2}}}}} \right. \kern-0em} {\partial {{z}^{2}}}}}}} \right. \kern-0em} {\partial {{y}^{2}} + {{{{\partial }^{2}}} \mathord{\left/ {\vphantom {{{{\partial }^{2}}} {\partial {{z}^{2}}}}} \right. \kern-0em} {\partial {{z}^{2}}}}}}$ – оператор Лапласа, z – координата, вдоль которой происходит фокусировка пучка, τ = tz/c0 – время в сопровождающей системе координат, параметры c0, β, ρ0 и δ – скорость звука, коэффициент нелинейности, плотность среды и коэффициент термовязкого поглощения в среде, соответственно. Значения указанных физических параметров для воды были равны ρ0 = 997 кг/м3, c0 = 1485 м/с, β = 3.5, а для ткани говяжьей печени ρ0 = 1050 кг/м3, c0 = 1580 м/с, β = 4.0; коэффициент термовязкого поглощения в обеих средах был выбран одинаковым и составлял δ = 4.33 × 10–6 м2/с [9, 10].

В дополнение к термовязкому поглощению, в ткани печени также учитывались частотная зависимость коэффициента поглощения и соответствующая дисперсия, характерные для биологических тканей. В уравнении (1), записанном во временном представлении, эти эффекты описаны неявно с помощью оператора L(p). В частотном представлении оператора L(p) коэффициент поглощения в ткани печени α соответствовал линейной зависимости от частоты: ${{\alpha }}(f) = {{{{\alpha }}}_{0}}\frac{f}{{{{f}_{0}}}},$ где α0 = 7.2 м–1 – коэффициент поглощения на частоте f0 = 1.2 МГц, а дисперсия скорости звука c рассчитывалась с помощью локальных дисперсионных соотношений: $\frac{{c(f) - {{c}_{0}}}}{{{{c}_{0}}}} = \frac{{{{c}_{0}}{{\alpha }_{0}}}}{{{{\pi }^{2}}{{f}_{0}}}}\ln \left( {\frac{f}{{{{f}_{0}}}}} \right)$ [9, 11]. Отметим, что при указанных физических параметрах и достижимом уровне давления в ткани печени в фокусе решетки порядка ${{p}_{F}}$ = 90 МПа, значение обратного акустического числа Рейнольдса, ${\text{R}}{{{\text{e}}}^{{ - 1}}} = {{{{\pi \delta }}{{{{\rho }}}_{0}}{{f}_{0}}} \mathord{\left/ {\vphantom {{{{\pi \delta }}{{{{\rho }}}_{0}}{{f}_{0}}} {{{\beta }}{{p}_{F}}}}} \right. \kern-0em} {{{\beta }}{{p}_{F}}}}$, определяющего ширину ударного фронта относительно периода волны, составляет около 0.5 × 10–4. Отношение длины образования разрыва, ${{z}_{р}} = {{c_{0}^{3}{{{{\rho }}}_{0}}} \mathord{\left/ {\vphantom {{c_{0}^{3}{{{{\rho }}}_{0}}} {2{{\pi \beta }}{{p}_{F}}{{f}_{0}}}}} \right. \kern-0em} {2{{\pi \beta }}{{p}_{F}}{{f}_{0}}}}$= 1.5 мм, к длине поглощения на основной частоте, ${{z}_{п}} = {{\alpha }}_{0}^{{ - 1}}$ = = 140 мм, также мало, $~{{{{z}_{р}}} \mathord{\left/ {\vphantom {{{{z}_{р}}} {{{z}_{п}}}}} \right. \kern-0em} {{{z}_{п}}}}$ = 10–2, что говорит о существенном доминировании нелинейных эффектов над диссипативными процессами.

Для задания граничного условия в уравнении Вестервельта (1), соответствующего терапевтической решетке Sonalleve V2, использовался метод акустической голографии [12]. Метод основан на измерении двумерных распределений амплитуды и фазы акустического давления в плоскости, расположенной между излучателем и его фокальной плоскостью, а затем расчете обратного распространения акустического поля на поверхность источника. Акустическая голограмма реальной решетки была измерена ранее в дегазированной воде в плоскости (x, y, z = 99.8 мм) [7]. Измерения профилей давления проводились с помощью капсульного гидрофона HGL-0085, соединенного с предусилителем AH-2020 (Onda Corp., Sunnyvale, CA), и соответствовали мощности решетки 42 Вт. Положение гидрофона задавалось 3D системой позиционирования (Velmex NF90, Bloomfield, NY), сканирование поля производилось с шагом 0.6 мм по поперечным координатам x и y в пространственном окне размером 8 × 8 см. В каждой точке записывалось 128 периодов гармонической волны, по которым затем проводилось усреднение с целью уменьшения шумов и получения результирующего сигнала. Измеренные двумерные распределения акустического поля пересчитывались обратно до плоскости (x, y, z = 0 мм) с помощью метода углового спектра. Для обеспечения рассматриваемой в данной работе мощности решетки (645 Вт) полученное распределение амплитуды давления в плоскости (x, y, z = 0 мм) домножалось на соответствующий коэффициент (3.9), равный квадратному корню из соотношения требуемой и использованной в голографических измерениях мощностей. Полученные распределения амплитуды и фазы акустического давления использовались в качестве граничного условия в моделировании (рис. 2а, 2б). Методика похожего эксперимента и постановки граничного условия для решетки предыдущей версии системы Sonalleve (V1) подробно описаны в работе [12].

Рис. 2.

Граничное условие в плоскости z = 0 мм для моделирования ультразвукового поля решетки: двумерные распределения (а) – амплитуды и (б) – фазы акустического давления.

Для моделирования уравнения Вестервельта (1) использовался численный алгоритм, развитый в работе [13], а затем многократно апробированный нашей группой в расчетах мощных ультразвуковых полей различных медицинских излучателей [9, 12, 14, 15]. Алгоритм основан на методе расщепления по физическим факторам с использованием конечно-разностных методов и преимуществ как спектрального, так и временного представлений акустического поля.

Результаты моделирования уравнения (1) использовались для нахождения пространственного распределения плотности мощности тепловых источников Q(x, y, z) в ткани печени. Для этого профиль давления в каждой точке пространства сначала представлялся в виде разложения в ряд Фурье, затем в этой точке находилась полная интенсивность волны I как сумма интенсивностей всех гармоник (N = 800) с комплексными амплитудами давления. Мощность тепловыделения Q в среде за счет поглощения энергии волны считалась как скорость убыли интенсивности на каждом шаге сетки dz [16]:

(2)
$Q(x,y,z) = - \frac{{I(x,y,z + dz) - I(x,y,z)}}{{dz}}.$

Тепловые источники рассчитывались в узлах пространственной сетки с поперечными шагами dx = dy = 0.02 мм, продольным шагом dz = 0.1 мм и выводились в окне размерами [–16, 16] мм по осям x и y и [120, 160] мм по оси z.

Температурное поле

Температурное поле в образце ткани печени описывалось с помощью неоднородного уравнения теплопроводности:

(3)
$\frac{{\partial T}}{{\partial t}} = \chi \Delta T + \frac{Q}{{{{C}_{v}}}},$
где T = T(x, y, z, t) – температура в ткани, t – время, χ – коэффициент температуропроводности, Cv – теплоемкость единицы объема образца, Q(xy, z) – полученная на основе решения уравнения Вестервельта плотность мощности тепловых источников в ткани (2). Значения физических параметров в уравнении (3) соответствовали ткани печени и были равны χ = 1.93 × 10–7 м2/с, Cv = 5 × 106 Дж/(м3 °C) [9, 10, 17].

Алгоритм решения уравнения (3) подробно описан в работе и включал в себя использование аналитического решения уравнения теплопроводности в k-пространстве:

(4)
$\begin{gathered} \hat {T}({\mathbf{k}},t) = {{{\hat {T}}}_{0}}({\mathbf{k}})\exp \left( { - 4{{\pi }^{2}}{{{\mathbf{k}}}^{2}}\chi t} \right) + \\ + \,\,\frac{{\hat {Q}({\mathbf{k}})}}{{4{{\pi }^{2}}{{{\mathbf{k}}}^{2}}\chi {{C}_{v}}}}\left[ {1 - \exp \left( { - 4{{\pi }^{2}}{{{\mathbf{k}}}^{2}}\chi t} \right)} \right], \\ \end{gathered} $
где $\hat {T}({\mathbf{k}},t),\;{{\hat {T}}_{0}}({\mathbf{k}}),\;\hat {Q}({\mathbf{k}})$ – пространственные Фурье-спектры от соответствующих величин $T(x,y,z,t),$ ${{T}_{0}}(x,y,z),Q(x,y,z)$, а ${{T}_{0}}(x,y,z)$ – начальное распределение температуры в рассматриваемом объеме. Переходы между декартовыми координатами и k-пространством осуществлялись с помощью операций быстрого преобразования Фурье (БПФ), входящих в стандартную библиотеку FFTW.

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

Тестовые расчеты показали, что для корректного описания температурного поля в объеме ткани, решение для единичного воздействия должно быть взято в таком пространственном окне, чтобы при облучении по траектории добавляемое решение перекрывалось с самой удаленной точкой траектории. Таким образом, пространственное окно температурного распределения единичного фокуса должно быть, по крайней мере, в два раза больше диаметра наибольшей окружности траектории. Поэтому решение для температурного поля в единичном фокусе в случае траектории с внешним радиусом 4 мм бралось в поперечном окне [–8, 8] мм, а с радиусом 8 мм – [–16, 16] мм.

При расчетах объемного нагрева, чтобы избежать эффекта наложения частот при взятии операций БПФ, размеры пространственных окон подбирались таким образом, чтобы на границе окна не происходило повышения температуры. Расчет температуры при движении фокуса по траектории, состоящей из двух окружностей с радиусами 2 и 4 мм, проводился с шагами dx = dy = 0.02 мм и dz = 0.1 мм в окне размерами [–16, 16] мм по осям x и y и [120, 160] мм по оси z. Для траектории, состоящей из четырех окружностей, пространственные шаги были увеличены в два раза, dx = dy = 0.04 мм и dz = 0.02 мм, а размеры окна расчета составили [‒26, 26] мм по осям x и y и [100, 180] мм по оси z.

РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ

Акустическое поле и тепловые источники

Из результатов моделирования следует, что при начальной интенсивности на элементах терапевтической решетки равной 7.4 Вт/см2 (полная мощность 645 Вт) облучение образца печени происходит в ударно-волновом режиме. В профиле волны в фокусе формируется развитый разрыв с амплитудой ударного фронта 88 МПа (рис. 3а), при этом фокус решетки, определяемый как положение максимума плотности мощности тепловых источников, смещается из геометрического фокуса (z = F = 140.0 мм) в направлении к излучателю (z =138.4 мм).

Рис. 3.

Результаты численного моделирования акустического поля с помощью уравнения Вестервельта. (а) – Профиль давления с развитым разрывом Ap = 88 МПа в фокусе (z =138.4 мм) при номинальной интенсивности на элементах решетки 7.4 Вт/см2. (б, в) – Пространственные распределения плотности мощности тепловых источников в ткани в (б) – фокальной и (в) – аксиальной плоскостях. Местоположение фокуса здесь и далее соответствует положению максимума плотности мощности тепловых источников.

На распределениях плотности мощности тепловых источников Q в ткани печени (рис. 3б, 3в) хорошо видно, что поглощение энергии акустической волны, обусловленное присутствием ударных фронтов в ее профиле, локализовано в очень малом объеме: поперечные размеры теплового пятна, рассчитанные по уровню e–1 от пикового значения, составляют всего 0.24 × 0.18 мм (длины осей эллипса на рис. 3б), а продольный размер составляет 3.8 мм (рис. 3в). Форма фокальной области в виде эллипса вместо окружности объясняется тем, что экспериментально полученные распределения амплитуды и фазы акустического давления реальной решетки Sonalleve V2 3.0T не идеально симметричны (рис. 2). Пиковое значение плотности мощности тепловых источников Qmax =100 Вт/мм3 позволяет говорить о возможности локального достижения температуры кипения в ткани за миллисекунды. Оценить длительность импульсов tboil, нагревающих ткань до температуры кипения, без учета диффузии тепла, можно по формуле [18, 19]:

(5)
${{t}_{{{\text{boil}}}}} = \frac{{\Delta T{{C}_{v}}}}{{{{Q}_{{{\text{max}}}}}}},$
где ΔΤ – изменение температуры от начальной T0 = 23°C до температуры кипения (100°C). Результат оценки составляет tboil = 3.8 мс, что является характерным временем для режима гистотрипсии с кипением [18].

Температурное поле единичного воздействия

Численный расчет уравнения теплопроводности (3) с использованием решения (4) показал, что для единичного облучения образца печени от начальной температуры T0 = 23°C температура, близкая к температуре кипения (98.6°C) ткани, достигается за время theat = 4.1 мс. Незначительное превышение полученного значения theat над оценочным значением tboil = 3.8 мс (5) обусловлено малым размером области эффективного тепловыделения и диффузией тепла в окружающий объем в процессе нагрева. Для обеспечения коэффициента заполнения, характерного для режима гистотрипсии с кипением, время между началом последовательных облучений единичных фокусов было выбрано как tsingle = 100 × theat = 0.41 с.

В момент окончания нагрева t = theat= 4.1 мс пространственное распределение температуры единичного фокуса по форме и размерам полностью повторяет пространственное распределение плотности мощности тепловых источников Q, показанное на рис. 3б, 3в. К моменту начала повторного воздействия следующего импульса t = tsingle = 0.41 с тепловое пятно расплывается за счет эффектов диффузии тепла (рис. 4а). На рис. 4б показаны одномерные распределения температурного поля единичного воздействия в эти моменты времени: сплошная кривая соответствует времени окончания нагрева t = theat= 4.1 мс (температурная шкала слева), а пунктирная линия – моменту переключения на следующий фокус t = tsingle = 0.41 с (температурная шкала справа). К моменту повторного облучения образца печени t = tsingle тепловое пятно единичного воздействия успевает остыть с 98.6°C до 27.9°C (рис. 4б), при этом размеры теплового пятна увеличиваются от начальных размеров области тепловых источников 0.24 × 0.18 × 3.8 мм до размеров 1.3 × 1.3 × 5.3 мм, а форма пятна в фокальной плоскости меняется с эллиптической на аксиально симметричную. Таким образом, повышение температуры в образце печени в текущем единичном фокусе к моменту начала облучения следующего фокуса составляет около 5°C.

Рис. 4.

Результаты численного моделирования температурного поля при единичном воздействии. (а) – Пространственные распределения температуры единичного фокуса в момент времени t = 0.41 c в фокальной (слева) и аксиальной (справа) плоскостях. (б) – Распределения температуры в моменты времени t = 4.1 мc (сплошная линия, температурная шкала слева) и t = 0.41 с (пунктир, температурная шкала справа) в фокальной плоскости (слева) и вдоль оси пучка (справа).

Анализ дальнейшего расплывания теплового пятна единичного воздействия на временах t > tsingle показал, что остаточный нагрев становится меньше 1°C при t = 3tsingle (максимальный остаточный нагрев ткани 0.95°C), при этом размеры теплового пятна по уровню e–1 от пикового значения составляют 2.96 × 2.96 × 7.6 мм по осям x, y и z, соответственно. Однородное распределение температуры, при котором пиковый нагрев не превышает 0.3°C, достигается при t = 6tsingle (максимальный остаточный нагрев ткани 0.27°C), а тепловое пятно увеличивается до размеров 5.52 × 5.52 × 12.5 мм.

Облучение по траектории из дискретных фокусов

Для получения клинически значимого объема разрушения в ткани в методе гистотрипсии с кипением используются протоколы, в которых фокус излучателя перемещается по различным траекториям. Сравнение результатов облучения для двух характерных для метода гистотрипсии с кипением протоколов показало, что первый протокол с последовательным облучением каждого фокуса пятнадцатью импульсами подряд формирует неоднородное и несимметричное относительно центра объема ткани температурное пятно (рис. 5а, сверху). Пиковое значение температуры в этом режиме превышает начальную температуру образца на 24.8°C в момент окончания облучения t = 360tsingle. Таким образом, накопление тепла в одном фокусе при его последовательном облучении несколькими импульсами подряд может приводить к значительному проявлению тепловых эффектов.

Рис. 5.

Пространственные распределения изменения температуры в момент окончания облучения объема ткани в (а) – фокальной и (б) – аксиальной плоскостях для двух различных протоколов облучения: протокол 1 (сверху) с последовательным облучением каждой точки траектории пятнадцатью импульсами подряд и перемещением фокуса в следующую точку траектории; протокол 2 (снизу) с последовательным облучением каждой точки траектории по одному разу и пятнадцатью полными обходами по траектории. (в) – Зависимости изменения температуры от времени в центре объема ткани для протокола 1 (штрих-пунктирная линия) и протокола 2 (сплошная линия).

Во втором протоколе облучения, при котором каждая точка траектории облучалась один раз, но совершалось пятнадцать обходов по траектории, каждый фокус к моменту повторного облучения успевал остыть, в результате чего нагрев был более однородным, симметричным, и достигалось меньшее пиковое значение (13.4°C) увеличения температуры (рис. 5а, снизу). В этом режиме на температурных распределениях визуально остаются хорошо различимы локализации пяти последних единичных облучений, в то время как температурные пятна от более ранних воздействий успевают расплыться по объему образца (рис. 5а, снизу). Такая особенность пространственного распределения согласуется с полученным выше критерием на временной интервал 6tsingle, в течение которого нагрев от единичного воздействия успевает стать однородным.

В аксиальной плоскости различие температурных распределений менее выражено (рис. 5б), поскольку тепло от удаленных от оси излучателя точек траектории не успевает распространиться к указанной плоскости. Продольные размеры температурных пятен в аксиальной плоскости в двух режимах составляют порядка 20–25 мм по уровню e–1 от пикового значения.

Анализ изменения температуры в центральной точке образца, которая непосредственно не облучалась, но нагревалась за счет эффектов диффузии тепла, показал, что в первом протоколе облучения сначала происходит быстрый рост температуры, затем достигается максимум (ΔΤ  = 11°C), а после перехода по траектории облучения с внутренней окружности на внешнюю центр образца начинает остывать и, начиная с середины времени облучения, температура изменяется незначительно (рис. 5в). Полное время облучения для каждого протокола было одинаковым и составляло t = 360tsingle = 147.6 с. Облучение с использованием второго протокола, в котором временной интервал между повторным облучением одной и той же точки фокуса на 9.43 с больше, чем в первом протоколе, обеспечивает постепенное нагревание центра образца. На графике зависимости температуры центра от времени присутствуют периодические осцилляции с амплитудой порядка 0.5°C и периодом t = 24tsingle = 9.84 с, равным времени одного полного обхода по траектории из двух окружностей. Такой характер зависимости объясняется более эффективным нагревом центра от точек, расположенных на внутренней окружности траектории, и менее эффективным от точек с внешней окружности. Поскольку во втором протоколе в траектории облучения происходит последовательная смена фокусировки с внутренней окружности на внешнюю и обратно, то возникают упомянутые осцилляции.

Для исследования влияния размера траектории на равномерность и однородность конечного распределения температуры объемного разрушения был проведен численный эксперимент, в котором облучение происходило по второму протоколу для двух траекторий с внешними радиусами 4 и 8 мм. Первая траектория состояла из двух окружностей с радиусами 2 и 4 мм, вторая – из четырех с радиусами 2, 4, 6 и 8 мм. Каждая точка на обеих траекториях облучалась один раз за один обход по траектории, а всего было совершено тридцать полных обходов. Полное время облучения для случая большего объема составляло t = = 2400tsingle = 984 с, что в 3.3 раза дольше времени облучения меньшего объема (t = 720tsingle = 295.2 с), соответственно временной интервал между повторными воздействиями в одну и ту же точку был в такое же число раз больше (t = 80tsingle = 32.8 с для траектории с радиусом 8 мм и t = 24tsingle = 9.84 с для траектории с радиусом 4 мм).

Увеличение размера траектории, с одной стороны, приводит к большему интегральному количеству выделившегося в ткани тепла. С другой стороны, увеличивается облучаемая площадь, а значит, становится более эффективным процесс остывания образца. Поскольку количество тепла, выделившегося в ткани, увеличивается только в 3.3 раза с увеличением радиуса внешней окружности траектории с 4 до 8 мм, а эффективная площадь остывания увеличивается в 4 раза, то следует ожидать, что остаточный нагрев в случае облучения большего объема окажется меньше, чем для разрушения меньшего объема. Этот предварительный оценочный результат подтвердился в численном моделировании (рис. 6а, 6б): остаточный пиковый нагрев в конце облучения по траектории большего размера составил ΔΤ = 14.9°C, в то время как для траектории меньшего размера ΔΤ  = 16.4°C. На фокальных распределениях (рис. 6а) яркими круговыми пятнами прослеживаются локализации пяти последних фокусов каждой траектории. Интересно отметить, что с увеличением диаметра траектории облучения в два раза размер теплового пятна вдоль оси изменился не так сильно – с 31 до 38 мм по уровню e–1 от максимума (рис. 6б).

Рис. 6.

Пространственные распределения изменения температуры в момент окончания облучения объема ткани в (а) – фокальной и (б) – аксиальной плоскостях для траектории, состоящей из двух концентрических окружностей радиусами 2 и 4 мм (вверху) и траектории, состоящей из четырех концентрических окружностей радиусами 2, 4, 6 и 8 мм (внизу). (в) – Зависимости изменения температуры от времени в центре объема ткани для траектории радиусом 4 мм (сплошная линия) и 8 мм (штрих-пунктирная линия).

Анализ нагрева центра образца (рис. 6в) показал, что, несмотря на более длительное облучение и большее интегральное количество тепла в случае траектории радиусом 8 мм, ее центр также нагревается меньше (ΔΤ  =12.6°C в момент окончания протокола), чем центр траектории с радиусом 4 мм (ΔΤ = 13.8°C в момент окончания протокола). Осцилляции, возникающие на температурных кривых, по-прежнему обусловлены особенностями процесса теплопередачи: при облучении первых двенадцати фокусов траектории происходит эффективный нагрев центра, затем, при удалении облучаемых точек от центральной, последняя остывает. Период осцилляций в каждом случае равен времени облучения всех единичных фокусов траектории по одному разу.

Полученные результаты численного моделирования, несмотря на отсутствие учета вскипания ткани в численной модели, хорошо согласуются с известными экспериментальными данными [6], в которых нагрев центра образца ткани печени составлял 15–20°C.

ВЫВОДЫ

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

При рассмотренных в работе условиях облучения ткани, для получения однородного распределения температурного поля к моменту окончания облучения образца временной интервал между повторным воздействием на один и тот же фокус дискретной траектории должен составлять t = 6tsingle. При выборе такого критерия пиковый нагрев ткани от предыдущего единичного воздействия не будет превышать 0.3°C. Однако, если остаточный нагрев ткани в 1°C является незначительным по сравнению с фоновой температурой образца, то критерий можно смягчить и уменьшить время между облучением одной и той же точки траектории до t = 3tsingle. Эта оценка может меняться с изменением мощности решетки и глубины фокусировки в биологическую ткань.

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

ЗАКЛЮЧЕНИЕ

В работе показано, что при облучении биологической ткани в ударно-волновом режиме c целью получения объемных механических разрушений степень проявления тепловых эффектов сильно зависит от выбора траектории перемещения фокуса. При последовательном облучении каждой точки двух круговых концентрических траекторий температурное поле в ткани после пятнадцати полных обходов более симметрично и однородно, а пиковое значение температуры на 11°C ниже, чем в случае облучения каждого фокуса пятнадцатью импульсами подряд. Обеспечение большего временного интервала между повторным облучением каждого фокуса позволяет получить более однородную пространственную структуру температурного поля и уменьшить максимальную достигаемую температуру при получении объемных разрушений. При увеличении числа полных обходов по траектории, состоящей из двух окружностей, с пятнадцати до тридцати пиковый нагрев в момент окончания облучения увеличился с 13.4°C до 16.4°C. Увеличение размера траектории в два раза и включение в нее еще двух окружностей привело к уменьшению пикового нагрева и температуры центра образца в момент окончания протокола облучения примерно на 1°C.

Авторы выражают благодарность коллективу Лаборатории прикладной физики университета шт. Вашингтон (США) за предоставление экспериментальных данных по акустической голографии и Л.Р. Гаврилову за полезные обсуждения работы.

Работа выполнена при финансовой поддержке грантов РФФИ 20-32-70142 и 20-02-00210, Focused Ultrasound Foundation и студенческой стипендии фонда “БАЗИС”.

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

  1. Хилл К.Р., Бэмбер Дж. Ультразвук в медицине. Физические основы применения / Под ред. тер Хаар Г. Пер. с англ. М.: Физматлит, 2008.

  2. Бэйли М.Р., Хохлова В.А., Сапожников О.А., Каргл С.Г., Крам Л.А. Физические механизмы воздействия терапевтического ультразвука на биологическую ткань (обзор) // Акуст. журн. 2003. Т. 49. № 4. С. 437–464.

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

  4. 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.C. Histotripsy methods in mechanical desintegration of tissue: towards clinical applications // Int. J. Hyperthermia. 2015. V. 31. № 2. P. 145–162.

  5. Maxwell A.D., Sapozhnikov O.A., Bailey M.R., Crum L.A., Xu Z., Fowlkes B., Cain C., Khokhlova V.A. Disintegration of tissue using high intensity focused ultrasound: Two approaches that utilize shock waves // Acoust. Today. 2012. V. 8. № 4. P. 24–36.

  6. Khokhlova V.A., Maxwell A.D., Khokhlova T., Kreider W., Bailey M., Partanen A., Farr N., Sapozhnikov O. Generation of volumetric boiling histotripsy lesions in tissue using a multi-element array of a clinical HIFU system // Abstract Book of 14th Int. Symp. for Therapeutic Ultrasound. Las Vegas, Nevada, USA. 2014.

  7. Karzova M.M., Yuldashev P.V., Kreider W., Rosnitskiy P.B., Khokhlova T.D., Sapozhnikov O.A., Bawiec C., Partanen A., Khokhlova V.A. Comparison of Sonalleve V1 and V2 MR-HIFU systems for generating high-amplitude shock-wave fields // Proc. the 6th Int. Symp. on Focused Ultrasound. Reston, Virginia, USA. October 21–25, 2018.

  8. Köhler M.O., Mougenot C., Quesson B., Enholm J., Le Bail B., Laurent C., Moonen C.T.W., Ehnholm G.J. Volumetric HIFU ablation under 3D guidance of rapid MRI thermometry // Med. Phys. 2009. V. 36. № 8. P. 3521–3535.

  9. Yuldashev P., Shmeleva S., Ilyin S., Sapozhnikov O., Gavrilov L., Khokhlova V. The role of acoustic nonlinearity in tissue heating behind the rib cage using high intensity focused ultrasound phased array // Phys. Med. Biol. 2013. V. 58. № 8. P. 2537–2559.

  10. Duck F.A. Physical Properties of Tissue. London: Academic Press, 1990.

  11. Филоненко E.А., Хохлова В.А. Эффекты акустической нелинейности при терапевтическом воздействии мощного фокусированного ультразвука на биологическую ткань // Акуст. журн. 2001. Т. 47. № 4. С. 541–549.

  12. Kreider W., Yuldashev P.V., Sapozhnikov O.A., Farr N., Partanen A., Bailey M.R., Khokhlova V.A. Characterization of a multi-element clinical HIFU system using acoustic holography and nonlinear modeling // IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2003. V. 60. № 8. P. 1683–1698.

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

  14. Karzova M.M., Yuldashev P.V., Sapozhnikov O.A., Khokhlova V.A., Cunitz B.W., Kreider W., Bailey M.R. Shock formation and nonlinear saturation effects in the ultrasound field of a diagnostic curvilinear probe // J. Acoust. Soc. Am. 2017. V. 141. № 4. P. 2327–2337.

  15. Maxwell A.D., Yuldashev P.V., Kreider W., Khokhlova T.D., Schade G.R., Hall T.L., Sapozhnikov O.A., Bailey M.R., Khokhlova V.A. A prototype therapy system for transcutaneous application of boiling histotripsy // IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2017. V. 64. № 10. P. 1542–1557.

  16. Андрияхина Ю.С., Карзова М.М., Юлдашев П.В., Хохлова В.А. Ускорение тепловой абляции объемов биологической ткани с использованием фокусированных ультразвуковых пучков с ударными фронтами // Акуст. журн. 2019. Т. 65. № 2. С. 1–12.

  17. Karzova M.M., Kreider W., Partanen A., Sapozhnikov O.A., Khokhlova T.D., Yuldashev P.V., Khokhlova V.A. Mapping clinical HIFU thermal tissue ablation using simulation and MR-imaging // Abstract book of the 19th Int. Symp. of ISTU / 5th European Symp. of EUFUS. Barcelona, Spain, 2019. P. 235.

  18. Canney M., Khokhlova V., Bessonova O., Bailey M., Crum L. Shock-induced heating and millisecond boiling in gels and tissue due to high intensity focused ultrasound // Ultrasound Med. Biol. 2010. V. 36. P. 250–267.

  19. 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.

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