На правах рукописи

Чупикин Дмитрий Анатольевич

ДОЗИМЕТРИЧЕСКОЕ ПЛАНИРОВАНИЕ ДИСТАНЦИОННОЙ

ЛУЧЕВОЙ ТЕРАПИИ НА ОСНОВЕ МЕТОДА МОНТЕ-КАРЛО

Специальность

01.04.16 – Физика атомного ядра и элементарных частиц

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

Москва 2007 2

Работа выполнена в Государственном образовательном учреждении высшего профессионального образования «МИФИ» - Московском государственном инженерно-физическом институте (техническом университете)

Научный руководитель:

Московский государственный инженерно-физический институт, г. Москва Доктор физико-математических наук, Профессор Климанов Владимир Александрович

Научный консультант:

ВНИИ экспериментальной физики, г. Саров Кандидат физико-математических наук, доцент Донской Евгений Николаевич

Официальные оппоненты:

Обнинский Государственный Технический Университет Атомной Энергетики доктор физико-математических наук, профессор Андросенко Петр Александрович НИИ экспериментальной ядерной физики МГУ, г. Москва кандидат физико-математических наук, доцент Кэбин Эдуард Йоханнесович

Ведущая организация:

ГНЦ РФ Институт теоретической и экспериментальной физики РАН, г. Москва

Защита состоится « 01 » ноября 2007 г. в 14 часов на заседании диссертационного совета Д 501.001.65 при Биологическом факультете МГУ им. М. В. Ломоносова – Московском Государственном Университете им. М.В. Ломоносова по адресу: 119992, Россия, Москва, Ленинские горы, д. 1, стр. 12, Биологический факультет МГУ, ауд. _.

С диссертацией можно ознакомиться в библиотеке Биологического факультета МГУ им. М. В. Ломоносова Отзыв на автореферат в 2-х экз. просим направлять по указанному адресу.

Автореферат разослан «»2007 г.

Ученый секретарь диссертационного совета, кандидат биологических наук Веселова Т.В.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность работы.

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

В результате в большинстве клиник стали использоваться системы 3-х мерного дозиметрического планирования с соответствующими алгоритмами 3-х мерного расчета дозы.

В настоящее время в таких передовых системах применяются в основном два метода: 1) метод свертки/суперпозиции, использующий дозовые ядра дифференциального тонкого луча; 2) метод тонкого луча в различных модификациях. Эти методы пришли на смену эмпирическим и полуэмпирическим методам. Они обладают достаточной быстротой и обеспечивают в случае однородности среды высокую точность расчета. Однако при наличии гетерогенности в области расчета, а также вблизи поверхности тела пациента погрешность расчета с помощью этих методов сильно возрастает. Принято считать, что единственным методом, который в условиях сложной геометрии расчета может обеспечить рекомендуемую международными комиссиями точность, является вероятностный метод Монте-Карло. Важно, что метод не требует записи соответствующего уравнения переноса, фактически здесь нет проблем с формулировкой дискретной модели. В обоих случаях реализация вычислительного алгоритма не требует упрощений. Метод максимально приспособлен для использования на ЭВМ. Единственным ограничением являются затраты машинного времени для получения результата с необходимой точностью. В настоящее время эти проблемы не являются решенными. Но в последнее время и в этом направлении ведутся целенаправленные работы и имеются определенные достижения.

Таким образом, метод Монте-Карло является наиболее точным методом расчета доз в лучевой терапии, и поэтому, ускорение расчетов путем модификации метода Монте-Карло является актуальной научно-технической задачей.

Целью работы является усовершенствование одной из модификаций метода Монте-Карло, а именно, метода PL-оценок потока в точке, а также разработка новых моделей источников для использования их в методе Монте-Карло.

При этом решаются следующие основные задачи:

1.Поиск возможных способов ускорения расчетов поглощенной дозы с использованием метода PL-оценок. На основе исследований предлагаются методы ускорения расчетов.

2. Метод PL-оценок потока в точке обобщается на случай расчета доз в неоднородной среде.

3. Разработка корректных математических моделей источников излучения радиотерапевтических ускорителей и облучателей с Co с целью использования их в методе PL-оценок.

4. Разработка моделей фотонных источников излучения радиотерапевтических установок для использования их в расчетах различными алгоритмами метода Монте-Карло в лучевой терапии.

Теоретическая база исследования. Для решения поставленных в работе задач использовалась теория переноса излучения, теория вероятностей и математическая статистика.

Научная новизна настоящей работы заключается в следующем:

1. Впервые предложены методы ускорения метода PL-оценок потока в точке.

2.Впервые показана возможность использования метода PL-оценок потока в точке при расчетах доз в химически неоднородных средах.

3.Усовершенствованы модели источников используемых в методе PLоценок потока в точке.

4. Предложен оригинальный метод формирования функций распределения вероятности, пригодных для выборки при моделировании фотонных источников для расчетов поглощенных доз методом Монте-Карло.

5. Предложен метод восстановления энергетического спектра дистанционного радиотерапевтического облучателя на основе пространственного распределения энергопоглощений.

Практическая значимость результатов работы состоит прежде всего в том, что предложенные в работе быстрые модификации метода PL-оценок с корректным учетом неоднородностей, а также новые математические модели источников дистанционной лучевой терапии, позволили создать на основе метода Монте-Карло методику расчета доз в лучевой терапии за приемлемое для клиник расчетное время.

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

Основные положения, выносимые на защиту:

1. Способы ускорения расчетов при использовании метода PL-оценок.

2. Методика учета гетерогенностей при расчете доз в лучевой терапии при использовании метода PL-оценок.

3. Новые модели фотонных источников излучения для использования их в методе PL-оценок.

4. Методика реконструкции эффективного спектра тормозного излучения на основании обработки пространственных дозных распределений, измеренных при вводе в эксплуатацию радиотерапевтических ЛУЭ и кобальтовых облучателей.

Апробация работы. Основные материалы и положения работы докладывались и обсуждались на следующих научных конференциях и семинарах:

Научная конференция в Московском инженерно-физическом институте «Научная сессия МИФИ 2004, 2005». – Москва, 2004г, 2005 г.

II Евразийский конгресс по медицинской физике и инженерии «Медицинская физика 2005», Москва, 2005 г.

Публикации.

Материалы диссертации опубликованы в 5 печатных работах: 1 статье в научных журналах, 4 научных трудах и тезисах конференций. 1 статья принята к изданию в журнале «Атомная энергия».

Структура и объем диссертации. Диссертация состоит из введения, трех глав, заключения, списка литературы и приложения. Общий объем составляет 131 страницу печатного текста, включая 50 рисунков. Список литературы включает 42 наименования, из них 30 иностранных.

СОДЕРЖАНИЕ РАБОТЫ

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

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

Описаны существующие программы, реализующие метод Монте-Карло для лучевой терапии. Для каждой программы показаны её преимущества перед стандартной реализацией метода Монте-Карло. Приведен список основных производителей систем дозиметрического планирования лучевой терапии, у которых запланирована (или уже введена) реализация метода Монте-Карло.

Показаны существующие методы ускорения расчетов методом МонтеКарло. В главе рассмотрены наиболее часто встречающиеся способы ускорения расчетов, показаны их условия применения и уменьшение времени расчета вследствие их применения. Использование данных способов напрямую влияет на ускорение расчетов методом Монте-Карло, однако их применение также косвенно влияет на ускорение расчетов с помощью метода PL-оценок.

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

Вторая глава настоящей работы посвящена исследованию метода PLоценок. Суть PL-оценок состоит в том, что внутри традиционной схемы метода Монте-Карло используются новые оценки для вычисления вклада в рассматриваемые функционалы. При этом преодолеваются два существенных недостатка, присущие методу Монте-Карло: медленная сходимость метода и значительная флуктуация результатов расчетов. При использовании PLоценок скорость счета повышается и исчезают флуктуации в результатах расчетов, так как методу PL-оценок присуще такое замечательное свойство как внутренне сглаживание расчетных данных. Расчеты с использованием PL-оценок также можно комбинировать с другими быстрыми приближенными методами расчета доза, основанными либо на использовании экспериментальных данных, либо на методе суперпозиции/свертки и применять их только для вычисления только в тех точках, где приближенные методы дают недопустимую погрешность.

Расчет поглощенной дозы с использование PL-оценок делится на два этапа:

• построение дерева траекторий в бесконечной однородной водной • PL обработка дерева траекторий.

При построении дерева траектории используются следующие правила:

частица источника – фотон, фазовые координаты для всех частиц источника (r0, E0 ). Непосредственно в точке источника происходит взаимодействие часr тицы с веществом, т.е. построение траектории происходит для источника первых столкновений.

Построение дерева траекторий может производиться любой программой, реализующей метод Монте-Карло для фотон-электрон-позитронных каскадов.

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

рии совпадают Eisimilar = Eiprimary = Ei 2. Оптические длины каждого звена основной и подобной траектории совпадают.

Для построения обратной подобной траектории, проходящей через точr ку детектора rD, на текущем i-м звене основной траектории выбираем случайную точку на текущем звене траектории. Затем совмещаем выбранную точку с точкой детектора.

В случае воксельной геометрии уравнение поиска расстояния до узловой точки li выглядит следующим образом N – число пересеченных вокселей, k – плотность водной среды в k-м вокселе, lk –расстояние пройденное в k-м вокселе, вдоль направления (ui ).

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

C pGy –переводной множитель из МэВ/г в пГр, wi -вес звена на траектории.

Функционал которого перенос заряженных частиц не учитывается, 0 если текущее звено – трек гамма-кванта, для которого перенос заряженных частиц учитывается;

текущего звена траектории, Esource = Ep(E )dE –средняя энергия источника в MэВ. – площадь виртуального источника. Функция источника h r0similar, E0 представляется в виде дой, l расстояние в от точки первого взаимодействия r0 similar до границы воксельной геометрии в направлении обратном движению первичного гампГр см Исходя из этой модели, метод обобщается на случай гетерогенной среды следующим образом. Во-первых, при построении обратных подобных траекторий использовать не плотность вещества, а «электронно-эквивалентную»

плотность, которая вычисляется исходя из соотношения:

Поглощенная доза, пГр*см^2/МэВ Рисунок 1 Глубинное распределение поглощенной дозы от источника 6 MV соответственно в воде. Расположение химической неоднородности {2 см < z< 5 см}, состав неоднородности – алюминий.

Во-вторых, для учета поглощения излучения веществами с Z отличными от воды, необходимо в формулу (2) ввести следующий коэффициент • µ aj – массовый коэффициент поглощения вещества в детекторной точке • µ awater – массовый коэффициент поглощения в воде • E0 – энергия источника. Для энергетически распределенных источников используется средняя энергия спектра.

По этой методике рассчитаны глубинные дозовые распределения, результаты сравнивались с результатами, посчитанными по программе NovaXYZ. (рис.1) Для того чтобы наиболее эффективно применять метод PL-оценок, в настоящей работе предложено использовать способ расчета, который является комбинацией электронного и фотонного способов. Введем понятие уровня взаимодействия для вторичной частицы.

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

Ниже приведена зависимость относительного времени расчета дозы от выбора граничного уровня взаимодействия.

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

ОТНОСИТЕЛЬНОЕ ВРЕМЯ РАСЧЕТА

УРОВЕНЬ ВЗАИМОДЕЙСТВИЯ

уровнях взаимодействия для источника излучения с энергией 5 МэВ Пусть имеется три подряд идущих электронных вершины с координатами {Ei, ri } i=1, 2, 3, которые являются вершинами двух подряд идущих электронных треков. Расстояния между вершинами определяется по формуле li j = ri rj. Для того чтобы исключить из выборки среднюю вершину, т.е.

рассматривать два последовательных трека как один, необходимо выполнение двух условий:

Поглощенная энергия на сглаженном треке l13 вычисляется исходя из формулы:

Были выполнены расчеты дозовых распределений для различных источников излучения с применением описанной методики для k=[1...1,10] c шагом 0,1. Зависимость времени расчета от k показана на рис. 3.

Относительное время расчета Рисунок 3 Зависимость времени расчета от критерия сглаживания для энергии E= 6MV Т.к. в методе PL-оценок используется набор траекторий, рассчитанный от мононаправленного источника, то угол рассеяния после первого взаимодействия фотона со средой будет отличаться от истинного угла рассеяния, который зависит помимо всего прочего от направления начального излучения.

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

Рисунок 4 Дозовые профили, рассчитанные с использованием технологии смещения угла и без нее. Источник 1,25 МэВ, SSD 80 см, размер поля 10x10 см.

С применением описанной технологии учета расходимости первичного излучения были рассчитаны профили дозовых распределений от различных источников излучения. Для полей 10x10 см и SSD 100 см различие между профилями, рассчитанными с использованием сдвига и без него, незначительны. Разница растет с уменьшением SSD, т.к. поток энергии все более и более расходится. В качестве примера на рис. 4 приведены дозовые профили от источника 1,25 МэВ, SSD 80 см, размера поля 10x10. Как и следовало ожидать, применение описанной технологии «поднимает» дозу на краях плато дозового профиля.

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

Рассмотрим сначала метод восстановления спектров тормозного излучения, который осуществляется на основе комбинированного использования нескольких глубинных дозных профилей в плоскостях перпендикулярных оси падения пучка. Если измерения Pi k проведены в плоскостях на глубинах k = 1,2… и известны соответствующие функции отклика Fi k (E ), где E - значения энергий искомого спектра, то для каждого к – го уравнения Pi k = ( E ) Fi k ( E ) dE или в дискретном представлении Дозовые профили для моноэнергетических фотонов – ядра Fi k ( En ) уравнения (9) - рассчитывались методом Монте-Карло помощью кода DPM в диапазоне энергий 0.120 МэВ для геометрических параметров, указанных ниже. Статистика расчета ядер составляла 50·106 историй, статистическая погрешность результатов не превосходила 12%. Исходные профили Pi k так же рассчитывались по той же программе.

Расчеты проводились для трех плоскостей k = 1,2,3 (водный фантом 303030 см3, глубины 5см, 15см, 25см, SSD = 100 см и площадь облучения S = 1010 см2); полученные для первичного спектра начального приближения и для каждой итерации три спектра усреднялись k = ( 1 + 2 + 3 ) / 3 и итерационный спектр k являлся исходным для +1 итерации. Естественным обоснованием этой процедуры является независимость формы восстановленного спектра от используемого глубинного дозного профиля с индексом k.

После проведения расчета комплекса итераций для всех трех глубин дозных профилей определялось относительное отклонение R(к) 10 центральных точек исходного профиля Р(к,i) от соответствующих значений S(к,i), вычисляемых в ходе итерационного процесса:

n(E),отн.ед.

Рисунок 5. Результаты восстановления спектра тормозного излучения с граничной энергией 6 MV (нормированы по площади).

Выбиралась величина максимального отклонения R(к)max и с учетом знака отклонения изменялся параметр начального приближения в сторону уменьшения или увеличения (на 10%); формировалось новое начальное приближение и вся расчетная процедура повторялась до достижения минимума значения R(к)max. Результаты восстановления спектра показаны на рис. Суть другого метода моделирования источника состоит в следующем.

Мы имеем набор координат частиц (x,y,E,u,v,w), которые должна получить программа метода Монте-Карло от программы моделирования вторичного источника, однако для определения направляющих косинусов (u,v,w) достаточно двух углов и (полярный и азимутальный углы, соответственно).

Итак, необходимо сформировать функции распределения вероятности для пяти переменных (x,y,E,,). В настоящем случае обработки файлов фазового пространства в функции распределения вероятности порядок в иерархии переменных подразумевает, что для данного интервала переменной x суммируются все попавшие в интервал фотоны (независимо от y, E,, ); для пары интервалов переменных (x,y) суммируются все попавшие в интервал фотоны (независимо от E,, ); для тройки интервалов (x,y,E) суммируются все попавшие в интервал фотоны (независимо от, ); и так далее.

Для каждого считанного фотона определяются номера ячеек, номер энергетической группы. Выполняется преобразование направляющих косинусов (u,v,w) (u’,v’,w’) для уменьшения анизотропии по углу.

По окончании чтения фазового пространства производится обработка накопленных данных, которая заключается в суммировании и нормировке накопленных данных, в результате чего получается набор функций плотности вероятности (в групповом виде) для соответствующих переменных. Одновременно, еще до нормировки данных, производится накапливающееся суммирование с целью получения функций распределения вероятности в групповом виде. По окончании получения функций плотности вероятности и функций распределения вероятности последние преобразуются в равновероятные интервалы распределения с тем же количеством интервалов для каждой переменной. На рис.6 показаны результаты работы описанного алгоритма Рисунок 6. "Дифференциальные" энергетические распределения фотонов установки Rocus, идущих через ячейки (iX,iY) поверхности запоминания фазового пространства:

черный – (iX,iY) = (1,1); красный – (iX,iY) = (4,4); синий – (iX,iY) = (5,5).

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

ЗАКЛЮЧЕНИЕ

В настоящей работе исследован метод PL-потока в точке и предложены его улучшения, направленные на увеличение скорости и точности расчетов.

Отметим, что скорость счета с использованием PL-оценок имеет почти линейную зависимость от количества детекторных точек. Это позволяет практически сразу рассчитать дозу в любой критической точке за время от 0,04 до 1,5 сек., в зависимости от условий расчета. Использование обычных схем в подобных ситуациях практически невозможно, т.к. зависимость скорости от числа детекторных точек незначительна. Эта особенность делает метод уникальным для проведения расчетов связанных с небольшим количеством точек. Обычно такие расчеты проводятся с исследовательскими или научными целями, где нужно получить глубинные дозовые распределения или профили дозовых распределений.

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

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

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

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

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

ВЫВОДЫ

1. Предложены методы ускорения метода PL-оценок потока в точке, основанные на обработке дерева траекторий.

2. Показана возможность использования метода PL-оценок потока в точке, при расчетах доз в средах с неоднородностью по атомному номеру