МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ РЕГИСТРИРУЕМЫХ СИГНАЛОВ В МЕДИЦИНСКОЙ ЛАЗЕРНОЙ НЕИНВАЗИВНОЙ ФЛЮОРЕСЦЕНТНОЙ ДИАГНОСТИКЕ
УДК 535.243 + 615.471 + 681.7
Математическое моделирование регистрируемых сигналов в медицинской лазерной неинвазивной флюоресцентной диагностике
© 2013 г. Д. А. Рогаткин, доктор техн. наук; О.Д. Смирнова
Московский областной научно-исследовательский клинический институт (МОНИКИ) им. М. Ф. Владимирского, Москва
E-mail: rogatkin@medphyslab.com
На основе модифицированной авторами двухпотоковой модели Кубелки–Мунка, позволяющей в одномерных задачах получать точные аналитические выражения для потоков излучения на границе мутной среды, и решения Кохановского для потока излучения флюоресценции рассматриваются вопросы моделирования спектров вынужденной эндогенной флюоресценции биологических тканей применительно к задачам неинвазивной медицинской диагностики. Приводится аналитическое выражение для функции искажения спектров, зависящей от рассеивающих и поглощающих свойств клеточных биотканей и крови. Показано хорошее совпадение модельных спектров с экспериментальными данными.
Ключевые слова: флюоресценция, медицинская лазерная диагностика, неинвазивная флюоресцентная диагностика, модель Кубелки–Мунка.
Коды OSIC: 170.6280, 170.6510.
Поступила в редакцию 21.01.2013.
Введение
Неинвазивная прижизненная (in vivo, in situ) лазерная флюоресцентная диагностика (ЛФД) появилась на свет как новое научнодиагностическое направление в медицине примерно в середине 1980-х годов [1, 2]. Сегодня многие научные группы по всему миру предпринимают интенсивные попытки вывода ЛФД уже на уровень практического здравоохранения, включая создание систем визуализации внутренних структур биоткани с использованием излучения вынужденной эндогенной флюоресценции (автофлюоресценции) биотканей и построения на этой основе так называемых диффузионных флюоресцентных томографов [3]. Однако ряд проблем, связанных с точностью таких измерений, с количественным определением по интенсивности флюоресценции концентрации активных флюорофоров в зоне обследования, интерпретацией результатов суммарных спектров флюоресценции и т.п., не позволяют пока в полной мере реализовать эти планы [4]. Одной из серьезных проблем в ЛФД является проблема учета искаже-
ний спектров флюоресценции, когда регистрируемый диагностическим прибором in vivo спектр не всегда соответствует исходному спектру флюоресценции эндогенных флюорофоров [5]. Такая ситуация возникает в подавляющем большинстве реальных случаев ввиду наличия в тканях кроме излучающих свет флюорофоров еще и большого количества веществ, интенсивно поглощающих и рассеивающих свет (например, меланина в коже, гемоглобина в крови и т.д.).
Рядом авторов ранее неоднократно предпринимались попытки построения теоретических моделей на основе линейной теории переноса (ТП) для расчета спектров эндогенной флюоресценции биотканей и разработки алгоритмов коррекции спектров с учетом кровенаполнения ткани и некоторых других факторов [см., например, 6–10]. В частности, было показано, что в общем случае регистрируемый in vivo с биоткани сигнал флюоресценции зависит не только от квантового выхода флюоресценции и концентрации флюоресцирующей субстанции в биоткани, но также имеет сложную, причем существенно нелинейную зависимость от общих
54 “Оптический журнал”, 80, 9, 2013
оптико-физических свойств биоткани, формулируемых в терминах ТП – транспортного коэффициента рассеяния (μs), коэффициента поглощения (μa) и т.д. [7, 8]. Однако до последнего времени все подобные работы строились на тех или иных приближенных методах решения транспортной задачи – методе моментов [7], диффузионном приближении [3], упрощенных эвристических алгоритмах [10] или на основе широко известного статистического метода случайных блужданий фотона в среде (метода Монте–Карло). Каждый из этих подходов имеет определенные недостатки. Эвристические алгоритмы, как правило, применимы лишь для выбранной конструкции диагностического устройства. Приближенные решения обладают невысокой точностью, а численные методы, как и метод Монте–Карло, требуют длительных вычислений и не позволяют получать решение в виде замкнутого аналитического выражения, которое могло бы быть легко проанализировано на предмет влияния того или иного параметра на конечный диагностический результат.
Недавно А. Кохановский предложил строгий аналитический метод решения задачи флюоресценции для классического двухпотокового приближения Кубелки–Мунка (КМ) [11]. А в работе [12] авторами был показан основной источник ошибок модели КМ и предложен вариант их устранения, позволяющий в одномерных задачах (1D) получать сегодня точные значения для потоков излучения на границе среды (для регистрируемых прибором потоков). Это открывает перспективы построения простых аналитических алгоритмов реального времени, позволяющих легко анализировать искажения спектров флюоресценции, не прибегая к сложным многоступенчатым вычислениям и численным методам. Цель данной статьи – на основе совмещения модифицированной двухпотоковой модели КМ и решения Кохановского показать в рамках модели 1D возможность получения точного аналитического выражения для функции искажения спектров флюоресценции, а также проанализировать это решение на соответствие типовым экспериментальным данным.
Теоретическая постановка и решение задачи
Рассматривается известная общая постановка фотометрической 1D задачи КМ в среде с флюоресценцией [11], когда при освещении левой границы среды (x = 0) внешним пото-
“Оптический журнал”, 80, 9, 2013
ком излучения Φ0 с длиной волны λ1 внутри среды распространяются вдоль оси x два разнонаправленных потока i(x) и j(x) с исходной длиной волны λ1. Вследствие наличия поглощения и рассеяния света в среде эти потоки частично поглощаются, конвертируются друг в друга, а также вызывают вынужденное излучение флюоресценции находящихся в среде флюорофоров, что приводит к образованию в среде дополнительных аналогичных потоков флюоресценции I(x) и J(x) с длиной волны λ2>λ1 (рис. 1). Никакие волновые свойства излучения в данной модели ввиду ее простейшей фотометрической (энергетической) формулировки и одномерной расчетной схемы во внимание не принимаются [12].
Рис. 1. Постановка задачи распространения излучения в среде с флюоресценцией. Поясне-
ния в тексте.
При диагностических измерениях с передней (освещенной) поверхности биоткани прибором регистрируются выходящие из среды потоки j(0) и J(0), которые и будут нас интересовать в плане получения для них точных аналитических выражений. В такой постановке задачи для потоков i(x), j(x), I(x) и J(x) можно записать две связанные системы линейных дифференциальных уравнений 1-го порядка [11, 12]:
{ diλ1(x) / dx = –β1(λ1)iλ1(x) + β2(λ1)jλ1(x)
djλ1(x) / dx = β1(λ1)jλ1(x) – β2(λ1)iλ1(x)
(1)
{и
dIλ2(x) dJλ2(x)
/ /
dx dx
= =
–ββ1(1λ(λ2)2J)Iλ2λ2((xx)) –+ ββ2(2λ(λ2)2I)Jλ2λ(2x(x) )–+ FF121(2x(x) ,)
(2)
где β1(λ) и β2(λ) называются соответственно коэффициентом ослабления (экстинкции) и коэффициентом обратного рассеяния излучения.
55
Они определяют оптические свойства среды распространения излучения (биоткани). Согласно [12] в модифицированном варианте, в отличие от стандартной модели КМ, коэффициенты β1(λ) и β2(λ) представляют собой до-
статочно сложные функции коэффициентов поглощения μa(λ), плотности неоднородностей в среде μρ и коэффициента отражения Френеля на границах неоднородностей внутри среды R(λ):
β1
=
ω
μа(λ)
–
μρln(1
–
R(λ)) + μρln(1 – ω(λ) + √ ω(λ)2 – R(λ)2e–2μа(λ)
√ ω(λ)2
/ μρ
–
R(λ)2e–2μа(λ)
/
μρ
,
β2
(λ)
=
R(λ)e–μа(λ)
/
μρ
μа(λ)
–
μρln(1
–
R(λ)) + μρln(1 – ω(λ) + √ ω(λ)2 √ ω(λ)2 – R(λ)2e–2μа(λ) / μρ
–
R(λ)2e–2μа(λ)
/
μρ
,
(3)
1 – (1 – R(λ))e–2μа(λ) / μρ
где обозначено ω(λ) =
.
2
Средняя плотность неоднородностей μρ в среде определяется [12] через среднее количество неоднородностей N, на которых происходит рассеяние света, по отношению к толщине рассматриваемого участка биоткани H (рис. 1):
μρ = N / H.
Связь между системами уравнений (1) и (2) осуществляется через функцию F12(x), описывающую формирование излучения флюоресценции внутри среды, –
а решение для jλ1(0) – как
jλ1(0) = Φ0r∞λ1, где в решениях обозначено r∞λ1 = αλ1 = √ β21(λ1) – β22(λ2).
(6)
β2(λ1)
и
β2(λ1) + αλ1
Осуществляя подстановку (5)→(4)→(2) и решая систему (2) относительно Jλ2(x) стандартным методом дифференцирования, получим
для Jλ2(x) следующее дифференциальное уравнение второго порядка:
[ ]1
F12(x)
=
— 2
iλ1(x) + jλ1(x) μaf(λ1)φ12Φ0,
(4)
где μaf(λ1) – коэффициент поглощения исходного излучения активным флюорофором внутри биоткани на длине волны λ1, φ12 – квантовый выход флюоресценции на длине волны λ2 при возбуждении флюоресценции в линии λ1 (коэффициент преобразования излучения λ1→λ2).
Система уравнений (1)–(4) при соответствующих простейших граничных условиях (левая граница (x = 0) iλ1(0) = Φ0 и Iλ2(0) = 0; правая граница (x = H) jλ1(H) = Jλ2(H) = 0) позволяет аналитически решить поставленную задачу относительно потоков jλ1(0) и Jλ2(0).
В случае приближения полубесконечной среды (H→∞) решение системы (1) для суммы потоков iλ1(x) и jλ1(x) может быть записано в виде [11]
iλ1(x) + jλ1(x) = Φ0(1 + r∞λ1)e-αλ1x,
(5)
56
d2Jλ2(x) dx2
–
αλ2
dJλ2(x) dx
= –(β1(λ2) + β2(λ2)) ×
dF(x) × F(x) – dx .
(7)
Как показал Кохановский [11], уравнение (7) легко аналитически решается относительно Jλ2(x). Однако в его решении использовались другие граничные условия, так как в его задаче среда освещалась одновременно и на длине волны λ1, и на λ2, поэтому мы не можем здесь напрямую использовать это решение. В нашей же постановке задачи и в случае потока Jλ2(0), выходящего наружу из биоткани с ее передней поверхности, решение методом Кохановского с указанными выше граничными условиями в приближении полубесконечной среды (H→∞) примет вид
Jλ2(0) = Φ0μaf(λ1)φ12
(1 + r∞λ1) (1 + r∞λ2) 2(αλ1 + αλ2)
.
(8)
“Оптический журнал”, 80, 9, 2013
Произведение Φ0μaf(λ1)φ12 определяет исходный, неискаженный спектр флюоресценции. Таким образом, искажение спектра определяется коэффициентом γ, равным
γ = (1 + r∞λ1)(1 + r∞λ2) , 2(αλ1 + αλ2)
где r∞λ2 =
β2(λ2) β1(λ2) + αλ2
и αλ2 = √ β21(λ2) – β22(λ2).
(9)
Этот коэффициент является функцией длины волны излучения и сложным образом зависит от оптических свойств β1(λ) и β2(λ) среды распространения излучения, так что никакая простейшая нормировка Jλ2(0) на jλ1(0) или Φ0 в (8) не позволяет избавиться от него. Только в идеальном случае «серой» биоткани, если все ее оптические свойства не зависят от длины волны (r∞λi=consti и αλj=constj), для двух длин волн λ2 и λ3 и двух флюорофоров f1 и f2 отношение регистрируемых сигналов флюоресценции будет равно отношению сигналов для неискаженных спектров:
Jλ2(0) = μaf1(λ1)φ12 . Jλ3(0) μaf2(λ1)φ13
(10)
В общем же случае реальных клинических ситуаций для получения исходного (неискаженного) спектра необходимо точное знание γ во всем рабочем диапазоне длин волн, ибо зависимость γ от оптических свойств среды принципиально не линейна и не монотонна.
В качестве иллюстрации этих последних слов на рис. 2 показана зависимость коэффициента γ по (9) от плотности неоднородностей (рассеивателей) μρ внутри среды, а на рис. 3 представлена зависимость γ от отношения μa(λ1)/μρ, также рассчитанная по (9) с использованием (3). Обе зависимости, как мы видим, носят нелинейный характер и далеко не очевидны априори.
Моделирование экспериментальных спектров Преимущества решений (6) и (8) заключа-
ются в их наглядной аналитической форме, которая легко позволяет рассчитывать полный
Рис. 3. Зависимость параметра γ от отношения μa(λ1)/μρ. R(λ1) = 0,03, R(λ2) = 0,02; μρ = 50 мм–1.
Рис. 2. Зависимость параметра γ от плотности рассеивателей. R(λ1) = 0,03; R(λ2) = 0,02, μaλ1 = 4 мм–1.
“Оптический журнал”, 80, 9, 2013
Рис. 4. Экспериментальные спектры эндогенной флюоресценции, снятые in vivo на диагностическом комплексе «ЛАКК–М» при возбуждении флюоресценции узкополосным излучением УФ светодиода с максимумом в линии λ = 365 нм. 1 – кожа хвоста крысы, 2 – слизистая оболочка мягкого неба человека. I – область флюоресценции коллагена, эластина, NADH, кератина и т.д.; II – область флюоресценции флавопротеинов и липофусцина, III – область флюоресценции порфиринов.
57
Рис. 5. Результаты моделирования экспериментальных спектров, показанных на рис. 4. Для порфирина моделировался только один главный максимум протопорфирина IX в рай-
оне λ = 640 нм, StO2 = 95%. Суммарная плотность рассеивателей в ткани μρ = 100 мм–1.
Рис. 6. Зависимость спектров флюоресценции от объемного кровенаполнения ткани в зоне
обследования и степени тканевой сатурации
оксигемоглобина.
Рис. 7. Экспериментально зафиксированное появление инцизуры в спектре обратно рассе-
янного излучения He–Ne-лазера.
Рис. 8. Моделирование инцизуры. Параметры фильтра соответствуют цветному оптическому
стеклу КС-17 по ГОСТ 9411–75 [15].
регистрируемый прибором в процессе диагностики спектр флюоресценции и спектр обратно рассеянного биотканью излучения во всем диапазоне длин волн регистрации и возбуждения флюоресценции. Рисунок 4 демонстрирует такие экспериментально зарегистрированные типовые спектры. Похожие по форме спектры, но с более слабым разрешением по длине волны, приведены и в [13]. В этом смысле они достаточно показательны, поэтому мы специально их выбрали в качестве примера, так как их интерпретация, как мы увидим далее, может быть ошибочна. Спектры рис. 4 зарегистрированы in vivo на диагностическом комплексе «ЛАКК-М» при возбуждении флюоресценции узкополосным излучением УФ светодиода с максимумом в линии λ1 = 365нм [14]. Область флюоресцен-
58
ции I соответствует флюоресценции таких эндогенных флюорофоров в биоткани, как эластин, коллаген, кератин, NADH (никотинамид) и др. В области II наиболее активно флюоресцируют NADH, флавины, липофусцин. Область III наиболее ярко отражает наличие в биоткани и флюоресценцию протопорфирина IX, хлорофилла и ряда других веществ.
Для выравнивания амплитуд обратно рассеянного излучения (левый пик на графиках) и излучения флюоресценции (в более длинноволновой области) в оптической схеме диагностической системы “ЛАКК-М” используются ослабляющие оптические фильтры из цветного оптического стекла [15] марок ОС–13, КС–15 и т.п., ослабляющие регистрируемое обратно рассеянное излучение примерно в 103 раз [14]. По-
“Оптический журнал”, 80, 9, 2013
этому на графиках рис. 4 мы видим уменьшен- функции ослабления оптического фильтра в
ные амплитуды сигналов обратно рассеянного приборе и оптических свойств тестируемого
излучения по сравнению с их реальными зна- объекта, можно на спектре обратно рассеянно-
чениями, и теоретическое моделирование спек- го излучения наблюдать появление инцизуры
тров по (8) проводилось в нашей работе с уче- (выемки, надреза). Пример экспериментально
том спектральной характеристики используе- регистрируемой инцизуры в спектре представ-
мого фильтра. Оптические свойства биоткани лен на рис. 7 для излучения He–Ne-лазера при
при моделировании задавались типовыми зна- использовании в приборе обрезающего оптиче-
чениями μa(λ), μρ и R(λ) с учетом известных спектральных данных по оптическим свойствам ок-
ского фильтра из цветного оптического стекла марки КС-17. На рис. 8 представлен вариант
сигенированного и восстановленного гемогло- теоретической модельной задачи, по результа-
бина крови. Учитывались также разные воз- там вычислений в которой также появляется
можные значения объемного кровенаполнения подобная инцизура. Таким образом, сравнение
Vb биоткани и тканевой сатурации оксигемоглобина StO2 в крови [16]. Сама биоткань представлялась в виде макрооднородной полубесконеч-
экспериментальных и теоретически вычисленных спектров флюоресценции и обратно рассеянного излучения источника возбуждения
ной сплошной светорассеивающей и светопо- флюоресценции позволяет говорить о том, что
глощающей среды-основы, в которой равномер- предложенная модель и алгоритм расчета спек-
но по всему объему распределены в разных кон- тров могут достаточно адекватно описывать ре-
центрациях кровь и различные флюорофоры.
альные ситуации.
Результаты моделирования эксперимен-
тальных спектров, показанных на рис. 4,
представлены на рис. 5. Для порфирина моде-
Заключение
лировался только один главный максимум про-
топорфирина IX в районе λ2 = 640 нм. Тканевая сатурация оксигемоглобина для обоих случаев
В представленной работе, в рамках модифицированной авторами двухпотоковой мо-
выбрана равной StO2 = 95%. Суммарная плотность рассеивателей в ткани μρ = 100 мм–1. R(λ) рассчитывался для каждой длины волны как
дели Кубелки–Мунка и решения Кохановского для излучения флюоресценции, предложено аналитически строгое решение задачи теоре-
коэффициент отражения Френеля на границе тического моделирования спектров вынужден-
вода/воздух. Оба рисунка демонстрируют весь- ной эндогенной флюоресценции биологических
ма хорошее как качественное, так и количе- тканей применительно к задачам неинвазив-
ственное совпадение спектров при выбранных ной медицинской диагностики. Продемонстри-
значениях Vb. Важно отметить, что небольшое видимое увеличение амплитуды в спектре 2 на
ровано адекватное описание экспериментально зарегистрированных in vivo спектров флю-
рис. 4 в районе λ2 = 604 нм можно трактовать как повышенное содержание в ткани флюо-
оресценции реальных биотканей. Полученное решение (9) для функции искажения спектров
ресцирующих флавинов или липофусцина. И в некоторых работах иногда встречается такая
показывает, что эта функция (коэффициент γ) принципиально нелинейно и немонотонно за-
трактовка. Однако это «увеличение» – лишь висит от оптических свойств среды распростра-
следствие поглощения излучения флюоресцен- нения излучения. Поэтому в общем случае ре-
ции кровью. На рис. 6 показаны модельные альных клинических диагностических обследо-
спектры флюоресценции биоткани для разных ваний для получения исходного (неискаженно-
значений Vb и StO2 крови при неизменном содержании флюорофоров. С увеличением кровена-
го) спектра необходимо точное знание γ во всем рабочем диапазоне длин волн для каждого па-
полнения и тканевой сатурации оксигемогло- циента и каждой его области исследования. Со-
бина в крови как раз и проявляется искажение истинного спектра в виде общего нелинейного
ответственно, на методы определения γ и их погрешности и должны быть направлены сегодня
уменьшения амплитуды регистрируемых сиг- все исследования, ориентированные на реше-
налов и появления характерных гемоглобино- ние задачи количественной обработки спектров
вых «провалов» в районе λ2=530–580 нм. Иногда, при определенных сочетаниях
флюоресценции в практической медицине. Работа выполнена в рамках гранта РФФИ №
спектра излучения источника, спектральной 13–08–12188.
*****
“Оптический журнал”, 80, 9, 2013
59
ЛИТЕРАТУРА 1. Loschenov V.B., Konov V.I., Prokhorov A.M. Photodynamic Therapy and Fluorescence Diagnostics // Laser
Physics. 2000.V. 10. №. 6. P. 1188–1207. 2. Handbook of biomedical fluorescence / Ed. by Mycek M.A., Pogue B.W. New York: Marcel Dekker Inc., 2003.
665 p. 3. Leblond F., Davis S., Valdes P., Pouge B. Pre-clinical whole-body fluorescence imaging: review of instruments,
methods and applications // J. Photochem. Photobiol. 2010. V. 98. P. 77–94. 4. Bonfert-Taylor P., Leblond F., Holt R., Tichauer K., Pouge B., Taylor E. Information loss and reconstruction in
diffuse fluorescence tomography // J. Opt. Soc. Am. A. 2012. V. 29. №. 3. P. 321–330. 5. Рогаткин Д.А. Инструментальные и методические погрешности измерений в неинвазивной медицин-
ской спектрофотометрии // Материалы III Евразийского конгресса по медицинской физике и инженерии «Медицинская физика-2010». Т. III. М.: МГУ, 2010. С. 38–41. 6. Putten W.J. M., Gemert M.J.C. A modeling approach to the detection of subcutaneous tumours by haematoporphyrin-derivative fluorescence // Phys. Med. Biol. 1983. V. 28. №. 6. P. 639–645. 7. Rogatkin D., Svirin V., Hachaturyan G. The theoretical model for fluorescent field calculation in nonhomogenious and scattering biological tissues // Proc. SPIE. 1998. V. 3563. P. 125–136. 8. Rogatkin D.A., Tchernyi V.V. Mathematical simulation as a key point of the laser fluorescence diagnostic technique in oncology // Proc. SPIE. 2000. V. 4059. P. 73–78. 9. Baraghis E., Devor A., Fang Q., Srinivasan V., Wu W., Boas D., Sakadzic S., Lesage F., Ayata C., Kasischke K. Two-photon microscopy of cortical NADH fluorescence intensity changes: correcting contamination from the hemodynamic response // J. Biomed. Opt. 2011. V. 16. 106003. 10. Kanick S., Robinson D., Sterenborg C., Amelink A. Extraction of intrinsic from single fiber fluorescence measurements on a turbid medium // Optics Letters. 2012. V. 37. №.5. P. 948–950. 11. Kokhanovsky A.A. Radiative properties of optically thick fluorescent turbid media // J. Opt. Soc. Am. A. 2009. V. 26. №. 8. P. 1896–1900. 12. Рогаткин Д.А. Об особенности в определении оптических свойств мутных биологических тканей и сред в расчетных задачах медицинской неинвазивной спектрофотометрии // Медицинская техника. 2007. № 2. С. 10–16. 13. Субботин А.Р., Савельева Т.А., Горяйнов С.А. Алгоритм обработки спектров флуоресценции протопорфирина IX и эндогенных флуорофоров в глиальных опухолях головного мозга // Сб. материалов V Троицкой конф. «Медицинская физика и инновации в медицине». 2012. Т. 1. Троицк. С.268–269. 14. Rogatkin D.A., Lapaeva L.G., Petritskaya E.N., Sidorov V.V., Shumskiy V.I. Multifunctional laser noninvasive spectroscopic system for medical diagnostics and some metrological provisions for that // Proc. SPIE. 2009. V. 7368. 73681Y. 15. Стекло оптическое цветное. ГОСТ 9411–75. М.: Изд-во Стандартов. 1976. 50 с. 16. Рогаткин Д.А. Физические основы оптической оксиметрии // Медицинская физика. 2012. №2. С. 97– 114.
60 “Оптический журнал”, 80, 9, 2013
Математическое моделирование регистрируемых сигналов в медицинской лазерной неинвазивной флюоресцентной диагностике
© 2013 г. Д. А. Рогаткин, доктор техн. наук; О.Д. Смирнова
Московский областной научно-исследовательский клинический институт (МОНИКИ) им. М. Ф. Владимирского, Москва
E-mail: rogatkin@medphyslab.com
На основе модифицированной авторами двухпотоковой модели Кубелки–Мунка, позволяющей в одномерных задачах получать точные аналитические выражения для потоков излучения на границе мутной среды, и решения Кохановского для потока излучения флюоресценции рассматриваются вопросы моделирования спектров вынужденной эндогенной флюоресценции биологических тканей применительно к задачам неинвазивной медицинской диагностики. Приводится аналитическое выражение для функции искажения спектров, зависящей от рассеивающих и поглощающих свойств клеточных биотканей и крови. Показано хорошее совпадение модельных спектров с экспериментальными данными.
Ключевые слова: флюоресценция, медицинская лазерная диагностика, неинвазивная флюоресцентная диагностика, модель Кубелки–Мунка.
Коды OSIC: 170.6280, 170.6510.
Поступила в редакцию 21.01.2013.
Введение
Неинвазивная прижизненная (in vivo, in situ) лазерная флюоресцентная диагностика (ЛФД) появилась на свет как новое научнодиагностическое направление в медицине примерно в середине 1980-х годов [1, 2]. Сегодня многие научные группы по всему миру предпринимают интенсивные попытки вывода ЛФД уже на уровень практического здравоохранения, включая создание систем визуализации внутренних структур биоткани с использованием излучения вынужденной эндогенной флюоресценции (автофлюоресценции) биотканей и построения на этой основе так называемых диффузионных флюоресцентных томографов [3]. Однако ряд проблем, связанных с точностью таких измерений, с количественным определением по интенсивности флюоресценции концентрации активных флюорофоров в зоне обследования, интерпретацией результатов суммарных спектров флюоресценции и т.п., не позволяют пока в полной мере реализовать эти планы [4]. Одной из серьезных проблем в ЛФД является проблема учета искаже-
ний спектров флюоресценции, когда регистрируемый диагностическим прибором in vivo спектр не всегда соответствует исходному спектру флюоресценции эндогенных флюорофоров [5]. Такая ситуация возникает в подавляющем большинстве реальных случаев ввиду наличия в тканях кроме излучающих свет флюорофоров еще и большого количества веществ, интенсивно поглощающих и рассеивающих свет (например, меланина в коже, гемоглобина в крови и т.д.).
Рядом авторов ранее неоднократно предпринимались попытки построения теоретических моделей на основе линейной теории переноса (ТП) для расчета спектров эндогенной флюоресценции биотканей и разработки алгоритмов коррекции спектров с учетом кровенаполнения ткани и некоторых других факторов [см., например, 6–10]. В частности, было показано, что в общем случае регистрируемый in vivo с биоткани сигнал флюоресценции зависит не только от квантового выхода флюоресценции и концентрации флюоресцирующей субстанции в биоткани, но также имеет сложную, причем существенно нелинейную зависимость от общих
54 “Оптический журнал”, 80, 9, 2013
оптико-физических свойств биоткани, формулируемых в терминах ТП – транспортного коэффициента рассеяния (μs), коэффициента поглощения (μa) и т.д. [7, 8]. Однако до последнего времени все подобные работы строились на тех или иных приближенных методах решения транспортной задачи – методе моментов [7], диффузионном приближении [3], упрощенных эвристических алгоритмах [10] или на основе широко известного статистического метода случайных блужданий фотона в среде (метода Монте–Карло). Каждый из этих подходов имеет определенные недостатки. Эвристические алгоритмы, как правило, применимы лишь для выбранной конструкции диагностического устройства. Приближенные решения обладают невысокой точностью, а численные методы, как и метод Монте–Карло, требуют длительных вычислений и не позволяют получать решение в виде замкнутого аналитического выражения, которое могло бы быть легко проанализировано на предмет влияния того или иного параметра на конечный диагностический результат.
Недавно А. Кохановский предложил строгий аналитический метод решения задачи флюоресценции для классического двухпотокового приближения Кубелки–Мунка (КМ) [11]. А в работе [12] авторами был показан основной источник ошибок модели КМ и предложен вариант их устранения, позволяющий в одномерных задачах (1D) получать сегодня точные значения для потоков излучения на границе среды (для регистрируемых прибором потоков). Это открывает перспективы построения простых аналитических алгоритмов реального времени, позволяющих легко анализировать искажения спектров флюоресценции, не прибегая к сложным многоступенчатым вычислениям и численным методам. Цель данной статьи – на основе совмещения модифицированной двухпотоковой модели КМ и решения Кохановского показать в рамках модели 1D возможность получения точного аналитического выражения для функции искажения спектров флюоресценции, а также проанализировать это решение на соответствие типовым экспериментальным данным.
Теоретическая постановка и решение задачи
Рассматривается известная общая постановка фотометрической 1D задачи КМ в среде с флюоресценцией [11], когда при освещении левой границы среды (x = 0) внешним пото-
“Оптический журнал”, 80, 9, 2013
ком излучения Φ0 с длиной волны λ1 внутри среды распространяются вдоль оси x два разнонаправленных потока i(x) и j(x) с исходной длиной волны λ1. Вследствие наличия поглощения и рассеяния света в среде эти потоки частично поглощаются, конвертируются друг в друга, а также вызывают вынужденное излучение флюоресценции находящихся в среде флюорофоров, что приводит к образованию в среде дополнительных аналогичных потоков флюоресценции I(x) и J(x) с длиной волны λ2>λ1 (рис. 1). Никакие волновые свойства излучения в данной модели ввиду ее простейшей фотометрической (энергетической) формулировки и одномерной расчетной схемы во внимание не принимаются [12].
Рис. 1. Постановка задачи распространения излучения в среде с флюоресценцией. Поясне-
ния в тексте.
При диагностических измерениях с передней (освещенной) поверхности биоткани прибором регистрируются выходящие из среды потоки j(0) и J(0), которые и будут нас интересовать в плане получения для них точных аналитических выражений. В такой постановке задачи для потоков i(x), j(x), I(x) и J(x) можно записать две связанные системы линейных дифференциальных уравнений 1-го порядка [11, 12]:
{ diλ1(x) / dx = –β1(λ1)iλ1(x) + β2(λ1)jλ1(x)
djλ1(x) / dx = β1(λ1)jλ1(x) – β2(λ1)iλ1(x)
(1)
{и
dIλ2(x) dJλ2(x)
/ /
dx dx
= =
–ββ1(1λ(λ2)2J)Iλ2λ2((xx)) –+ ββ2(2λ(λ2)2I)Jλ2λ(2x(x) )–+ FF121(2x(x) ,)
(2)
где β1(λ) и β2(λ) называются соответственно коэффициентом ослабления (экстинкции) и коэффициентом обратного рассеяния излучения.
55
Они определяют оптические свойства среды распространения излучения (биоткани). Согласно [12] в модифицированном варианте, в отличие от стандартной модели КМ, коэффициенты β1(λ) и β2(λ) представляют собой до-
статочно сложные функции коэффициентов поглощения μa(λ), плотности неоднородностей в среде μρ и коэффициента отражения Френеля на границах неоднородностей внутри среды R(λ):
β1
=
ω
μа(λ)
–
μρln(1
–
R(λ)) + μρln(1 – ω(λ) + √ ω(λ)2 – R(λ)2e–2μа(λ)
√ ω(λ)2
/ μρ
–
R(λ)2e–2μа(λ)
/
μρ
,
β2
(λ)
=
R(λ)e–μа(λ)
/
μρ
μа(λ)
–
μρln(1
–
R(λ)) + μρln(1 – ω(λ) + √ ω(λ)2 √ ω(λ)2 – R(λ)2e–2μа(λ) / μρ
–
R(λ)2e–2μа(λ)
/
μρ
,
(3)
1 – (1 – R(λ))e–2μа(λ) / μρ
где обозначено ω(λ) =
.
2
Средняя плотность неоднородностей μρ в среде определяется [12] через среднее количество неоднородностей N, на которых происходит рассеяние света, по отношению к толщине рассматриваемого участка биоткани H (рис. 1):
μρ = N / H.
Связь между системами уравнений (1) и (2) осуществляется через функцию F12(x), описывающую формирование излучения флюоресценции внутри среды, –
а решение для jλ1(0) – как
jλ1(0) = Φ0r∞λ1, где в решениях обозначено r∞λ1 = αλ1 = √ β21(λ1) – β22(λ2).
(6)
β2(λ1)
и
β2(λ1) + αλ1
Осуществляя подстановку (5)→(4)→(2) и решая систему (2) относительно Jλ2(x) стандартным методом дифференцирования, получим
для Jλ2(x) следующее дифференциальное уравнение второго порядка:
[ ]1
F12(x)
=
— 2
iλ1(x) + jλ1(x) μaf(λ1)φ12Φ0,
(4)
где μaf(λ1) – коэффициент поглощения исходного излучения активным флюорофором внутри биоткани на длине волны λ1, φ12 – квантовый выход флюоресценции на длине волны λ2 при возбуждении флюоресценции в линии λ1 (коэффициент преобразования излучения λ1→λ2).
Система уравнений (1)–(4) при соответствующих простейших граничных условиях (левая граница (x = 0) iλ1(0) = Φ0 и Iλ2(0) = 0; правая граница (x = H) jλ1(H) = Jλ2(H) = 0) позволяет аналитически решить поставленную задачу относительно потоков jλ1(0) и Jλ2(0).
В случае приближения полубесконечной среды (H→∞) решение системы (1) для суммы потоков iλ1(x) и jλ1(x) может быть записано в виде [11]
iλ1(x) + jλ1(x) = Φ0(1 + r∞λ1)e-αλ1x,
(5)
56
d2Jλ2(x) dx2
–
αλ2
dJλ2(x) dx
= –(β1(λ2) + β2(λ2)) ×
dF(x) × F(x) – dx .
(7)
Как показал Кохановский [11], уравнение (7) легко аналитически решается относительно Jλ2(x). Однако в его решении использовались другие граничные условия, так как в его задаче среда освещалась одновременно и на длине волны λ1, и на λ2, поэтому мы не можем здесь напрямую использовать это решение. В нашей же постановке задачи и в случае потока Jλ2(0), выходящего наружу из биоткани с ее передней поверхности, решение методом Кохановского с указанными выше граничными условиями в приближении полубесконечной среды (H→∞) примет вид
Jλ2(0) = Φ0μaf(λ1)φ12
(1 + r∞λ1) (1 + r∞λ2) 2(αλ1 + αλ2)
.
(8)
“Оптический журнал”, 80, 9, 2013
Произведение Φ0μaf(λ1)φ12 определяет исходный, неискаженный спектр флюоресценции. Таким образом, искажение спектра определяется коэффициентом γ, равным
γ = (1 + r∞λ1)(1 + r∞λ2) , 2(αλ1 + αλ2)
где r∞λ2 =
β2(λ2) β1(λ2) + αλ2
и αλ2 = √ β21(λ2) – β22(λ2).
(9)
Этот коэффициент является функцией длины волны излучения и сложным образом зависит от оптических свойств β1(λ) и β2(λ) среды распространения излучения, так что никакая простейшая нормировка Jλ2(0) на jλ1(0) или Φ0 в (8) не позволяет избавиться от него. Только в идеальном случае «серой» биоткани, если все ее оптические свойства не зависят от длины волны (r∞λi=consti и αλj=constj), для двух длин волн λ2 и λ3 и двух флюорофоров f1 и f2 отношение регистрируемых сигналов флюоресценции будет равно отношению сигналов для неискаженных спектров:
Jλ2(0) = μaf1(λ1)φ12 . Jλ3(0) μaf2(λ1)φ13
(10)
В общем же случае реальных клинических ситуаций для получения исходного (неискаженного) спектра необходимо точное знание γ во всем рабочем диапазоне длин волн, ибо зависимость γ от оптических свойств среды принципиально не линейна и не монотонна.
В качестве иллюстрации этих последних слов на рис. 2 показана зависимость коэффициента γ по (9) от плотности неоднородностей (рассеивателей) μρ внутри среды, а на рис. 3 представлена зависимость γ от отношения μa(λ1)/μρ, также рассчитанная по (9) с использованием (3). Обе зависимости, как мы видим, носят нелинейный характер и далеко не очевидны априори.
Моделирование экспериментальных спектров Преимущества решений (6) и (8) заключа-
ются в их наглядной аналитической форме, которая легко позволяет рассчитывать полный
Рис. 3. Зависимость параметра γ от отношения μa(λ1)/μρ. R(λ1) = 0,03, R(λ2) = 0,02; μρ = 50 мм–1.
Рис. 2. Зависимость параметра γ от плотности рассеивателей. R(λ1) = 0,03; R(λ2) = 0,02, μaλ1 = 4 мм–1.
“Оптический журнал”, 80, 9, 2013
Рис. 4. Экспериментальные спектры эндогенной флюоресценции, снятые in vivo на диагностическом комплексе «ЛАКК–М» при возбуждении флюоресценции узкополосным излучением УФ светодиода с максимумом в линии λ = 365 нм. 1 – кожа хвоста крысы, 2 – слизистая оболочка мягкого неба человека. I – область флюоресценции коллагена, эластина, NADH, кератина и т.д.; II – область флюоресценции флавопротеинов и липофусцина, III – область флюоресценции порфиринов.
57
Рис. 5. Результаты моделирования экспериментальных спектров, показанных на рис. 4. Для порфирина моделировался только один главный максимум протопорфирина IX в рай-
оне λ = 640 нм, StO2 = 95%. Суммарная плотность рассеивателей в ткани μρ = 100 мм–1.
Рис. 6. Зависимость спектров флюоресценции от объемного кровенаполнения ткани в зоне
обследования и степени тканевой сатурации
оксигемоглобина.
Рис. 7. Экспериментально зафиксированное появление инцизуры в спектре обратно рассе-
янного излучения He–Ne-лазера.
Рис. 8. Моделирование инцизуры. Параметры фильтра соответствуют цветному оптическому
стеклу КС-17 по ГОСТ 9411–75 [15].
регистрируемый прибором в процессе диагностики спектр флюоресценции и спектр обратно рассеянного биотканью излучения во всем диапазоне длин волн регистрации и возбуждения флюоресценции. Рисунок 4 демонстрирует такие экспериментально зарегистрированные типовые спектры. Похожие по форме спектры, но с более слабым разрешением по длине волны, приведены и в [13]. В этом смысле они достаточно показательны, поэтому мы специально их выбрали в качестве примера, так как их интерпретация, как мы увидим далее, может быть ошибочна. Спектры рис. 4 зарегистрированы in vivo на диагностическом комплексе «ЛАКК-М» при возбуждении флюоресценции узкополосным излучением УФ светодиода с максимумом в линии λ1 = 365нм [14]. Область флюоресцен-
58
ции I соответствует флюоресценции таких эндогенных флюорофоров в биоткани, как эластин, коллаген, кератин, NADH (никотинамид) и др. В области II наиболее активно флюоресцируют NADH, флавины, липофусцин. Область III наиболее ярко отражает наличие в биоткани и флюоресценцию протопорфирина IX, хлорофилла и ряда других веществ.
Для выравнивания амплитуд обратно рассеянного излучения (левый пик на графиках) и излучения флюоресценции (в более длинноволновой области) в оптической схеме диагностической системы “ЛАКК-М” используются ослабляющие оптические фильтры из цветного оптического стекла [15] марок ОС–13, КС–15 и т.п., ослабляющие регистрируемое обратно рассеянное излучение примерно в 103 раз [14]. По-
“Оптический журнал”, 80, 9, 2013
этому на графиках рис. 4 мы видим уменьшен- функции ослабления оптического фильтра в
ные амплитуды сигналов обратно рассеянного приборе и оптических свойств тестируемого
излучения по сравнению с их реальными зна- объекта, можно на спектре обратно рассеянно-
чениями, и теоретическое моделирование спек- го излучения наблюдать появление инцизуры
тров по (8) проводилось в нашей работе с уче- (выемки, надреза). Пример экспериментально
том спектральной характеристики используе- регистрируемой инцизуры в спектре представ-
мого фильтра. Оптические свойства биоткани лен на рис. 7 для излучения He–Ne-лазера при
при моделировании задавались типовыми зна- использовании в приборе обрезающего оптиче-
чениями μa(λ), μρ и R(λ) с учетом известных спектральных данных по оптическим свойствам ок-
ского фильтра из цветного оптического стекла марки КС-17. На рис. 8 представлен вариант
сигенированного и восстановленного гемогло- теоретической модельной задачи, по результа-
бина крови. Учитывались также разные воз- там вычислений в которой также появляется
можные значения объемного кровенаполнения подобная инцизура. Таким образом, сравнение
Vb биоткани и тканевой сатурации оксигемоглобина StO2 в крови [16]. Сама биоткань представлялась в виде макрооднородной полубесконеч-
экспериментальных и теоретически вычисленных спектров флюоресценции и обратно рассеянного излучения источника возбуждения
ной сплошной светорассеивающей и светопо- флюоресценции позволяет говорить о том, что
глощающей среды-основы, в которой равномер- предложенная модель и алгоритм расчета спек-
но по всему объему распределены в разных кон- тров могут достаточно адекватно описывать ре-
центрациях кровь и различные флюорофоры.
альные ситуации.
Результаты моделирования эксперимен-
тальных спектров, показанных на рис. 4,
представлены на рис. 5. Для порфирина моде-
Заключение
лировался только один главный максимум про-
топорфирина IX в районе λ2 = 640 нм. Тканевая сатурация оксигемоглобина для обоих случаев
В представленной работе, в рамках модифицированной авторами двухпотоковой мо-
выбрана равной StO2 = 95%. Суммарная плотность рассеивателей в ткани μρ = 100 мм–1. R(λ) рассчитывался для каждой длины волны как
дели Кубелки–Мунка и решения Кохановского для излучения флюоресценции, предложено аналитически строгое решение задачи теоре-
коэффициент отражения Френеля на границе тического моделирования спектров вынужден-
вода/воздух. Оба рисунка демонстрируют весь- ной эндогенной флюоресценции биологических
ма хорошее как качественное, так и количе- тканей применительно к задачам неинвазив-
ственное совпадение спектров при выбранных ной медицинской диагностики. Продемонстри-
значениях Vb. Важно отметить, что небольшое видимое увеличение амплитуды в спектре 2 на
ровано адекватное описание экспериментально зарегистрированных in vivo спектров флю-
рис. 4 в районе λ2 = 604 нм можно трактовать как повышенное содержание в ткани флюо-
оресценции реальных биотканей. Полученное решение (9) для функции искажения спектров
ресцирующих флавинов или липофусцина. И в некоторых работах иногда встречается такая
показывает, что эта функция (коэффициент γ) принципиально нелинейно и немонотонно за-
трактовка. Однако это «увеличение» – лишь висит от оптических свойств среды распростра-
следствие поглощения излучения флюоресцен- нения излучения. Поэтому в общем случае ре-
ции кровью. На рис. 6 показаны модельные альных клинических диагностических обследо-
спектры флюоресценции биоткани для разных ваний для получения исходного (неискаженно-
значений Vb и StO2 крови при неизменном содержании флюорофоров. С увеличением кровена-
го) спектра необходимо точное знание γ во всем рабочем диапазоне длин волн для каждого па-
полнения и тканевой сатурации оксигемогло- циента и каждой его области исследования. Со-
бина в крови как раз и проявляется искажение истинного спектра в виде общего нелинейного
ответственно, на методы определения γ и их погрешности и должны быть направлены сегодня
уменьшения амплитуды регистрируемых сиг- все исследования, ориентированные на реше-
налов и появления характерных гемоглобино- ние задачи количественной обработки спектров
вых «провалов» в районе λ2=530–580 нм. Иногда, при определенных сочетаниях
флюоресценции в практической медицине. Работа выполнена в рамках гранта РФФИ №
спектра излучения источника, спектральной 13–08–12188.
*****
“Оптический журнал”, 80, 9, 2013
59
ЛИТЕРАТУРА 1. Loschenov V.B., Konov V.I., Prokhorov A.M. Photodynamic Therapy and Fluorescence Diagnostics // Laser
Physics. 2000.V. 10. №. 6. P. 1188–1207. 2. Handbook of biomedical fluorescence / Ed. by Mycek M.A., Pogue B.W. New York: Marcel Dekker Inc., 2003.
665 p. 3. Leblond F., Davis S., Valdes P., Pouge B. Pre-clinical whole-body fluorescence imaging: review of instruments,
methods and applications // J. Photochem. Photobiol. 2010. V. 98. P. 77–94. 4. Bonfert-Taylor P., Leblond F., Holt R., Tichauer K., Pouge B., Taylor E. Information loss and reconstruction in
diffuse fluorescence tomography // J. Opt. Soc. Am. A. 2012. V. 29. №. 3. P. 321–330. 5. Рогаткин Д.А. Инструментальные и методические погрешности измерений в неинвазивной медицин-
ской спектрофотометрии // Материалы III Евразийского конгресса по медицинской физике и инженерии «Медицинская физика-2010». Т. III. М.: МГУ, 2010. С. 38–41. 6. Putten W.J. M., Gemert M.J.C. A modeling approach to the detection of subcutaneous tumours by haematoporphyrin-derivative fluorescence // Phys. Med. Biol. 1983. V. 28. №. 6. P. 639–645. 7. Rogatkin D., Svirin V., Hachaturyan G. The theoretical model for fluorescent field calculation in nonhomogenious and scattering biological tissues // Proc. SPIE. 1998. V. 3563. P. 125–136. 8. Rogatkin D.A., Tchernyi V.V. Mathematical simulation as a key point of the laser fluorescence diagnostic technique in oncology // Proc. SPIE. 2000. V. 4059. P. 73–78. 9. Baraghis E., Devor A., Fang Q., Srinivasan V., Wu W., Boas D., Sakadzic S., Lesage F., Ayata C., Kasischke K. Two-photon microscopy of cortical NADH fluorescence intensity changes: correcting contamination from the hemodynamic response // J. Biomed. Opt. 2011. V. 16. 106003. 10. Kanick S., Robinson D., Sterenborg C., Amelink A. Extraction of intrinsic from single fiber fluorescence measurements on a turbid medium // Optics Letters. 2012. V. 37. №.5. P. 948–950. 11. Kokhanovsky A.A. Radiative properties of optically thick fluorescent turbid media // J. Opt. Soc. Am. A. 2009. V. 26. №. 8. P. 1896–1900. 12. Рогаткин Д.А. Об особенности в определении оптических свойств мутных биологических тканей и сред в расчетных задачах медицинской неинвазивной спектрофотометрии // Медицинская техника. 2007. № 2. С. 10–16. 13. Субботин А.Р., Савельева Т.А., Горяйнов С.А. Алгоритм обработки спектров флуоресценции протопорфирина IX и эндогенных флуорофоров в глиальных опухолях головного мозга // Сб. материалов V Троицкой конф. «Медицинская физика и инновации в медицине». 2012. Т. 1. Троицк. С.268–269. 14. Rogatkin D.A., Lapaeva L.G., Petritskaya E.N., Sidorov V.V., Shumskiy V.I. Multifunctional laser noninvasive spectroscopic system for medical diagnostics and some metrological provisions for that // Proc. SPIE. 2009. V. 7368. 73681Y. 15. Стекло оптическое цветное. ГОСТ 9411–75. М.: Изд-во Стандартов. 1976. 50 с. 16. Рогаткин Д.А. Физические основы оптической оксиметрии // Медицинская физика. 2012. №2. С. 97– 114.
60 “Оптический журнал”, 80, 9, 2013