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

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

МАТЕМАТИЧЕСКИЕ МОДЕЛИ ПРОЦЕССА ПОГЛОЩЕНИЯ ТЕРАПЕВТИЧЕСКИХ ПУЧКОВ В ТКАНЕЭКВИВАЛЕНТНЫХ СРЕДАХ

05.13.18 — Математическое моделирование, численные методы и комплексы

программ

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

Гордеев Дмитрий Федорович

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

2 8 НОЯ ¿№

Санкт-Петербург

2013

005540743

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

Научный руководитель: доктор физико-математических

наук, профессор Овсянников Дмитрий Александрович

Официальные оппоненты: доктор физико-математических

наук, профессор Жабко Алексей Петрович (Санкт-Петербургский государственный университет, заведующий кафедрой теории управления)

кандидат физико-математических наук, доцент Елизарова Марина Владиславовна (Научно-исследовательский институт онкологии имени H.H. Петрова)

Ведущая организация: Федеральное государственное

унитарное предприятие «Научно-исследовательский институт электрофизической аппаратуры им. Д.В. Ефремова»

Защита состоится 18 декабря 2013 г. в 16 часов на заседании диссертационного совета Д 212.232.50 по защите диссертаций па соискание ученой степени кандидата наук, на соискание ученой степени доктора наук при Санкт-Петербургском государственном университете по адресу: 198504, Санкт-Петербург, Петродворец, Университетский проспект, 35, ауд. 327.

С диссертацией можно ознакомиться в Научной библиотеке имени М. Горького Санкт-Петербургского государственного университета но адресу: 199034, Санкт-Петербург, Университетская наб., 7/9. Автореферат размещен на сайте www.spbu.ru.

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

Ученый секретарь диссертационного совета, доктор физ.-мат. наук, профессор

Курбатова Г. И.

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

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

Появление новых технических решений, таких как мпоголепсстко-вый коллиматор, томотераиия (tomotherapy), ротационная интенсивно-модулированная радиотерапия (intensity-modulated arc therapy) требуют от систем для реализации лучевой терапии нового подхода к решению задач планирования. Однако в России в настоящее время все еще не создана завершенная трехмерная система планирования лучевой терапии, поддерживающая современные технологии облучения пациентов.

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

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

1. Построение математической модели дозного распределения в ткансэквивалснтной среде.

2. Разработка методов решения прямой и обратной задачи планирования лучевой терапии.

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

4. Разработка и тестирование программного обеспечения, предназначенного для планирования лучевой терапии.

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

'Поглощенной ДОЗОЙ, или просто дозой, в лучевой терапии называют отношение поглощенной энергии ионизирующего излучения к массе поглощающего вещества. За единицу измерении поглощенной дозы п системе СИ принят грей (Гр).

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

Практическая значимость. Разработанная модель дозпого распределения в тканеэквивалсптной среде, а также методы и алгоритмы расчета поглощенной дозы могут быть использованы при создании новых версий системы планирования лучевой терапии СКАНПЛАН (СПбГУ).

Методы исследования. В работе основными методами исследования являются методы математического и компьютерного моделирования, численные методы.

Результаты, выносимые на защиту.

1. Математическая модель дозпого распределения в ткансэквивалептной

среде.

2. Алгоритм решения прямой задачи планирования лучевой терапии.

3. Алгоритм решения обратной задачи планирования лучевой терапии.

4. Комплекс программ для планирования лучевой терапии.

Апробация работы. Основные результаты диссертационной работы докладывались на 41-ой, 43-ей и 44-ой международных научных конференциях аспирантов и студентов «Процессы управления и устойчивость» (СПбГУ -2010, 1012 и 2013 гг.), доклад па 23-ей всероссийской конференции по ускорителям заряженных частиц ЯиРАС 2012 (Санкт-Петербург, 24 - 28 сентября 2012).

Публикации. Матсри<1лы диссертации опубликованы в 4-х работах [1— 4], из которых 1 [1] является статьей в журнале, входящем в Перечень ведущих рецензируемых научных журналов и изданий ВАК.

Структура и объём работы. Диссертация состоит из введения, пяти глав, заключения и списка литературы, включающего 75 наименований. Объём составляет 93 страницы.

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

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

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

§1.1 посвящен физике поглощения излучения. Процесс поглощения энергии заключается в столкновении между фотоном и каким-либо электроном

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

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

Во второй главе представлены математические модели дозного распределения в ткансэквивалентпой среде. Базовой моделью является модель точечного мононаправленного источника (§2.1). Данная модель описывает распространение в нолубсскопсчной среде мопопанравленпого потока частиц, падающего нормально к поверхности.

Модель точечного моноиаправлеииого источника. Введем декартову систему координат (x',y',z'). Пусть fl = {(.т',?/, г') | х! 6 (-оо;+оо),?/ е [0;+оо),г' 6 (—оо;+оо)} - облучаемый объем. Тогда, разбив поверхность S = {{x',y',z ) | х' 6 (—оо; + оо), у' = 0, z' € (-оо;+оо)} облучаемого объема на элементарные площадки AS, представим дозу в точке (х1, у', z') € fl в виде суммы доз элементарных источников. Дозовую функцию / точечного мононаправленного источника можно определить как

f(x',y',z',x,z) = \unQf(x\y',z',x,z,AS), (1)

где f(x',y',z',x,z,AS) - доза созданная в точке (x',y',z') в П частицами, прошедшими через элементарную площадку AS с центром в точке (х, 0, z) 6 S.

Умножим функцию f(x',y',z',x,z) на плотность падающего па поверхность среды потока частиц j{x,z) и проинтегрируем по всей облучаемой поверхности 5, в результате получим D - распределение доз в точке (а;', у', z') €

П:

D(x\ у', z') = jj j{x, z)f(x', y', z', x, z)dS. (2)

s

При равномерном стационарном распределении частиц по поверхности (т.е. при j{x,z) = jo) в силу однородности облучаемой среды функция D(x', у', z') не должна зависеть от неременных х' и z'. С учетом этого для бесконечно широкого нучка получаем:

A»(2/')=jo Jj f(y',x,z)dS. (3)

5=00

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

Модель дозпого распределения в случае прямоугольного поля. Поскольку расчет дозовой функции точечного мопопаправлеппого источника / на основе физических представлений о взаимодействии терапевтического пучка с веществом требует существенно больше времени на компьютерные вычисления, чем это допустимо в клинической практике, в работе предлагается определение функции / с помощью предварительно обработанных экспериментальных данных (§2.2). Осуществить это позволяет выражение дозовой функции / точечного мононаправленного источника через экспериментальное глубинное (зависящее только от у') распределение доз достаточно широкого терапевтического пучка Dœ(y'), измеренное в реальных условиях формирования пучка излучения. Кроме того, анализ результатов фотометрии узких пучков, а также хорошее согласие экспериментальных данных и результатов расчета дозных полей в широком диапазоне условий облучения, показали, что поперечную компоненту дозовой функции / точечного мононаправленного источника можно представить в виде функции Гаусса. Таким образом, дозовая функция f точечного мононаправленного источника имеет вид:

J2 (x'-zt)2 + (z'-iS)2

/Or', у', z1, x, z) = ——D^e , (4)

7ТГ*{у')

где 5 - геометрический фактор ослабления, а г (у') - параметр функции Гаусса, характеризующий дозиое распределение.

Отсюда следует, что дозное распределение для прямоугольного поля размером 2А х 2В при равномерном облучении поверхности имеет вид:

D{x',y',z') = D°°4(y ^ {Ф(А(Л<5 + х')) + Ф(А(А<5 - ж'))} х

х{Ф(А(35+2')) + Ф(А(Ш-г'))}, (5)

где Ф(х) = f0x е edt, А = 0-jÇy.

В §1.1 влияние излучения на дозиое распределение представлено в виде суммы двух компонент, первичной и рассеянной. Выражение (5) описывает первичную компоненту излучения. Однако, анализ экспериментальных данных показал, что это выражение также хорошо подходит для описания рассеянной компоненты излучения. Таким образом, общая формула, учитывающая первичную и рассеянную компоненту поля будет иметь вид:

D(x', у\ z') = £^11{Ф(Х(А0 + х')) + Ф(\(А6 - i'))} х х{Ф(А(б<5 + z')) + Ф(Л(BÔ - г'))} + + О£ооЬЛ{ф{к{м + ^ + ф(х,(А6 - ®'))} х

х{Ф(А.,(М + z')) + Ф{Xs(BS - г'))}, (6)

где As = у^, a rs(y') - параметр функции Гаусса, характеризующий дозиое распределение рассеянной компоненты.

Модель дозиого распределения в тхаиеэквивалеитной среде. В §2.3 используется полученное выражение (6) для моделирования нолей произвольной формы. Для этого область облучения в горизонтальной плоскости, касающейся поверхности срсды, разбивается на конечное число прямоугольников со сторонами, параллельными координатным осям. Каждому такому прямоугольнику со сторонами 2А х 2В соответствует прямоугольное ноле (6). Таким образом, целый пучок представляется в виде набора составных пучков, называемых в работе карандашными. Однако, центральные оси карандашных пучков пересекают касательную плоскость под некоторым углом. Пусть г'-ый карандашный пучок пересекает касательную плоскость в точке с координатами (xí,0,Zí). Положим:

X ' Z ' /

at = arctg-£, /?; = arctg-^, с = \Jx} + zf + R, (7)

где R - расстояние источник-поверхность. Тогда координаты (х, у, z) точки M € П в системе, связанной с г-м карандашным пучком, вычисляются но формулам:

X = {х' - Xi) COS ni - у' Sino-i, Z = (2' - Zi) cos fa - y'sin fa,

y = -{x'xi + z'zí + R(y' + R)). (8)

с

Кроме того, исходный прямоугольник со сторонами, равными А и В, при проецировании на плоскость, ортогональную центральной оси г-го_карандашного пучка, переходит в параллелограмм со сторонами, равными A¡ и Bf.

Ài = A cos a-^j, Bi = В cos (9)

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

7 П

0(х',у',г') = ^20{(х1,у'>г'), (10)

¿=1

где т - количество карандашных пучков, а Д(ж', у', г') - доза в точке х', у', г', созданная г-м карандашным пучком.

Перечислим основные преимущества предложенной модели (1) - (10) перед известными в настоящее время моделями дозного распределения в тка-неэквивалентной среде. Можно выделить три типа математических моделей дозного распределения в тканеэквивалентной среде. К первому типу относятся модели, основанные на расчете поглощенной дозы путем интерполяции большого объема дозиметрических данных, полученных для широкого спектра условий облучения. Модели первого типа плохо пригодны для моделирования полей, сформированных МЛК, из-за технических трудностей в процессе дозиметрии узких пучков.

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

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

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

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

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

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

В §3.3 предложен алгоритм прямого планирования, основанный па модели из §2.3. Для того чтобы получить дозное распределение в облучаемом объеме, необходимо построить сетку, состоящую из прямоугольных параллелепипедов одинакового размера. При этом грани этих параллелепипедов параллельны соответствующим координатным плоскостям системы координат Ох'у'г'.

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

Рис. 1 — Проекции лепестков коллиматора па плоскость Охг

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

В §3.4 представлен алгоритм обратного планирования. Данный алгоритм разработай для решения следующей задачи. Сеанс лечения осуществляется шагами, на каждом шаге производится несколько последовательных облучений с разной конфигурацией МЛК. Между шагами происходит поворот гаитри терапевтического аппарата на некоторый угол. Пусть количество шагов и углы поворота гантри заданы. Необходимо найти такие конфигурации МЛК, чтобы минимизировать количество облучений на каждом шаге и при этом получить заданное дозное распределение.

Пусть I - матрица интенсивности карандашных пучков, элементы которой являются натуральными числами. Рассмотрим первый частный случай, в котором матрица / состоит из нулей и единиц. Обозначим многоугольник, состоящий из проекций карандашных пучков ненулевой интенсивности, через Р (рис. 2). Многоугольник Р ограничен вертикальными и горизонтальными ребрами, а также может содержать «дыры». Назовем вертикальное ребро е многоугольника Р левым пиком, если оба примыкающих к нему горизонтальных ребра находятся слева от е. Аналогично правым пиком назовем вертикальное ребро многоугольника Р, у которого примыкающие горизонтальные ребра находятся справа. Само наличие по крайней мере одного правого или левого пика свидетельствует о невозможности создания поля, обладающего формой Р, с помощью одной конфигурации лепестков МЛК.

(а! <Ъ>

Рис. 2 — Многоугольник Р, состоящий из проекций карандашных пучков ненулевой интенсивности

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

1. Построим направленный граф <3, узлами которого являются вершины проекций карандашных пучков, лежащие внутри многоугольника Р или на его пиках. Ребрами графа С будут стороны проекций карандашных пучков, соединяющие узлы С. Направления ребер С? определяются по правилам:

(a) Ребро, левый конец которого принадлежит левому пику, имеет па-правление от левого пика.

(b) Ребро, правый конец которого принадлежит правому пику, имеет направление к правому пику.

(c) Направление горизонтальных ребер, концы которых не принадлежат пикам, имеют направление слева направо.

(с1) Вертикальные ребра, концы которых не принадлежат пикам, являются двунаправленными.

2. Для всех узлов графа б, принадлежащих одному левому (правому) пику, введем узел-исток (узел-сток) и ребра соединяющие этот узел-исток (пиковый узел) с пиковым узлом (узлом-стоком).

3. Присвоим пропускную способность каждому узлу и ребру графа б.

4. Вычисляем максимальный поток па направленном графе (7 от узлов-истоков до узлов-стоков, на основе полученного максимального потока соединяем между собой левые и правые пики Р.

5. Удаляем оставшиеся несвязные пики.

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

Во втором частном случае максимальный элемент матрицы интенсивности I достаточно мал, т.е. не больше наперед заданного числа, которое является параметром алгоритма. Пусть Б = {,%} - будет искомое множество конфигураций МЛК (множество сегментов), составляющих матрицу интенсивности I. Так как максимальная интенсивность каждого сегмента равна 1, то пересечение Б, с любым столбцом С,- матрицы / представляет собой непрерывный блок интенсивности 1. Назовем такой непрерывный блок полевым открытием, и обозначим его как В^

Рассмотрим множество полевых открытий {В„•}, образованное всеми сегментами на столбце Су Назовем такое множество {Ву} вариантом доставки соответствующего профиля интенсивности. Таким образом, для второго частного случая необходимо сначала найти оптимальные варианты доставки для каждого профиля интенсивности, а затем «сшить» вместе нолевые открытия, принадлежащие соседним столбцам матрицы интенсивности, чтобы сформировать сегменты.

Общий алгоритм решения обратной задачи планирования состоит из следующих шагов:

1. Пусть I - матрица интенсивности карандашных пучков, элементы которой гпт, п = 1,..., Л/, т = 1,..., М являются натуральными числами, а

1тах ~ максимальный уровень интенсивности, т.е. гпт < 1тах. Если 1тах достаточно мал, то выполняем шаг 5а.

2. Выберем уровень интенсивности й, такой что 1 < й < 1тах.

3. Построим матрицу интенсивности /, элементы которой 1ПТП, п= 1,..., т = 1,..., М удовлетворяют условию:

Т _ /о, гпт < й

Чт - 1пт><1

4. Применяем алгоритм первого частного случая для матрицы I.

5. Вычтем матрицу / из исходной матрицы /, полученную матрицу / = 1—1 следует рекурсивно обработать согласно условию:

(a) Если максимальный элемент матрицы I достаточно мал, то следует применить алгоритм для второго частного случая.

(b) Если максимальный элемент матрицы I не достаточно мал, то следует повторить шаги 2-5 теперь уже для матрицы I.

В четвертой главе описывается программное обеспечение для планирования лучевой терапии, созданное на основе разработанных моделей и алгоритмов. Программное обеспечение написано на С+ I- и имеет модульную архитектуру (рис. 3).

Программный модуль подготовки дозиметрических данных «Инициализация» является рабочим местом медицинского физика и предназначен для получения и обработки физических данных терапевтических установок и их дозиметрических характеристик (§4.1). Модуль предполагает работу с двумя видами представления дозиметрических данных - глубинными распределениями и профилями. Глубинные распределения отражают зависимость поглощаемой дозы от глубины на центральной оси пучка. Для каждого размера поля (ширины пучка) строится свое глубинное распределение. Профили, в свою очередь, представляют собой зависимости поглощаемой дозы от расстояния до оси пучка на фиксированной глубине. Для каждой глубины и каждого размера поля строится свой профиль.

Программный модуль подготовки анатомо-топометрической информации «Анатомия» предназначен для ввода и редактирования анкетных данных и анатомических срезов по пациентам (§4.2).

Программный модуль подготовки планов облучения «План» предназначен для создания, редактирования и расчета планов облучения по пациентам (§4.3). Создание планов облучения производится следующим образом: на полученную в модуле «Анатомия» группу анатомических срезов (анатомию) накладывается схема облучения, т.е. задаются форма, расположение и другие параметры облучающих пучков, после чего производится расчет дозных

Данные 1 | Анатомии | ; Параметры

пациента 1 установки '

Рис. 3 — Схема взаимодействий программных модулей

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

щ-

7

Рис. 4 — Изодозные кривые: ротационный режим, сектор = 180°, поле 10 х 10 см, РИО = 100 см

Для ввода анатомических срезов и сохранения готового плана облучения используется стандарт Б1СОМ. В §4.4 приводится описание стандарта Б1СОМ и его поддержки в разработанном программном обеспечении.

С помощью разработанного программного обеспечения были рассчитаны дозные распределения пучков, созданных линейным ускорителем и гамма-

Рис. 5 — Дозный профиль: многолепестковый коллиматор, глубина - 5 см

аппаратом. При расчете использовались различные конфигурации облучения: открытые поля, пучки с экранирующими блоками и клиновидными фильтрами, ротационная терапия, пучки с МЛК (рис. 4, 5).

В пятой главе приводятся результаты экспериментальной проверки разработанных моделей, алгоритмов и программного обеспечения. В качестве экспериментальных данных были использованы дозные распределения для широкого спектра условий облучения, измеренные в водном фантоме. Измерения и расчет были осуществлены как для линейного ускорителя ЭЛЛУС-6М, разработанного НИИЭФА им. Ефремова, с энергией электронов 6 МэВ, так и для гамма-аппарата Theratron Equinox, выпускаемого компанией Best Theratronics (Канада). Отклонения рассчитанных дозных распределений от полученных экспериментально лежат в допустимом интервале.

Публикации автора по теме диссертации

Статьи в журналах, рекомендованных ВАК

1. Гордеев Д.Ф. Моделирование и расчет дозного распределения в тка-нсэквивалентной среде // Вестник С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2013. Вып. 3. С. 142—149.

Другие публикации

2. Гордеев Д.Ф. Сравнение доз для фотонов, рассчитанных системами лучевого планирования СКАНПЛАН и R.OCS // Процессы управления и устойчивость: Труды 41-й международной научной конференции аспирантов и студентов / Под ред. Н. В. Смирнова, Г. Ш. Тамасяна. СПб.: Издат. Дом С.-Петерб. гос. ун-та, 2010. С. 316-319.

3. Гордеев Д.Ф. Обработка дозиметрической информации для инициализации системы планирования лучевой терапии СКАНПЛАН // Процессы управления и устойчивость: Труды 43-й международной научной конференции аспирантов и студентов / Под ред. А. С. Ерёмина, Н. В. Смирнова. СПб.: Издат. Дом С.-Петерб. гос. ун-та, 2012. С. 277—282.

4. Гордеев Д.Ф. Метод расчета дозного распределения в тканеэквивалеит-ной среде // Процессы управления и устойчивость: Труды 44-й международной научной конференции аспирантов и студентов / Под ред. Н. В. Смирнова, Т. Е. Смирновой. СПб.: Издат. Дом С.-Петерб. гос. ун-та, 2013. С. 329—334.

Подписано к печати 12.11.13. Формат 60 х 84 '/¡6. Бумага офсетная. Гарнитура Тайме. Печать цифровая. Печ. л. 1,00. Тираж 100 экз. Заказ 5917.

Отпечатано в Отделе оперативной полиграфии химического факультета СПбГУ 198504, Санкт-Петербург, Старый Петергоф, Университетский пр., 26 Тел.: (812) 428-4043,428-6919

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

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

г;/ тл 1 счпт и-гьу I -т^^иьи

Гордеев Дмитрий Федорович

МАТЕМАТИЧЕСКИЕ МОДЕЛИ ПРОЦЕССА ПОГЛОЩЕНИЯ ТЕРАПЕВТИЧЕСКИХ ПУЧКОВ В ТКАНЕЭКВИВАЛЕНТНЫХ СРЕДАХ

05.13.18 - Математическое моделирование, численные методы и комплексы

программ

ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

Научный руководитель доктор физико-математических наук, профессор Овсянников Д.А.

Санкт-Петербург 2013

Оглавление

Введение................................................................................................................3

Глава 1. Физика и анализ дозного распределения............................................6

1.1. Физика поглощения излучения.................................................................6

1.2.Анализ дозного распределения и рассеяния..........................................10

Глава 2. Математическое моделирование дозного распределения в тканеэквивалентной среде.................................................................................23

2.1. Модель точечного мононаправленного источника..............................23

2.2. Моделирование дозного распределение в случае прямоугольных полей облучения..............................................................................................27

2.3. Моделирование дозного распределение в случае полей облучения произвольной формы......................................................................................30

Глава 3. Алгоритмы...........................................................................................33

3.1. Прямая и обратная задачи планирования лучевой терапии................33

3.2. Алгоритмы расчета дозы.........................................................................33

3.3. Алгоритм прямого планирования...........................................................39

3.4. Алгоритм обратного планирования.......................................................40

Глава 4. Программное обеспечение.................................................................50

4.1. Обработка дозиметрической информации............................................52

4.2. Обработка анатомо-топометрической информации............................59

4.3. Подготовка плана облучения и его расчет............................................62

4.4. Поддержка стандарта DICOM................................................................65

Глава 5. Экспериментальная проверка............................................................74

5.1. Требования к точности планирования лучевой терапии.....................74

5.2. Линейный ускоритель СЛ-75/ЭЛЛУС...................................................74

5.3. Гамма-аппарат Theratron Equinox...........................................................78

Заключение.........................................................................................................86

Библиографический список..............................................................................87

Введение

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

Появление новых технических решений, таких как много лепестковый коллиматор, томотерапия (tomotherapy), ротационная интенсивно-модулированная радиотерапия (intensity-modulated arc therapy) требуют от систем для реализации лучевой терапии нового подхода к решению задач планирования. Однако в России в настоящее время все еще не создана завершенная трехмерная система планирования лучевой терапии, поддерживающая современные технологии облучения пациентов.

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

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

• построение математической модели дозного распределения в тканеэквивалентной среде;

• разработка методов решения прямой и обратной задачи планирования лучевой терапии;

• разработка методов поиска оптимальных параметров облучения;

• разработка и тестирование программного обеспечения, предназначенного для планирования лучевой терапии;

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

Методы исследования. В работе основными методами исследования являются методы математического и компьютерного моделирования, численные методы.

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

Практическая значимость. Разработанная модель дозного распределения в тканеэквивалентной среде, а также методы и алгоритмы расчета поглощенной дозы могут быть использованы при создании новых версий системы планирования лучевой терапии СКАНПЛАН (СПбГУ).

Положения, выносимые на защиту.

1. Математическая модель дозного распределения в тканеэквивалентной среде.

2. Алгоритм решения прямой задачи планирования лучевой терапии.

3. Алгоритм решения обратной задачи планирования лучевой терапии.

4. Комплекс программ для планирования лучевой терапии.

Апробация работы. Основные результаты диссертационной работы докладывались на 41-ой, 43-ей и 44-ой международных научных конференциях аспирантов и студентов «Процессы управления и устойчивость» (СПбГУ - 2010, 1012 и 2013 гг.), доклад на 23-ей всероссийской конференции по ускорителям заряженных частиц ЯиРАС 2012 (Санкт-Петербург, 24 - 28 сентября 2012). По результатам диссертационных исследований опубликовано 4 работы.

Глава 1. Физика и анализ дозного распределения

1.1. Физика поглощения излучения

1.1.1. Рассеяние энергии

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

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

1.1.2. Линейный коэффициент ослабления

Предположим, что детектор, который может регистрировать число фотонов, проходящих через него, помещен в точку Р. Пусть число зарегистрированных фотонов равно N. Если на пути пучка фотонов поместить слой вещества толщиной Ах, то количество фотонов, достигающих Р, уменьшится на величину АТУ. Число фотонов АТУ, вышедшее из пучка, непосредственно зависит от числа фотонов, присутствующих первоначально в пучке. Если N увеличить вдвое, то вероятность взаимодействия также увеличится вдвое. Величина А/У непосредственно зависит и от толщины Ах. Если увеличить вдвое Ах, т.е. число атомов вещества на пути излучения, то удвоится и вероятность взаимодействия. Величина АЫ изменяется, как произведение N на Ах. Это можно записать в виде

А/У = -/¿/УАх,

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

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

1.1.3. Экспоненциальный закон ослабления

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

1.1.4. Фотоэлектрическое поглощение

Предположим, что фотон с энергией hv взаимодействует с атомом так, что выбивается электрон с K-оболочки. Этот электрон будет иметь энергию, равную ho —Ек, где Ек - энергия связи K-оболочки. Этот электрон называется фотоэлектроном. При этом виде поглощения вся энергия фотона передается фотоэлектрону, удаляемому из атома. Таким образом, фотон исчезает, а вместо него появляется быстрый электрон, и атом приходит в возбужденное состояние. Через короткий промежуток времени другой электрон заполнит K-оболочку, и при этом произойдет испускание характеристического излучения.

Энергия связи на K-оболочке в ткани составляет порядка 500 эВ, и поэтому фотоэлектрон получит энергию, сопоставимую с энергией первичного фотона, т.е. порядка 1-10 МэВ. Таким образом, если фотон поглощается за счет фотоэффекта, то вся энергия передается фотоэлектрону. Фотоэлектрическое поглощение происходит на связных электронах.

1.1.5. Рассеяние излучения

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

1.1.6. Резонансное поглощение первичных фотонов

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

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

1.1.7. Комптоновское рассеивание

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

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

В биологических тканях в области энергии от 100 кэв до 10 МэВ комптоновское поглощение гораздо более существенно, чем фотоэлектрическое и процесс образования пар.

1.1.8. Образование пар

Когда энергия падающего фотона больше 1,02 Мэв, фотон может поглотиться посредством механизма образования пар. Если фотон проходит около ядра атома, то под действием поля ядра фотон исчезает и образуется пара частиц - электрон позитрон. Поскольку масса покоя электрона равна 0,511 Мэв и образуются две частицы, минимальная энергия фотона должна быть равна 1,022 Мэв. В этом процессе не создается заряд, так как электрон и позитрон имеют противоположные заряды. Если фотон имеет энергию свыше 1,022 Мэв, этот излишек энергии распределяется между позитроном и электроном.

При образовании пар происходит взаимодействие фотона и электрического поля ядра:

• Порог этого процесса равен 1,02 Мэв;

• Вероятность процесса растет с ростом энергии выше этого порога;

• Процесс быстро возрастает с увеличением атомного номера;

• Поглощенная энергия меньше, чем энергия падающего фотона, на 1,02 Мэв;

• При аннигиляции создаются два фотона каждый с энергией 0,511 Мэв.

1.1.9. Относительный вклад различных видов поглощения

Всякий раз, когда происходит взаимодействие между пучком фотонов и тканеэквивалентным поглотителем, один фотон выходит из пучка и один электрон (фотоэлектрон, комптоновский электрон или электронная пара) приобретает кинетическую энергию. Можно рассчитать относительное число этих электронов. При энергии 20 кэв 70% электронов являются фотоэлектронами, а 30% комптоновскими. При 26 кэв образуется по 50% каждого типа электронов, а при 100 кэв - 99% комптоновские, остальные - фотоэлектроны. При 24 Мэв 50% комптоновских электронов и 50 % электрон-позитронных пар.

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

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

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

Таким образом, при энергии фотонов до 50 кэв существенно фотоэлектрическое поглощение, при 60 - 90 кэв фотоэлектрический и комптоновские процессы одинаково важны. При энергии фотонов 200 кэв - 2 Мэв доминирует комптоновское поглощение. При 5-10 Мэв существенным становится эффект образования пар, а при 50-100 Мэв этот эффект - наиболее важный тип поглощения.

1.2. Анализ дозного распределения и рассеяния

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

1.2.1. Фантомы

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

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

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

1.2.2. Глубинное распределение

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

Для вычисления глубинного распределения дозы вдоль центральной оси пучка было введено множество величин, таких как процентная глубинная доза [39], отношение ткань-воздух (tissue-air ratios, TAR) [31, 37, 50, 51], отношение ткань-фантом (tissue-phantom ratios, TPR) [38, 54, 64], отношение ткань-максимум (tissue-maximum ratios, TMR) [38, 56]. Эти величины обычно получаются из измерений, сделанных в водных фантомах, используя маленькие ионизационные камеры.

1.2.3. Процентная глубинная доза

Один из способов определения дозового распределения на центральной оси это нормализация дозы на произвольной глубине (I относительно дозы на заданной глубине с10. Таким образом, процентная глубинная доза Р вычисляется следующим образом:

Р = —^-хЮО,

где - значение дозы на произвольной глубине с!, �