CC BY
732
УДК 57.017
Вестник СПбГУ. Сер. 11. 2013. Вып. 4
И. В. Красников, В. Е. Привалов, А. Ю. Сетейкин, А. Э. Фотиади
РАСПРОСТРАНЕНИЕ ОПТИЧЕСКОГО ИЗЛУЧЕНИЯ В БИОЛОГИЧЕСКИХ ТКАНЯХ
Оптические методы диагностики и терапии находят в настоящее время все большее применение в биологии и медицине [1-7]. Основное преимущество этих методов заключается в их неинвазивности: низкоинтенсивное лазерное излучение в ближнем ИК диапазоне не оказывает вредного воздействия на биологическую среду. Кроме того, помимо структурной информации оптические методы позволяют получать функциональную информацию о биологическом объекте, например: проводить анализ гемодинамики и метаболических процессов (мозговое кровообращение, объем крови, оксигенация мышечной ткани), определять локализацию неодно-родностей (раковых опухолей, разрушения зубной эмали), осуществлять диагностику заболеваний и т. д. [4-10].
В данном обзоре рассмотрены лишь некоторые проблемы взаимодействия света с биологическими тканями, связанные с поглощением и рассеянием света случайными и квазиупорядоченными структурами. Рассмотренные эффекты показывают необходимость детального изучения оптических свойств биотканей с разнообразной структурой.
1. Поглощение света биологическими тканями
Известно, что ослабление электромагнитного излучения при прохождении его через среду определяется двумя основными процессами: поглощением и рассеянием. Поглощательная способность среды определяется как отношение поглощенной и падающей интенсивностей. Поглощение является следствием частичного перехода световой энергии в тепловое движение или колебательную энергию поглощающих молекул. Полностью прозрачная среда не поглощает свет. Полная световая энергия, вошедшая в такую среду и вышедшая из нее, одинаковы. Среди биологических тканей почти прозрачными для видимого света можно считать роговицу и хрусталик глаза. Структуры же, в которых падающее излучение практически полностью ослабляется, называют непрозрачными [8].
Способность вещества поглощать электромагнитное излучение зависит, главным образом, от электронной структуры поглощающих его атомов или молекул (по-
Красников Илья Владимирович — кандидат физико-математических наук, доцент, Амурский государственный университет
Привалов Вадим Евгеньевич — доктор физико-математических наук, профессор, Санкт-Петербургский государственный политехнический университет; e-mail: [email protected]
Сетейкин Алексей Юрьевич — кандидат физико-математических наук, доцент, Санкт-Петербургский государственный политехнический университет
Фотиади Александр Эпаминондович — доктор физико-математических наук, профессор, зав. кафедрой, Санкт-Петербургский государственный политехнический университет; e-mail: Fotiadi@rphf. SPbstu.ru
© И. В. Красников, В. Е. Привалов, А. Ю. Сетейкин, А. Э. Фотиади, 2013
глощающих центров), их концентрации и температуры, длины волны излучения, а также от толщины поглощающего слоя. При взаимодействии падающего света с молекулами среды поглощается квант, равный разности энергий их электронных состояний. Вследствие этого поглощение является спектрально зависимым. В случае, если рассеяние в среде незначительно, то интенсивность излучения (1Т), проходящего через среду, убывает экспоненциально с расстоянием, в соответствии с законом Бугера—Ламберта—Бера, т. е.:
1Т (X) = !0(Х) х ехр(-М), (1.1)
где ца — коэффициент поглощения, равный вероятности поглощения на единицу длины. Коэффициент поглощения, по существу, — это сумма вкладов всех хромофоров среды, т. е.:
мад = 1еД)[С, ], (1.2)
I
где £ и С — коэффициент экстинкции и концентрация ¡-го хромофора. В биологических тканях рассеяние, как правило, гораздо сильнее, чем поглощение, следовательно, закон Бугера—Ламберта—Бера должен быть модифицирован [8].
В биологических тканях основными поглощающими центрами (хромофорами) являются молекулы воды или макромолекулы [10]. Здесь и далее под хромофорами понимаются ненасыщенные группы атомов, обуславливающие поглощающее электромагнитное излучение.
У белков хромофорами являются фрагменты аминокислот, которые поглощают свет преимущественно в ультрафиолетовой области спектра (от 200 до 300 нм). В этом же диапазоне длин волн поглощают нуклеиновые кислоты (их хромофоры — ароматические и гетероциклические кольца азотистых оснований). Клетки биологических тканей содержат сотни хромофоров, поглощающих свет в видимой и ближней ультрафиолетовой областях спектра, среди которых основными являются витамины, флавины, флавиновые ферменты, НАД^Н, гемоглобин, каротиноиды, фико-билины, фитохромы и др. В инфракрасной области спектра все биомолекулы имеют достаточно интенсивные колебательные полосы поглощения. Начиная с X = 1500 нм и более, спектр поглощения тканей в основном определяется спектром поглощения воды.
Рис.1. Коэффициенты поглощения основных поглощающих хромофоров в биологических тканях.
Коэффициенты поглощения рассчитаны для 20 цМ дезоксигемоглобина (НЬ), 80 цМ оксигемоглобина (НЬ02), 70% воды (Н2О) и 15% жира [10].
Гемоглобин. Гемоглобин преобладает в сосудистой ткани. Существует две формы гемоглобина: оксигемоглобин и дезоксигемоглобин. Их спектральные характеристики различны, как показано на рисунке 1, и могут быть использованы для оценки оксигенации тканей [11, 12]. Оксигенация тканей является ценным инструментом при оценке гипоксии, что может быть использовано в онкологии, поскольку вид опухоли можно определить по недостатку кислорода.
Вода. Вода поглощает в основном в ИК области (рис. 1). Вода присутствует в изобилии по всему телу и обычно представлена водной фракцией (% от чистой воды) [13].
Липиды. Липиды (от греч. Проз — жир) — это название многих различных видов жироподобных веществ, входящих в состав всех живых клеток и играющих важную роль в жизненных процессах. Например, липиды играют важную роль в качестве строительных блоков в клеточных мембранах. Подобно воде, липиды поглощают в области более длинных волн (рис. 1). По данным липид-спектроскопии [14], в злокачественных опухолях молочной железы содержание липидов сокращается, о чем свидетельствует сравнительно малый их показатель поглощения относительно показателя поглощения молекул воды.
Меланин. Меланин является основным пигментом кожи и, безусловно, самым главным хромофором эпидермиса. Его коэффициент поглощения монотонно возрастает по всему видимому диапазону спектра с уменьшением длины волны. Меланин выступает в качестве защитного слоя от солнечного УФ-излучения. Таким образом, он находится в поверхностном слое кожи. Меланин отвечает за цвет кожи, но также может быть найден и в других частях тела, например, мозге [15].
Миоглобин. Миоглобин — белок, сходный по строению и функциям с гемоглобином и содержащийся в скелетной и сердечной мышцах, имеет аналогичные с гемоглобином характеристики поглощения. Существуют окси- и дезоксимиоглобин [16].
Цитохромы. Цитохромы являются веществами с общими характеристиками, они содержатся в митохондриях. Спектроскопическая оценка таких веществ может дать информацию о метаболизме исследуемой среды [17-19].
Основной особенностью всех биомолекул является их комплексная структура полосы между 400 нм и 600 нм. Так как ни макромолекулы, ни вода не поглощают сильно в ближнем ИК диапазоне, то «терапевтическое окно» заключено приблизительно между 600 нм и 1200 нм. В этом спектральном диапазоне излучение проникает в биологические ткани с наименьшими затруднениями, что делает возможным лечение достаточно глубоких тканевых структур.
2. Рассеяние света
Ткани состоят из клеток, в которых внутриклеточные органеллы и внеклеточные структуры образуют комплексную матрицу, макроскопически управляющую анатомией органов. Очевидно, существует множество типов тканей с различным составом. Например, соединительные ткани содержат коллаген и эластин, тогда как мышечные ткани состоят из клеток, заполненных миозином и актином. Независимо от того, какая ткань рассматривается, рассеяние возникает за счет изменения показателя преломления среды [20-22]. На рассеяние влияет также размер составляющих ткани частиц (от нм — белки — до мкм — клетки). Этот факт позволяет исполь-
зовать рассеяние как основу диагностических методов для исследования морфологических изменений в ткани [23, 24].
Различают неупругое и упругое рассеяние, в зависимости от того, изменяется в процессе рассеяния начальная энергия фотона или нет [25]. Одним из характерных и важных типов неупругого рассеяния является рассеяние Бриллюэна. Оно возникает при распространении через среду акустических волн, вызывающих неоднородности показателя преломления. Рассеяние Бриллюэна происходит для света с более высокими (или более низкими) частотами, так как рассеивающие частицы движутся навстречу (или удаляются) относительно источника света. Таким образом это может быть рассмотрено как оптический эффект Доплера, когда частота фотонов увеличивается или уменьшается. При взаимодействии лазерного излучения с тканью рассеяние Бриллюэна становится значительным только во время образования ударной шоковой волны.
Отдельным видом упругого рассеяния является рэлеевское рассеяние. Оно имеет место при условии, что длина волны падающего излучения намного превосходит
размеры рассеивающих частиц. В этом слу- -
чае интенсивность 10 рассеиваемого средой г
света обратно пропорциональна 4-й степени *__^
ДЛИНЫ ВОЛНЫ А. падающего Света (1° ~ Аг4) (за- Падающее излучение
кон Рэлея) [8]. На рисунке 2 показана про- <-.
стая геометрия рэлеевского рассеяния. _
Если размер рассеивающих частиц ста- Рассеивающая среда
новится соизмеримым С длиной волны пада- Рис. 2. Геометрия рэлеевского рассеяния [8] ющего излучения, как в случае клеток крови,
закон Рэлея становится неприменим и имеет место другой тип рассеяния, называемый рассеянием Ми. Необходимо особо отметить два важных отличия между рассеянием Ми и рассеянием Рэлея. Во-первых, рассеяние Ми показывает более слабую зависимость от длины волны (~Х-х) по сравнению с рассеянием Рэлея (~Х-4). Во-вторых, рассеяние Ми происходит предпочтительно в направлении вперед, тогда как для рассеяния Рэлея интенсивности света, рассеянного вперед и назад, одинаковы [10].
В большинстве биологических тканей фотоны рассеиваются предпочтительно в направлении вперед. Это явление не может быть объяснено с помощью рассеяния Рэлея. С другой стороны, наблюдаемая зависимость от длины волны более сильная, чем допускает рассеяние Ми. Таким образом, ни рассеяние Рэлея, ни рассеяние Ми не могут полностью описать рассеяние в тканях. Поэтому удобно ввести функцию вероятности р того, что фотон рассеется на угол 0, который может быть подобран по экспериментальным данным. Если р не зависит от 0, говорят об изотропном рассеянии. В противном случае имеет место анизотропное рассеяние.
Характеристикой анизотропии рассеяния является фактор анизотропии g, в случае g = 1 рассеяние происходит только вперед, g = -1 — рассеяние только назад, и если g = 0 — изотропное рассеяние. В полярных координатах фактор анизотропии g определяется как
g =
(2.1)
| рш®
где р(0) — функция вероятности и dw = sin 0d 0d — элементарный телесный угол. По определению, фактор анизотропии g представляет собой средний косинус угла рассеяния 0. Для большинства биологических тканей g лежит в диапазоне от 0,7 до 0,99. Отсюда соответствующие углы рассеяния наиболее часто равны 8-45°. Важной величиной в выражении (2.1) является функция вероятности р. Она также называется фазовой функцией и обычно нормируется следующим образом:
1 г
— I p(9)dro= 1. (2.2)
4л/
4 п
Некоторые теоретические фазовые функции известны как функции Хени— Гринштейна [26], Рэлея—Ганса [27-29], Дельта—Эддингтона и Рейнольда [29-31]. Среди них, в соответствии с экспериментальными наблюдениями, наилучшей является первая. Она была введена Хейни и Гринштейном (1941) и записывается так:
1 - g2
Р(9) = "-~-^ (2.3)
(l + g2 - 2g cos e)
Эта фазовая функция математически очень удобна для использования, так как она эквивалентна представлению
p(9) = Х(2/ +1) giPi (cos 9), (2.4)
i=0
где Pi — полиномы Лежандра. Хотя, в некоторых случаях, сложная функция изотропной величины u и функции Хени—Гринштейна лучше соответствуют экспериментальным данным. Эта модифицированная функция может быть записана следующим образом:
1 u + (1 - u)(l - g2)
4П (l + g2 - 2g cos 9)
3. Распространение лазерного излучения в мутных средах
Математическое описание процессов распространения излучения в мутных средах может быть проведено двумя способами — с помощью аналитической теории и с помощью теории переноса [10]. Первая основывается на уравнениях Максвелла и в принципе является наиболее фундаментальным подходом. Однако его использование ограничено сложностью получения точных аналитических решений. С другой стороны, теория переноса в основном рассматривает перенос фотонов через поглощающие и рассеивающие среды, не основываясь на уравнениях Максвелла. Она имеет эвристический характер, и ей не хватает строгости аналитических теорий. Тем не менее, теория переноса широко используется для описания взаимодействий лазерного излучения с тканью, и экспериментально подтверждено, что во многих случаях ее прогнозы являются достаточными [32-37].
3.1. Теория переноса излучения
Теория переноса, называемая также теорией переноса излучения, берет свое начало с работы Шустера 1903 г. [25]. Теория оперирует непосредственно переносом энергии в среде, содержащей частицы. Сама по себе она не включает дифракционных эффектов. Предполагается, что при суммировании полей отсутствует корреляция между ними, так что складываются интенсивности, а не сами поля. Основное дифференциальное уравнение этой теории называется уравнением переноса или уравнением транспорта и эквивалентно уравнению Больцмана, используемому в кинетической теории газов. В теории переноса можно учесть поляризационные эффекты. Однако в большинстве случаев из соображений математического удобства поляризацией пренебрегают.
Основное уравнение теории переноса излучения формулируется через функцию распределения частиц у(г, 5, £), т. е. число частиц в единице объема при единичном угле, двигающихся в направлении 5 в точке г в момент времени £ [32]. При этом принимается предположение, что энергия частицы при взаимодействии в среднем не изменяется. Таким образом, мы ограничиваем анализ изучением распространения монохроматичных фотонов. В области оптики тканей предпочтительно выражать распределение фотонов в радиометрических (детектируемых) величинах. Энергетическая яркость Ц(г, 5, £) часто используется и определяется умножением функции распределения фотонов на энергию и скорость движения фотонов, V, т. е. скорость света в вакууме, деленную на показатель преломления среды, в которой происходит распространение. Уравнение переноса излучения выглядит следующим образом [25]:
1 Э1(г, 5, £) V д£
+ 5 • У1(г, 5, £) + 5 + ¡а )Ь(г, 5, £) = || Л Цг, 5 ', £) р(5,5 ')ЭО'+д(г, 5, £).
(3.1)
4 п
При рассмотрении бесконечно малого элемента объема среды становится очевидным тот факт, что уравнение переноса излучения выражает условие энергетического баланса в рассматриваемом элементе объема (рис. 3).
Поглощение
новым направлениям
\
Выход фотонов чер€з боковые грани
/ Источник света
5' \ V
/ Рассеивание из / направления 5' по
направлению 5
Рис. 3. Бесконечно малый элемент объема ткани
Действительно, в левой части уравнения (3.1) первое слагаемое описывает изменение во времени числа фотонов в направлении 5; второе слагаемое соответствует уменьшению числа фотонов вследствие покидания ими рассматриваемого элемента объема через боковые грани. Последнее слагаемое учитывает потерю фотонов в результате поглощения и рассеяния средой. В правой части уравнения переноса излучения первое слагаемое определяет фотоны, рассеиваемые из направления 5' в направлении 5, т. е. фотоны, приходящие в элемент объема; второе слагаемое — не что иное, как функция источника 0.(т,5,1), определяющая число фотонов, генерируемых источником света. Физический смысл уравнения переноса излучения состоит в том, что число фотонов, покинувших элемент объема среды, равняется числу фотонов, пришедших в среду, т. е. данное уравнение описывает закон сохранения энергии.
Главная проблема, с которой имеет дело теория переноса, — определение диффузной составляющей лучевой интенсивности, так как рассеяние фотонов носит случайных характер [8]. Поэтому применяются различные приближения, в соответствии с которыми доминирующим процессом ослабления света является либо поглощение, либо рассеяние [10]. Наиболее часто используемыми являются следующие методы: теория Кубелки—Мунка, диффузионное приближение и метод Монте-Карло. Рассмотрим коротко эти приближения.
3.2. Теория Кубелки—Мунка
Главным допущением данной теории является то, что лучевая интенсивность является диффузной, т. е. Ьр = 0 [10]. Внутри ткани диффузный поток разделен на два: поток Ь1 в направлении падающего излучения и поток, рассеянный назад Ь2, (соответственно, в обратном направлении). Для поглощения и рассеяния диффузного излучения вводятся два коэффициента Кубелки—Мунка, соответственно АКМ и 8КМ. С использованием указанных обозначений можно записать два дифференциальных уравнения:
= -8кмЬ - Акм^1 + 8км^2 > (3.2)
ах йЬ2
= ~8КМ^2 - АКМ^2 + 8КМ^1, (3.3)
ах
где г: определяет среднее направление падающего излучения. Эти уравнения утверждают, что лучевая интенсивность в каждом направлении два раза испытывает потери вследствие поглощения и рассеяния и один раз усиливается вследствие рассеяния фотонов с противоположного направления.
Коэффициенты А и 8 в величинах ца и ц5 записываются следующим образом
[10]: КМ КМ
АКМ = 2М<а, 8КМ = М-5 . (3.4)
Теория Кубелки—Мунка — это частный случай так называемой многопотоковой теории, где уравнение переноса превращается в матричное дифференциальное уравнение, учитывающее лучевую интенсивность в направлении многих отдельных телесных углов. Однако данная теория имеет дело только с диффузной компонентой лучевой интенсивности и ограничена случаями, когда рассеяние во много раз пре-
вышает поглощение. Другим неудобством теории Кубелки—Мунка является то, что она может быть применена только для одномерной геометрии системы.
3.3. Метод Монте-Карло
Численное решение уравнения переноса основывается на методе Монте-Карло. Вообще метод Монте-Карло — это численный метод решения математических задач (систем алгебраических, дифференциальных, интегральных уравнений) [38, 39]. Он позволяет осуществлять прямое статистическое моделирование (физических, химических, биологических, экономических, социальных процессов) при помощи получения и преобразования случайных чисел. Первая работа по использованию метода Монте-Карло была опубликована Холлом [38] в 1873 г. Первая работа, где этот метод излагался систематически, опубликована в 1949 г. Метрополисом и Уламом
[40]. В ней метод Монте-Карло применялся для решения линейных интегральных уравнений. В работе явно угадывалась задача о прохождении нейтронов через вещество. В нашей стране работы по методам Монте-Карло стали активно публиковаться после Международной Женевской конференции по применению атомной энергии в мирных целях. Одной из первых можно считать работу Владимирова и Соболя
[41].
Общая схема метода Монте-Карло основана на центральной предельной теоре-
N
ме теории вероятности, утверждающей, что случайная величина Y = ^ Хг, равная
г=1
сумме большого количества N, произвольных случайных величин Хг с одинаковыми математическими ожиданиями т и дисперсиями а2 всегда распределена по нормальному закону с математическим ожиданием ^т и дисперсией №а2.
Общие свойства методов Монте-Карло: ^
• абсолютная сходимость к решению как —;
• зависимость погрешности от числа испытаний как (для уменьшения
погрешности на порядок необходимо увеличить количество испытаний на два порядка);
• основным методом уменьшения погрешности является максимальное уменьшение дисперсии;
• погрешность не реагирует на размерность задачи. В конечно-разностных методах, при переходе от одномерной задачи к трехмерной, количество вычислений увеличивается на два порядка, в то время как в методах Монте-Карло количество вычислений остается того же порядка;
• простая структура вычислительного алгоритма N раз повторяющиеся однотипные вычисления реализаций случайной величины);
• конструкция случайной величины, вообще говоря, может основываться на физической природе процесса и не требовать обязательной, как в регулярных методах, формулировки уравнения, что для современных проблем становится все более актуальным.
С точки зрения решения уравнения переноса излучения в мутных средах метод Монте-Карло заключается в компьютерном моделировании случайного блуждания N числа фотонов [10]. Для получения приемлемой аппроксимации необходимо рас-
сматривать большое количество фотонов, поскольку точность результатов пропорциональна
Главной идеей метода является учет явлений поглощения и рассеяния на всем оптическом пути фотона через непрозрачную среду. Расстояние между двумя столкновениями выбирается из логарифмического распределения, используя случайное число, генерируемое компьютером. Для учета поглощения каждому фотону присваивается статистический вес, и при распространении через среду этот вес постоянно уменьшается. Если имеет место рассеяние, выбирается новое направление распространения в соответствии с фазовой функцией и другим случайным числом. Эта процедура продолжается до тех пор, пока фотон не выйдет из рассматриваемого объема или его статистический вес не достигнет определенной величины. Метод Монте-Карло включает в себя пять основных шагов: генерация источника фотона, генерация траектории, поглощение, ликвидация, регистрация [10]. Рассмотрим кратко каждый из них.
(I) Генерация источника фотона. Фотоны генерируются на поверхности рассматриваемой среды. Их пространственное и угловое распределение соответствует распределению падающего излучения (например, Гауссов пучок).
(II) Генерация траектории. После генерации фотона определяется расстояние до первого столкновения. Предполагается, что поглощающие и рассеивающие частицы случайно распределены в непрозрачной среде.
Следовательно, величина свободного пробега равна 1/ро5, где р — плотность частиц и о5 — их сечение рассеяния. Случайное число 0 < ^ < 1 генерируется компьютером, и расстояние Ц(^1) до следующего столкновения рассчитывается из выражения
ДУ = (3.5)
1 Р^5
Поскольку 11п Ъ)1йЪ)1 =-1, средняя величина Ц(^1) действительно равна 1/ро5.
0
Отсюда получают рассеивающую точку. Угол рассеяния определяется вторым случайным числом в соответствии с некой фазовой функцией, например функцией Хени—Гринштейна. Соответствующий полярный угол Ф определяется выражением Ф = 2п^3, где — третье случайное число между 0 и 1.
(III) Поглощение. Для учета поглощения каждому фотону присваивается собственный вес. На входе в непрозрачную среду вес фотона равен 1. Вследствие поглощения (в более точных программах также вследствие отражения) вес уменьшается в соответствии с выражением ехр [-|ЯЬ )] . Как альтернатива присвоению веса может быть введено четвертое случайное число между 0 и 1. Затем предполагается, что рассеяние имеет место, только если < а, где а — оптическое альбедо,
которое определяется в соответствии с выражением а = —I5— . Для > а фотон
1а +15
поглощается, что является аналогом шага (IV).
(IV) Ликвидация. Этот шаг используется только в случае присвоения веса каждому фотону в шаге 3. Когда этот вес достигает определенной величины отсечки, фотон ликвидируется. Затем запускается новый фотон, и программа продолжается с шага 1.
(V) Регистрация. После повторения шагов 1-^ для достаточного количества фотонов карта траекторий рассчитывается и накапливается в компьютере. Таким образом может быть получен статистический отчет о порции падающих фотонов, поглощенных средой, а также пространственное и угловое распределение фотонов, вышедших из нее.
Рассмотрим один из вариантов реализации построения алгоритма метода Монте-Карло.
Моделируемая среда задается следующими параметрами: толщиной Ьср, коэффициентами рассеяния ц5 и поглощения ца, средним косинусом угла рассеяния g, относительным показателем преломления п. Среда представляется совокупностью рассеивающих и поглощающих фотоны центров.
Рассмотрим алгоритм (рис. 4). Падающий импульс состоит из одного миллиона фотонов, входящих в среду вдоль оси ^ перпендикулярно ее поверхности (х, у) в точке с координатами (0,0,0). Все расчеты производятся в трехмерной декартовой системе координат. После входа фотона в образец определяются длина свободного пробега фотона в среде и углы рассеяния, 0 и ф. Угол рассеяния 0 задается фазовой функцией рассеяния. В общем случае
р(5,5') = р(9) р(ф), (3.6)
где 5 — направление падения, 5 — направление рассеяния фотона.
Расчет длины свободного пробега и направления движения фотона
1
Фотон поглощен 1
— Фотон достиг
границы слоя
да
Фотон отразился да н __ от границы слоя
(^^Стоп
Рис. 4. Схема алгоритма метода Монте-Карло [10]
Считаем, что частицы среды, на которых происходит рассеяние и поглощение, являются сферически симметричными. Такое приближение часто используется в аналогичных случаях и основано на том, что в процессе прохождения через среду с сильным рассеянием фотон взаимодействует с частицами под разными углами. Поэтому можно применять усредненную индикатрису рассеяния. Использование данной модели и сравнение численных расчетов с экспериментальными результатами показали, что данное приближение удовлетворительно описывает свойства большинства биологических тканей [32-36]. 1
Таким образом, если используется данное приближение, то в нем р(ф) = — .
2п
В случае ткани с сильным рассеянием в качестве фазовой функции рассеяния можно применить фазовую функцию Хени—Гринштейна р(0), откуда получаем выражение для угла 0:
0 = cos
-1
1 + g2 -
1-
g
л2
1 + g -2gRandom
2 g
(3.7)
где Random — случайное равномерно распределенное число из диапазона (0,1).
На каждом шаге угол 0 определяется относительно «старого» направления распространения, угол ф — в плоскости, перпендикулярной «новому» направлению движения.
Длина свободного пробега фотона определяется функцией плотности вероятности:
p(L) =
/ \ L 1
(3.8)
где средняя длина свободного пробега фотона определяется как
lph
1
^ a s
Поскольку
J p(L)dL = 1,
(3.9)
(3.10)
то для расчета длины свободного пробега берется случайное число е (0,1):
Ц
^ = J p(l)dl.
(3.11)
Число равномерно распределенное в интервале (0,1), выдается компьютерным генератором случайных чисел. Таким образом, длина свободного пробега фотона дается выражением
L = -lph ln(1 Ч).
(3.12)
После этого моделируется взаимодействие фотона с частицей среды, которая может быть либо поглощающим, либо рассеивающим центром. Вероятность рассеяния фотона на частице определяется как
ps = , (3.13)
аналогичным образом и вероятность поглощения:
p- =¡^=1 -p- (3Л4)
Если генератор выдает случайное число в диапазоне (0, ps), то считается, что фотон рассеян, в противном случае — поглощен. Весь слой среды вдоль оси z виртуально поделен на некоторое количество более тонких слоев одинаковой толщины, которым соответствуют массивы данных. В каждый из них записывается число поглощенных или рассеянных фотонов. Таким образом, пространственное разрешение по глубине образца составляет 1/Lcp. Если фотон рассеян, то рассчитываются его новое направление движения и координаты по формулам:
х = х0 + L sin 6 cos ф, (3.15)
y = y0 + L sin 6 sin ф, (3.16)
z = z0 + L cos 6. (3.17)
Здесь x0, y0) z0 — «старые» координаты фотона. Если фотон поглощен, то запускается следующий. Потом все координаты пересчитываются в первоначальную систему координат (оси х, y — на поверхности среды, ось z перпендикулярна им и направлена внутрь среды).
Расчет продолжается до тех пор, пока фотон либо в конце концов не будет поглощен, либо не выйдет за границы детектора, либо не попадет на него (z = 0 или z = Lcp). На границах среда—воздух полное внутреннее отражение учитывается при помощи соотношения
hp =sin 1
(3.18)
где n — показатель преломления среды.
3.4. Диффузионное приближение
Данное приближение предполагает, что диффузная интенсивность встречает много частиц и рассеивается на них почти равномерно во всех направлениях, поэтому его угловое распределение почти изотропно [25]. Но угловая зависимость не может сводиться к константе, так как поток при этом обращается в нуль и распространение мощности отсутствует. Поэтому диффузная компонента интенсивности должна быть немного больше для направления полного потока, чем для обратного направления.
Согласно [10] диффузная компонента освещенности может быть представлена в виде сферических гармоник полинома Лежандра. Рассматривая только первые два члена разложения в ряд, мы получим диффузное приближение, которое записывается следующим образом:
13 3
Ь, (г, 5) = — [ Ь, (г, + — [ Ь, (г,5 ')5 '• = Ь0 (г) + — F(r ) • 5, (3.19) 5 4я * 5 4п / 5 4п
где Ь0 (г) — средняя диффузная интенсивность, ¥(г) — вектор диффузного потока, ориентированный вдоль направления единичного вектора 5, [Вт/см2].
Для того чтобы получить точное диффузное уравнение для стационарного случая, необходимо выполнение условия соответствия этого уравнения балансному уравнению для диффузного потока и уравнению, выражающему суть закона сохранения энергии [25].
Первое из этих уравнений выражает закон Фика (плотность потока мощности пропорциональна градиенту освещенности), который описывает уменьшение или увеличение плотности потока мощности за счет поглощения и рассеяния коллими-рованной и диффузной компонент:
¥ (г) = --!- Уф, (г) + ^Е(г, 5о)5о, (3.20)
где = ця + (1 - g)ц5 — транспортный коэффициент затухания. Второе уравнение может быть представлено следующим выражением:
У¥ (г) = -цй ф, (г) + ц ,Е(г, 5о). (3.21)
Физически это уравнение означает, что выходящий из единичного объема поток ¥ равен мощности, излучаемой единицей объема, минус мощность, поглощаемая единицей объема.
Таким образом, в стационарном случае уравнение переноса в диффузионном приближении может быть записано следующим образом [42]:
У2ф, (г) - 3цй цг ф, (г) + 3ц, цг Е(г, 5о) - 3ц ^ У (Е(г, 5о )5о) = 0. (3.22)
Биоткани рассеивают свет преимущественно в направлении вперед. В результате диффузионное приближение не всегда является хорошей аппроксимацией теории переноса излучения вблизи источников или границ. Улучшением ситуации является включение дельта-функции в определение фазовой функции [42]:
р(5,5') = (1 - /) р'(5,5') + /8(1 - 5 • 5') — . (3.23)
2п
Это представление названо приближением Дельта—Эддингтона. Диффузионное уравнение при этом записывается с помощью новых переменных:
^ = Ца + ^5 ,
Ц5 ' = Ц5 (1 - /X
Р '(«,$'),
/ = Я2,
я'=-+г.
я +1
Указанные коэффициенты соответствуют представлению фазовой функции в приближении Хени—Гринштейна [10, 42].
Преобразование р ^ р' (р' — новая фазовая функция) является только математическим преобразованием. Изменения происходят в области источников и границ, что особенно важно для случая сильного рассеяния вперед. В этой ситуации интенсивность характеризуется сильной анизотропией вблизи границ и источников, а это не соответствует описанию интенсивности в диффузионном приближении.
Приближение Дельта—Эддингтона уменьшает степень направленности рассеяния (я' < я). Интенсивность становится менее анизотропной, что приводит к улучшению ситуации вблизи границ и источников.
Заключение
В настоящее время качественная картина распространения света в биотканях может быть описана достаточно полно, что позволяет реализовывать ту или иную стратегию оптической диагностики, терапевтического или хирургического воздействия. Однако в ряде случаев количественные оценки дозы облучения или диагностического параметра оказываются затруднительными из-за отсутствия надежных данных для оптических параметров биотканей или их изменения в процессе взаимодействия света с объектом.
Традиционные спектрофотометрические, угловые и поляризационные методы измерений оказываются полезными для изучения биотканей. Дальнейшее развитие этих методов в направлении биомедицинских приложений требует построения более совершенных моделей тканей, учитывающих пространственное распределение поглотителей и рассеивателей, их полидисперсность, оптическую активность и дву-лучепреломление материала рассеивателей и базового вещества.
Существует необходимость развивать методы решения обратных задач с учетом реальной геометрии образца и лазерного пучка, справедливые для произвольного отношения коэффициентов рассеяния и поглощения. Во многих случаях хорошие результаты дает инверсный метод Монте-Карло для быстрых расчетов, требующихся при реализации практических медицинских методик диагностики и дозиметрии, могут быть полезны методы, базирующиеся на приближенных решениях уравнений теории переноса излучения.
Дозиметрия лазерного излучения — очень важное направление исследований в фотобиологии, необходимое для практической реализации многих методов фототерапии. К настоящему времени многими авторами разработаны основные дозиметрические принципы.
Тем не менее, необходимо развивать 3-мерные дозиметрические модели, учитывающие индивидуальные особенности биотканей, технологию фототерапии и оптические характеристики терапевтических устройств.
Литература
1. Шифрин К. С. Рассеяние света в мутной среде. M.; Л., 1951.
2. Иванов А. П. Оптика рассеивающих сред. Минск, 1969.
3. Тучин В. В. Исследование биотканей методами светорассеяния // Успехи физических наук. 1997. Т. 167, № 5. С. 517-539.
4. Приезжее А. В., Тучин В. В., Шубочкин Л. П. Лазерная диагностика в биологии и медицине. М.: Наука, 1989.
5. Delpy D. Т., СореМ. Quantification in Tissue Near-Infrared Spectroscopy // Phil. Trans. R. Soc. Lond. В 352, 1997. P. 649-659.
6. Аниконов Д. С., Ковтанюк A. E., Прохоров И. В. Использование уравнения переноса в томографии. М.: Логос, 2004.
7. Barabanenkov Yu. N., Kargashin A. Yu. Diffusion calculation of change of back-scattered light beam intensity from turbid medium owing to the existence of an inhomogeneity // J. Mod. Opt. 1993. Vol. 40, N 11. P. 2243-2355.
8. Сетейкин А. Ю., Попов А. П. Взаимодействие света с биологическими тканями и наночастица-ми // LAP Lambert Academic Publishing. 2011. 212 с.
9. Sassaroli А., Fantini S. Comment on the modified Beer-Lambert law for scattering media // Physics in Medicine and Biology. 2004. Vol. 49. P. 255-257.
10. Тучин В. В. Лазеры и волоконная оптика в биомедицинских исследованиях. 2-е изд., испр. и доп. М.: ФИЗМАТЛИТ, 2010. 488 с.
11. Matcher S. J., Cooper C. E. Absolute quantification of deoxy-haemoglobin concentration in tissue near infrared spectroscopy // Physics in Medicine and Biology. 1994. Vol. 39. P. 1295-1312.
12. Bevilacqua F., Berger A. J., Cerussi A. E. et al. Broadband absorption spectroscopy in turbid media by combined frequency-domain and steady-state methods // Appl. Opt. 2000. Vol. 39. P. 6498-6507.
13. Hale G. M., Querry M. R. Optical constants of water in the 200-nm to 200-^m wavelength region // Appl. Opt. 1973. Vol. 12. P. 555-563.
14. Cerussi А., Shah N., Hsiang D. et al. In vivo absorption, scattering, and physiologic properties of 58 malignant breast tumors determined by broad-band diffuse optical spectroscopy // J. Biomed. Opt. 2006. Vol. 11. Р. 044005.
15. Fedorow H., Tribl F., Halliday G. et al. Neuromelanin in human dopamine neurons: comparison with peripheral melanins and relevance to Parkinson's disease // Progress in neurobiology. 2005. Vol. 75. P. 109124.
16. Bowen W. J. The absorption spectra and extinction coefficients of myoglobin // J. Biol. Chem. 1949. Vol. 179. P. 235-245.
17. Balaban W. J., Mootha V. K., Arai A. E. Spectroscopic determination of cytochrome С oxidase content in tissues containing myoglobin or hemoglobin // Anal. Biochem. 1996. Vol. 237. P. 274-278.
18. Prahl S. A. Optical absorption of hemoglobin. 2009. URL http://omlc.ogi.edu/spectra/hemoglobin/ index.html (дата обращения: 12.01.2011).
19. Van Veen R. L. P., Sterenborg H., Pifferi A. et al. Determination of VIS-NIR absorption coefficients of mammalian fat, with time and spatially resolved diffuse reflectance and transmission spectroscopy // Proceedings of Biomedical Topical Meetings. 2000. Paper SF4.
20. Mourant J. R., Fuselier T., Boyer J. et al. Predictions and measurements of scattering and absorption over broad wavelength ranges in tissue phantoms //Applied Optics. 1997. Vol. 36. P. 949-957.
21. Salzman G. C. Light scattering analysis of single cells in Cell Analysis / ed. by N. Catsimpoolas. Plenum. New York, 1982. Vol. 1. P. 111-143.
22. Fujime S., Takasaki-Oshito M., Miyamoto S. Dynamic light scattering from polydisperse suspensions of large spheres // Biophys. J. 1988. Vol. 54. P. 1179-1184.
23. Bigio I. J., Mourant J. R. Ultraviolet and visible spectroscopies for tissue diagnostics: fluorescence spectroscopy and elastic-scattering spectroscopy // Physics in Medicine and Biology. 1997. Vol. 42. P. 803-814.
24. Wang X., Pogue B. W., Jiang S. et al. Approximation of Mie scattering parameters in near-infrared tomography of normal breast tissue in vivo // Journal of Biomedical Optics. 2005. Vol. 10. Р. 051704.
25. Исимару А. Распространение и рассеяние волн в случайно-неоднородных средах. Т. 1: Однократное рассеяние и теория переноса. М.: Мир, 1981. 281 с.
26. Yoon G., Welch A. J., Motamedi M., van Gemert M. J. C. Development and application of threedimen-sional light distribution model for laser irradiated tissue // IEEE J. Quant. Electr. 1987. Vol. QE-23. P. 1721— 1733.
27. GraaffR., Aarnoudse J. G., de Mul F. F. M., Jentink H. W. Light propagation parameters for anisotropi-cally scattering media based on a rigorous solution of the transport equation // Appl. Opt., 1989. Vol. 28. P. 2273-2279.
28. GraaffR., Dassel A. C. M., Koelink M. H. et al. Optical properties of human dermis in vitro and in vivo // Appl. Opt. 1993. Vol. 32. P. 435-447.
29. Koelink M. H., de Mul F. F. M., Greve J. et al. Laser Doppler blood flowmetry using two wavelengths: Monte Carlo simulations and measurements // Appl. Opt. 1994. Vol. 33. P. 3549-3558.
30. Reynolds L. O., McCormick N. J. Approximate two-parameter phase function for light scattering // J. Opt. Soc. Am. 1980. Vol. 70. P. 1206-1212.
31. Kienle A., Patterson M. S., Ott L., Steiner R. Determination of the scattering coefficient and the anisot-ropy factor from laser Doppler spectra of liquids including blood // Appl. Opt. 1996. Vol. 35. P. 3404-3412.
32. Фадеев Д. А., Сетейкин А. Ю. Анализ многократного рассеяния лазерного излучения в биологических средах с пространственными флуктуациями оптических параметров // Научно-технические ведомости СПбГПУ Сер. «Физико-математические науки». 2010. Вып. 2. С. 102-106.
33. Сетейкин А. Ю., Красников И. В., Павлов М. С. Моделирование распространения оптического излучения методом Монте-Карло в биологических средах с замкнутыми внутренними неоднородно-стями // Оптический журнал. 2010. Вып. 77, № 10. С. 15-19.
34. Павлов М. С., Красников И. В., Сетейкин А. Ю. Трехмерная модель распространения в биологических тканях // Научно-технические ведомости СПбГПУ 2008. Т. 67, № 6. С. 120-123.
35. Krasnikov I., Seteikin A. The analysis of the thermal effects arising at interaction of laser radiation with the multilayered biomaterial by using Monte-Carlo method. Fundamental Problems of Optoelectronics and Microelectronics III // Proceedings of SPIE. Vol. 6595 (SPIE, Bellingham, 2007). Р. 47.
36. Красников И. В., Сетейкин А. Ю., Vogel N. Моделирование температурных полей с учетом распространения света в биоткани // Известия Вузов. Приборостроение. 2007. Т. 50, № 9. С. 24-27.
37. Кейз К., Цвайфель П. Линейная теория переноса. М.: Мир, 1972. 383 с.
38. Hall A. On an experiment determination of я // Messeng. Math. 1873. № 2.
39. Самарский А. А. Теория разностных схем. М., 1983.
40. Metropolis N., Ulam S. The Monte-Carlo method // J. Amer. Stat. Assos. 1949. Vol. 44, № 247.
41. Владимиров В. С., Соболь И. М. Расчет наименьшего характеристического числа уравнения Пайерлса методом Монте-Карло // Вычислит. математика. 1958. № 3. С. 130.
42. Star W.M. Diffusion Theory of Light Transport// Optical-Thermal Response of Laser-Irradiated Tissue / еds Welch A. J. and van Gemert M. J. C. New York, 1995. P. 131-206.
Статья поступила в редакцию 15 августа 2013 г.
