CC BY
5
Основанное на методе Монте-Карло моделирование временных функций рассеяния точки и функций чувствительности для мезоскопической время-разрешенной флуоресцентной молекулярной томографии
С.И. Самарин1, А.Б. Коновалов1, В.В. Власов1, И.Д. Соловьев 2, А.П. Савицкий 2, В.В. Тучин 2,3 1ФГУП "Российский федеральный ядерный центр - ВНИИ технической физики им. академика Е.И. Забабахина ", 456770, Россия, Челябинская обл., Снежинск, ул. Васильева, д. 13;
2 ФГУ "Федеральный исследовательский центр "Фундаментальные основы биотехнологии РАН," 119071, Россия, Москва, Ленинский проспект, д. 33, зд. 2;
3 Саратовский национальный исследовательский государственный университет им. Н.Г. Чернышевского,
410012, Россия, Саратов, ул. Астраханская, д. 83
Аннотация
В статье описан алгоритм программы TurbidMC, реализующей пертурбационный метод Монте-Карло и предназначенной для моделирования временных функций рассеяния точки и функций чувствительности для задач время-разрешенной флуоресцентной молекулярной томографии (FMT). Программа ориентирована на работу с конкретным ранее опубликованным методом FMT ([22] в списке литературы), определяющим специфику расчета функций чувствительности. Согласно этому методу обратная задача изначально решается относительно обобщенной функции распределения параметров флуоресценции, а затем уже выполняется разделение распределений коэффициента поглощения флуорофора и времени жизни флуоресценции. Корректность работы программы проверена сравнением тестовых расчетов флуоресцентных временных функций рассеяния точки с данными эксперимента по сканированию фантома с флуорофором трехканальным зондом в мезоскопическом режиме обратного рассеяния. Также приведен пример восстановления распределений параметров флуоресценции, подтверждающий корректность расчетов функций чувствительности.
Ключевые слова: программа TurbidMC, метод Монте-Карло, флуоресцентная молекулярная томография, временные функции рассеяния точки, функции чувствительности, коэффициент поглощения флуорофора, время жизни флуоресценции.
Цитирование: Самарин, С.И. Основанное на методе Монте-Карло моделирование временных функций рассеяния точки и функций чувствительности для мезоскопической время-разрешенной флуоресцентной молекулярной томографии / С.И. Самарин, А.Б. Коновалов, В.В. Власов, И.Д. Соловьев, А.П. Савицкий, В.В. Тучин // Компьютерная оптика. - 2023. -Т. 47, № 5. - С. 673-690. - DOI: 10.18287/2412-6179-CO-1295.
Citation: Samarin SI, Konovalov AB, Vlasov VV, Solovyev ID, Savitsky AP, Tuchin VV. Monte Carlo modeling of temporal point spread functions and sensitivity functions for mesoscopic time-resolved fluorescence molecular tomography. Computer Optics 2023; 47(5): 673-690. DOI: 10.18287/2412-6179-CO-1295.
Введение
В последние годы метод флуоресцентной молекулярной томографии времени жизни (fluorescence molecular lifetime tomography, FMLT) все чаще используется в молекулярном имиджинге мелких лабораторных животных [1 - 17]. Это связано, во-первых, с высокой информативностью такого параметра, как время жизни флуоресценции. Этот параметр является наиболее чувствительным к изменениям молекулярного окружения флуоресцентных биосенсоров и позволяет получить полезную информацию о пространственно-временных характеристиках процессов, происходящих в клетках и молекулах тканей мелких животных [18]. Во-вторых, хорошо развитый и широко применяемый в экспериментальной онкологии метод флуоресцентной микроскопии времени жизни [18 -21] оказывается недостаточно эффективным при визуализации глубоко залегающих в теле мыши флуоресцентных включений, и возникает необходимость в
3D томографической реконструкции пространственных распределений времени жизни. Совсем недавно в [22] нами описан оригинальный метод FMLT и представлены первые результаты его экспериментального тестирования на фантоме с флуорофором. Метод основан на асимптотической аппроксимации функции источника флуоресценции во временном домене, реконструкции обобщенной функции распределения параметров флуоресценции (fluorescence parameter distribution function, FPDF) и последующем разделении пространственных распределений коэффициента поглощения флуорофора и времени жизни флуоресценции также во временном домене. При этом мы использовали геометрию обратного рассеяния и мезо-скопический режим регистрации данных, который зачастую оказывается наиболее эффективным при исследованиях мелких лабораторных животных [23 -28]. Известно, что в ситуации, когда расстояние между источником и детектором менее 8- 10 мм, диффузионное приближение не работает [29, 30] и для рас-
четов функций и матриц чувствительности приходится использовать либо подходы, основанные на приближениях высокого порядка к уравнению переноса излучения [31 - 34], либо статистический метод Монте-Карло (Monte Carlo, MC) [2, 7, 22, 35 - 42]. В [22] мы использовали метод MC. Однако в указанной статье основной акцент был сделан на представление эксперимента по сканированию фантома с флуорофо-ром, метода решения обратной задачи FMLT, а также поиск эффективных стратегий формирования массивов измерительных данных и матриц чувствительности. Подход же, который мы использовали для моделирования временных функций рассеяния точки (temporal point spread functions, TPSF) и функций чувствительности, был только анонсирован и не был описан сколько-нибудь подробно. Настоящая статья задумана как своеобразное дополнение статьи [22] и призвана восполнить указанный информационный пробел. Здесь мы подробно излагаем алгоритм, который используется для MC моделирования и реализован в программном коде с рабочим названием TurbidMC.
Начало разработок алгоритмов и программ, реализующих метод Монте-Карло для целей биомедицинской оптики [43, 44], было связано с необходимостью верификации оптического метода и выбора для его описания адекватных аналитических моделей. В основе алгоритмов MC - отслеживание траекторий отдельно взятых фотонов с момента их входа в среду и до момента выхода из среды. При этом расчеты траекторий фотонов осуществляются путем последовательного моделирования элементарных событий: рассеяния, поглощения, отражения (или преломления) на границе, свободного пробега и флуоресценции. Первоначальная цель таких расчетов заключалась в моделировании оптических сигналов детекторов. Сегодня в мире существует достаточное количество специализированных программных комплексов, позволяющих решить эту задачу. Все они, как правило, обладают великолепным интерфейсом пользователя и расширенными возможностями распараллеливания вычислений. Примерами являются Geometry and Tracking (GEANT4) toolkit [45], Online Object Oriented MC (O3MC) tool [46], Molecular Optical Simulation Environment (MOSE) [47], MC Software ValoMC [48] и др. Ситуация несколько изменилась с началом использования метода MC для задач томографии. Методы диффузионной оптической томографии (diffuse optical tomography, DOT) и флуоресцентной молекулярной томографии (fluorescence molecular tomography, FMT) ставят целью решение обратной задачи относительно пространственных распределений оптических параметров или параметров флуоресценции соответственно. В большинстве случаев такая задача линеаризуется и сводится к решению системы линейных алгебраических уравнений. Матрица системы -суть матрица чувствительности, в которую записываются значения функций чувствительности, опреде-
ляющие весовые вклады каждой точки (малой ячейки дискретизации) объекта в сигналы, регистрируемые детекторами от источников. Оценка этих весовых вкладов по результатам MC моделирования как раз и является важнейшей задачей ламинарной DOT и ме-зоскопической FMT. Таким образом, алгоритмы MC, используемые для целей томографии, помимо моделирования траекторий, также должны реализовывать процедуры вычисления функций чувствительности. Прекрасный обзор таких алгоритмов дан, например, в [39]. Авторы указанной статьи выделяют три группы методов MC: прямо-сопряженные методы MC [35 -37, 40], пертурбационные методы MC [2, 7, 38, 41, 42] и так называемые средне-путевые (mid-way) методы MC [49, 50]. Наиболее часто в FMT применяются методы первых двух групп. Прямо-сопряженный метод вычисляет функцию чувствительности как произведение световых полей, рассчитанных в результате соответственно прямого и сопряженного моделирования. В случае прямого моделирования фотоны мигрируют непосредственно из источника, а в случае сопряженного моделирования условный источник помещается на место детектора. Этот метод хорош, когда не требуется высокое разрешение в пространстве и времени. В противном случае чаще применяется пертурбационный метод, обладающий более высокой вычислительной эффективностью [39]. Он вычисляет функции чувствительности сразу для нескольких детекторов (как правило, рассматривается вся совокупность связей источник-детектор для выбранного источника) за счет запоминания траекторий и вычисления длин отрезков их пересечений с ячейками (вокселями) исследуемого объекта.
В настоящей статье представлен алгоритм, который посредством кода TurbidMC реализует время-разрешенный пертурбационный метод MC. Конечно же, наш алгоритм во многом повторяет алгоритмы предшественников, поэтому многие отрезки изложения могут показаться читателю знакомыми. Однако мы используем оригинальный подход к изложению и исходим не только из общих принципов MC моделирования, но в первую очередь из конкретно поставленной расчетной задачи, связанной с обработкой данных сканирования и реконструкцией фантома с флуорофором [22]. Это дает возможность обратить повышенное внимание на те детали, которые вызвали наибольшие трудности в процессе реализации алгоритма. Так, например, в параграфе 1 более подробно описывается моделирование источника, а в приложении А - способ задания распределения интенсивности излучения внутри критического угла. В параграфе 2 приводятся формулы для расчета функций чувствительности, которые хорошо согласуются с конкретной моделью, выбранной в [22] для решения прямой задачи FMT и реконструкции FPDF. Параграф 3 посвящен представлению и анализу результатов моделирования. Мы сравниваем результаты мо-
делирования флуоресцентных TPSF (fluorescence TPSF, FTPSF) с экспериментальными флуоресцентными временными откликами (fluorescence temporal responses, FTR), полученными в результате сканирования фантома с флуорофором [22]. Показано, что результаты расчетов функций чувствительности зависят от места локализации флуорофора по отношению к паре источник-детектор, и с этим, по-видимому, следует считаться при реконструкции распределений коэффициента поглощения флуорофора и времени жизни флуоресценции. Корректность расчета функций чувствительности подтверждается результатами реконструкции фантома [22]. Пример реконструкции пространственных распределений параметров флуоресценции приводится в параграфе 4 настоящей статьи.
Поскольку программа TurbidMC разработана для решения задач время-разрешенной FMT, мы расширили стандартный подход (см., например, [44]) на случай работы во временном домене посредством введения в рассмотрение дополнительной фазовой координаты, времени. Это позволило количественные оценки всех моделируемых величин получать как функции времени задержки детектора, что нашло свое отражение в соответствующих формульных выражениях, приводимых в параграфах 1 и 2. Другой особенностью кода TurbidMC является то, что в нем реализован как стандартный метод аналогового моделирования, так и приближенный вариант метода неаналогового моделирования [51, 52], метод пробных частиц [53]. Этот метод при моделировании акта рассеяния определяет конус возможных направлений пробной частицы после рассеяния и моделирует траекторию такой частицы, но в отличие от метода неаналогового моделирования не рассматривает ветвящиеся траектории. В пункте 1.5 мы уделяем внимание краткому представлению метода пробных частиц и делаем предположение, что в дальнейшем он продемонстрирует свою эффективность в случае очень малых размеров источника и приемника. Тестирование же кода для заданных параметров показало, что метод пробных частиц является более точным, чем метод аналогового моделирования, но требует в значительной степени больше вычислительного времени. По этой причине в настоящей статье приводятся результаты, полученные с помощью только аналогового метода.
Завершая настоящий параграф, хотелось бы добавить, что программа TurbidMC создана на базе многолетнего опыта ФГУП «Российский федеральный ядерный центр - ВНИИ технической физики им. академика Е.И. Забабахина» (РФЯЦ-ВНИИТФ) в области разработок МК программ для расчета переноса различного вида излучений. Примеры: ПРИЗМА [52], IThMC [54], MCTools [55]. Особенно стоит отметить программный инструментарий MCTools, разработанный в среде С++11 специально для быстрого создания специализированных МК программ. Опыт и во
многом инструментарий МСТоок позволили нам существенно снизить трудозатраты на разработку новой программы и сделать ее более выгодной, чем доработка и адаптация для решения поставленных задач уже существующих программ.
1. Моделирование ТРБЖ 1.1. Постановка задачи и условия моделирования
Поскольку целью является моделирование ТРЗБб в геометрии сканирования фантома с флуорофором трехканальным световодом (см. [22] и пункт 3.1), уместно выбрать геометрию моделирования, с одной стороны, удобную для проведения расчетов, а с другой - легко приводимую к геометрии сканирования. Такая геометрия представлена на рис. 1.
z
Рис. 1. Геометрия моделирования
Задача моделирования формулируется следующим образом. Рассчитать временное распределение количества фотонов, испущенных из источника S в оптически мутную среду и попавших в детектор D под апертурным (критическим) углом меньше &ap = 8,20 (рис. 1). Источник S и детектор D расположены на поверхности раздела среда-воздух и представляют собой круги радиусом Rs=Rd = 0,2 мм, направленные по оси OZ с координатами центров: X = - Rsd /2, ys = 0, zs = 0) и (xd=Rsd /2, yd = 0, Zd = 0). В соответствии с конструктивными особенностями трехканального зонда [22] для Rsd рассмотрены значения 1,1; 2,2 и 3,3 мм. На рис. 1 среда занимает верхнее полупространство с границей раздела Z = 0 и определяется оптическими параметрами: коэффициентом поглощения | = 0,01 мм4, коэффициентом рассеяния |s=2,63 мм-1, показателем преломления n = 1,521 и показателем анизотропии g = 0,62. В среду помещен флуорофор, представляющий собой цилиндр с осью, параллельной оси OY, и радиусом Rf = 1,5 мм. Ось цилиндра определяется координатами (xf, y, Zf = 2,4 мм). При этом значение xf от расчета к расчету варьируется в пределах от - 4 мм до 4 мм с шагом 0,5 мм (всего 17 значений). На рис. 1 показан случай, когда xf = 0.
Флуорофор имеет такие же оптические параметры, что и среда, квантовый выход флуоресценции равен у = 0,2; а время жизни флуоресценции - т = 900 пс.
В настоящей статье мы полагаем, что флуоресцентный фотон мигрирует в среде с такими же оптическими параметрами, что и фотон возбуждающего излучения, но при этом он не способен вызвать повторное флуоресцентное излучение. Излучение флуоресцентного фотона происходит с временной задержкой tf, характеризуемой вероятностью
p(tf) = 1 exp I- f |.
(1)
Нижнее полупространство на рис. 1 имитирует воздух с коэффициентом преломления п = 1. Распространение излучения в нижнем полупространстве не моделируется. Источник фотонов описывается равномерным распределением интенсивности излучения на торцевой поверхности (равномерно в круге) и нормальным распределением каждой из двух угловых компонент $х и 8у, определяющих угол 9 с нормалью к поверхности ХОУ (рис. 2). Формулы нормального распределения угловых компонент, которые мы используем, имеют вид
p(Sx)=
p(S y ) =
1 exp i S2
^/Улст v . 2ct2.
1 eXp i ' Sy'
2ct2
(2)
2ln2
Более подробное обоснование выбора вида пространственно-углового распределения интенсивности излучения дано в приложении А.
Рис. 2. К описанию углового распределения источника
Направление фотона мы описываем направляющими косинусами, как это и принято при моделировании с использованием метода МС [43], т.е. вектором П = (Ох,ОУ,О,): О2 + О2 + О2 = 1. Связь между направляющими косинусами и углами:
Qx = cos Sx, Qy = cos S y, Qz = sj 1 - cos2 Sx - cos2 Sy .
(3)
При этом отслеживаются условия О2 + О у < 1 и ^2 +^У < . Детектором Б фиксируются все фо-
тоны, попадающие на поверхность детектора и имеющие угол к ее нормали меньше чем 9ар:
Qz > - cos Sa
(4)
Если условие (4) не выполняется, то фотон зеркально отражается от поверхности детектора. При пересечении границ раздела среда-воздух, среда-источник и среда-детектор фотон может отразиться от поверхности, а может пройти через нее в соответствии с законами френелевского рассеяния [56]. Взаимодействие фотона с границами раздела мы учитываем посредством расчета весов по формулам, которые приведены, например, в [43]. Фотоны, прошедшие через границы раздела, далее не моделируются.
1.2. Моделирование источника
В MC программе TurbidMC фотон при моделировании описывается следующими девятью параметрами: координатами x, y и z; направляющими косинусами Qx, Qy и Qz; временем t; статистическим весом w и типом частицы m: m = 1 в случае фотона возбуждающего излучения и m=2 в случае флуоресцентного фотона. Чтобы описать источник, необходимо определить начальные параметры каждого фотона источника. Это делается на основании знания следующих исходных данных: координат центра круга источника Xs, ys и Zs (Xs = - Rsd/2, ys = Zs = 0), радиуса круга источника Rs (Rs = 0,2 мм), единичного вектора ориентации круга n = (nx,ny,nz):nl + n2y + nZ = 1 (nx = ny = 0, nz = 1), апертурного угла Sap (Sap = 8,20), длительности импульса источника At (At = 6 пс).
Координаты точки старта фотона источника определяются с использованием метода Неймана [57] по следующему алгоритму:
Шаг 1. Задаются псевдослучайные числа и равномерно распределенные на интервале [0,1].
Шаг 2. Вычисляются величины x=Rs(2|x - 1), y=Rs(2|y - 1).
Шаг 3. Проверяется условие x2 + y2 < R2. Если условие не выполняется, возвращаемся на шаг 1. В противном случае переходим к шагу 4.
Шаг 4. Вычисляем координаты точки старта x, y и z по формулам x = x+xs, y = y + ys и z = 0.
Направление полета фотона источника также определяется с использованием метода Неймана по следующему алгоритму:
Шаг 1. Задаются псевдослучайные числа и ^, распределенные по Гауссу на интервале (0,1).
Шаг 2. Вычисляются величины Sx = x ,
S y =а£ y.
Шаг 3. Проверяется условие S;x + S2 < S^ . Если условие не выполняется, возвращаемся на шаг 1. В противном случае переходим к шагу 4.
Шаг 4. Вычисляем величины Qx = cosSx и Qv = cosSv.
S
ар
Шаг 5. Проверяется условие О2 + О5 ^ 1. Если условие не выполняется, возвращаемся на шаг 1. В противном случае переходим к шагу 6.
Шаг 6. Вычисляем направляющий косинус Оъ = ^ 1 - О: - О 2 и определяем направление полета фотона источника как О = (Ох, Оу, Ог).
Время старта фотона источника определяется в соответствии с формулой t = Д^, где ^ - псевдослучайное число, равномерно распределенное на интервале [0,1]. Также каждому фотону источника приписывается статистический вес м> = 1 и тип частицы т = 1.
1.3. Моделирование траектории
На рис. 1 схематически показана траектория, образованная одним фотоном источника. Фотон всегда стартует из круга источника и имеет направление внутри конуса апертурного угла. Узлы траектории на рис. 1 обозначают точки старта фотона, точки его взаимодействия с атомами среды и точки пересечения им границ раздела сред (далее точки взаимодействия). Регистрируются все фотоны, пересекающие круг детектора Б и при этом влетающие под апертур-ным углом к нормали поверхности. Фотоны, которые пересекают круг детектора под большими углами, отражаются зеркально.
По-разному моделируется взаимодействие фотона в зависимости от того, где оно происходит. Если взаимодействие происходит в среде вне цилиндра с флу-орофором, то координаты новой точки взаимодействия и новое значение времени частицы определяются в соответствии с формулами:
г,+1 = г, + • l, ti+1 = ti +-• l,
(5)
где г, и г« - радиус-векторы, определяющие координаты текущей и новой точек взаимодействия; О - единичный вектор, определяющий направление фотона в точке г,; ti и ti+\ - текущее и новое значение времени фотона, с - скорость света в вакууме, I - длина свободного пробега, определяемая по формуле [43]
l = -
log I
la + |
(6)
где % - псевдослучайное число, равномерно распределенное на интервале [0,1].
Вне флуорофора после розыгрыша точки нового взаимодействия происходит случайный выбор типа взаимодействия. С вероятностью p = | /(| + |s) фотон поглощается и далее не моделируется, в противном случае разыгрывается угол рассеяния 6 в соответствии с заданным показателем анизотропии g. Мы полагаем, что косинус угла рассеяния cos 6 распределен согласно формуле Хени-Гринштейна [58]:
1 1 - g2
P(cos 6) = — ---^-—. (7)
4л (1 - 2g • cos 6 + g2)3/2
Тогда косинус угла рассеяния определяется по формулам
cos 0 = 1 -2| и
cos 0 = — <¡1 + g2 -
2g I
1 -i
1 + g (2|-1)_
(8)
в случаях изотропного и анизотропного рассеяния соответственно.
После нахождения косинуса угла рассеяния происходит перерасчет нового направления фотона Он = (Ох,,+1, Оу,,+1, Ог,,н) с учетом разыгранной также величины азимутального угла ф. Если азимутальный угол распределен изотропно, что справедливо в случае неполяризованного излучения, то направляющие косинусы Ох,м, Оу,i+\ и Ог,,н вычисляются по следующим формулам [59]
Q r,+, = Q cos 0-
1 - cos 0
1-Q2 ,i
с (Q y,¡ sin ф - Qxí Qz ,¡ cos ф),
Q „ ,+1 = Q „, cos 0 +
1 - cos 0
(9)
1 - Q2 ,,
x(Q x ¡ sin ф + QyJQz, cos ф),
Qz,,+1 = Qz, cos 0-^/(1 - cos 0)(1 -Q2,) cos ф,
где cos ф = 1 - 2|, sinф = 2^|(1 -|) .
Если взаимодействие происходит внутри флуорофора, то схема моделирования меняется. В процессе моделирования взаимодействия фотона возбуждающего излучения всегда образуются два фотона: рассеянный и флуоресцентный. Для рассеянного фотона меняется направление движения и статистический вес. Направление определяется по изложенной выше схеме. Статистический вес рассеянного фотона умножается на вероятность рассеяния, т.е. w+1 = /(| + |). Тип рассеянного фотона не меняется (m = 1). Это означает, что он по-прежнему способен вызвать флуоресценцию. Этот фотон записывается в стек частиц. Флуоресцентному фотону приписывается вес Wj+1 = w¿Y|a /(| + |s) и время t,+1 = t¿+1 + tf, отличающееся от времени рассеянного фотона t¿+1 на величину tf. Временная задержка tf определяется из экспоненциального распределения (1) методом обратных функций [57]:
tf =-т log I. (10)
Направление образованного флуоресцентного фотона моделируется, исходя из предположения изотропного излучения, а тип частицы меняется на m = 2. Это означает, что этот фотон не способен вызвать флуоресценцию и его моделирование как внутри флуорофора, так и вне его далее происходит по схеме, описанной выше для фотона возбуждающего излучения.
Описанная схема моделирования взаимодействия фотонов подобна схеме с неявным поглощением. Она позволяет на траектории фотона с параметром т = 1 при его прохождении через флуорофор создать достаточное для набора статистики количество моделируемых флуоресцентных квантов с параметром т = 2. Однако эта схема, как и любая схема с неявным поглощением, требует введения предельного статистического веса, при котором необходим обрыв траектории фотона. Флуоресцентный фотон моделируется дальше, без записи в стек частиц.
Моделирование траектории фотона прекращается, когда:
1) фотон влетел в приемник с регистрацией вклада;
2) фотон пересек границу раздела среда-воздух или среда-источник;
3) фотон поглотился, что может произойти с фотоном с параметром т = 1 вне области флуорофора или с флуоресцентным фотоном;
4) статистический вес фотона меньше предельного значения. В случае зеркального отражения от границ раздела моделирование траектории продолжается.
Описанный алгоритм моделирования миграции фотонов в мутной среде реализован в МС программе TurbidMC. Программа написана на языке С++ стандарта С++11, использует стандартные возможности распараллеливания по потокам и программный инструментарий МСТоо1з [55], о котором упоминалось во Введении.
1.4. Выходные параметры при моделировании ТРБЕ
Промежуточными результатами счета по программе ТиГ^МС являются следующие оценки:
N
•Л1 6 тк, = 1)-w(AJ),
1 к,
N
•Л2=Х
^(к 6 О, т, = 1)-pk,)
Л1 =ХХ1(, 6 О, тк] = 2)-w(pk,),
, к,
N
Л2 = £
]
^1(к, 6 Б, т, = 2) - pk])
(11)
• ({Т}) = £ £ 1( 6 Б, т, = 1, ^ 6 {Т}) - ^^pkJ),
•з2({Т}) = £
] к, N
6 Б, тк, = 1, ^ 6 {Т}) - w()
•1({Т}) = ХХ1( 6 Б, т] = 2, ^ 6{Т})- ^^pk]),
] к,
N
•Л2({Т}) = £
]
6 Б, т^ = 2, г, 6 {Т}) - w(^к,)
В приведенных формулах 1(-) - индикаторная функция, принимающая значение 1 в случае наступления события, и значение 0 в противном случае, а w (pk]) -статистический вес фотона pk]■ в момент завершения моделирования траектории. Здесь ] и N - соответственно индекс и число историй, каждая из которых начинается вводом фотона в среду и заключается в моделировании всех траекторий, им образованных, к
- это индекс траектории ]-й истории. Событие pk]■ 6 Б означает, что фотон pk]■ внес свой вклад в приемник Б, т.е. пересек круг приемника в конусе апертурного угла. Значения тк]- = 1 или тк]- = 2 указывают на тип фотона pkJ
- фотон возбуждающего излучения или флуоресцентный фотон соответственно. Событие ^ 6 {Т} означает, что значение времени фотона должно попадать в выбранный временной отрезок Т, при этом расчеты проводятся для всех отрезков, на которые разбит исследуемый временной диапазон {Т }.
На основании оценок (11) рассчитываются: 1) интенсивность всех попавших в детектор фотонов возбуждающего излучения на один фотон источника
• =
л
N''
(12)
2) среднеквадратическое отклонение величины (12)
5, =
•?-( Л,1)2/ N N (N -1) :
(13)
3) интенсивность всех попавших в детектор флуоресцентных фотонов на один фотон источника
•»=#,
(14)
4) среднеквадратическое отклонение величины (14)
52 =
Л2 ( л1)2/ N
N (N -1) :
(15)
5) временное распределение интенсивности попавших в приемник фотонов возбуждающего излучения (ТР8Б) на один фотон источника
•3({Т })) =
•1({Т })
N '
(16)
6) среднеквадратические отклонения величин (16)
(17)
5з({Т}) =
лз2({Т})- [ Л3({Т})]2/ N
N (N -1)
7) временное распределение интенсивности попавших в приемник флуоресцентных фотонов (БТР8Р) на один фотон источника
•1({Т })
• 4({Т })) = •
N
(18)
8) среднеквадратические отклонения величин (18)
§4({^}) =
Л({Т}) -[Л({Г})] / ж
N (N -1)
(19)
Величины (12) - (19) записываются в т-файлы вывода формата МЛТЬЛБ.
1.5. Метод пробных частиц
Характерной особенностью расчета функции чувствительности с помощью пертурбационного метода МС (см. параграф 2) является то, что в ячейках объекта суммируются вклады только от модельных траекторий, которые вносят вклад в сигнал рассматриваемого детектора. Вероятность таких траекторий мала, вследствие чего возникает проблема ограниченности выборки траекторий, дающих вклад в функцию чувствительности. Решить эту проблему позволяет метод неаналогового моделирования [51, 52], который рассматривает ветвящиеся (расщепленные) траектории. Однако этот метод требует усложнения алгоритма моделирования: введения стека для точек ветвления, больших затрат на хранение данных и большего времени вычислений, необходимого для построения ветвящегося процесса. Поэтому в программе ТиГММС в качестве метода, альтернативного методу аналогового моделирования, реализован приближенный метод пробных частиц [53], который не приводит к необходимости рассмотрения ветвящихся траекторий.
Метод пробных частиц использует схему с неявным поглощением и заключается в следующем. Для каждого акта рассеяния определяется конус возможных направлений пробной частицы после рассеяния, как показано на рис. 3. Если рассеяние имеет место за пределами сферы детектора (рис. 3), происходит расщепление на две частицы: аналоговую и пробную.
Направление частицы
Точка рассеяния
Если рассеяние произошло внутри сферы детектора, работает обычная схема аналогового моделирования. Направление аналоговой частицы определяется в соответствии с формулой (7). Если это направление не попадает в конус возможных направлений пробной частицы, то траектория аналоговой частицы моделируется далее, а частице присваивается вес
Щ = щ .ра(Та),
(20)
где щ - вес частицы в точке 2 на рис. 3, ра ( Та) - весовой множитель аналоговой частицы, связанный с неаналоговым рассеянием. Если направление попадает в конус, то траектория аналоговой частицы далее не моделируется, а рассматривается пробная частица. Ее направление после рассеяния y¥t определяется из распределения
р(СОБ у) =
1
1 - СОБ у
(21)
где у - угол между направлением из точки рассеяния 2 на центр сферы детектора и направлением Вес частицы вычисляется по формуле
Щ = щ -р( (V,),
(22)
где pt ( - весовой множитель пробной частицы, связанный с неаналоговым рассеянием. Рассчитывается расстояние от точки 2 до точки пересечения круга детектора (рис. 3). Если пересечение возможно, оценивается вклад в детектор с весом
щ = щ .р((V,)• ехр[-(Ца + |(V,)], (23)
где ( - расстояние между точкой рассеяния 2 и точкой пересечения круга детектора. Далее траектория пробной частицы не прослеживается.
до рассеяния
Направление аналоговой частицы после рассеяния
Конус возможных направлений пробной частицы после рассеяния
Направление пробной частицы после рассеяния
Рис. 3. Схема построения траектории методом пробных частиц
Если траектория частицы пересекает круг детектора сразу после рассеяния, то моделирования аналоговой и пробной частиц не выполняется, а в показания детектора заносится результат с весом
= Щ • ехр (-цА^д ),
(24)
где - расстояние между точкой рассеяния и точкой пересечения круга детектора.
Следует заметить, что, несмотря на отсутствие необходимости рассматривать ветвящиеся траектории, метод пробных частиц является более трудоемким, нежели метод аналогового моделирования. Это связано с необходимостью вычислять вклады практически в каждом акте рассеяния и при этом рассчитывать оптические расстояния и экспоненты. Тестирование программы ТиГММС для заданных размеров источника и детектора показало, что метод пробных частиц выигрывает у метода аналогового моделирования в точности приблизительно в 1,5 раза, но значительно проигрывает ему во времени счета. По этой причине мы сделали акцент на метод аналогового моделирования и в параграф 3 приводим результаты, полученные именно с его помощью. Следует, однако, иметь в виду, что если, например, впоследствии в эксперименте (пункт 3.1) перейти от оптоволокна к бесконтактному режиму регистрации, то размеры источников и детекторов уменьшатся на порядок. В этом случае с большой долей вероятности метод пробных частиц окажется более эффективным, нежели аналоговый метод.
2. Моделирование функций чувствительности 2.1. Модель решения прямой и обратной задач ЕЫЬТ
Прежде чем рассмотреть, как моделируются функции чувствительности в программе ТиГ^МС, покажем, каким образом эти функции определяются применительно для выбранной модели реконструкции томографических изображений.
Время-разрешенный сигнал флуоресценции, возбуждаемой мгновенным источником, испускающим импульс из точки г, в момент времени 4 = 0, который регистрируется на границе среды в точке га в момент времени то есть РТР8Б, может быть записан в виде (см., например, [14, 22, 60])
Г1 (г,, г,, о ~ jфe (г-г, о ® Е(г, 0 ® О1 (г,-Г, ^ 3г ,(25)
V
где фe (г, Г) - плотность фотонов возбуждающего флуоресценцию излучения, Е (г, 0 - функция распределения флуорофора, О (г - г', (- Г) - функция Грина флуоресценции, ® - оператор временной свертки. Свертка фe (г, () ® Е (г, () представляет собой не что иное, как функцию источника флуоресценции
(г, I) = ф" (г, 0 ®Е (г, I) =
= Г'фe (г,,')ехр Г-^ 1йГ'.
т(г) ^ ' Я т(г))
(26)
где ц/ (г) и х (г) - пространственные распределения коэффициента поглощения флуорофора и времени жизни флуоресценции соответственно.
Проблема реализации БМЬТ во временном домене связана со сложной зависимостью функции источника флуоресценции (26) от распределения х (г). В результате не представляется возможным без специаль-
ных ухищрений построить линейную модель реконструкции относительно х (г). Проблему решают одним из следующих способов: 1) переходят в частотную область, где формулы оказываются более простыми [1, 4, 5, 8]; 2) применяют так называемый метод мультиплексирования [2, 3, 6, 7, 9 - 11]; 3) строят нелинейную модель и решают систему нелинейных уравнений [12 - 15]. Все эти методы имеют свои недостатки. В качестве альтернативного подхода в [22, 60] мы предложили аппроксимировать функцию источника (26) асимптотическим приближением, справедливым для фотонов переднего фронта ТР8Б:
S/ (г, ф
уца/ (г) - 4оа2 / п г |2 х(г) + 40^2/ п
фe (г, г),
(27)
где О - коэффициент диффузии в среде на длине волны возбуждающего излучения. Такой подход позволяет построить линейную модель реконструкции относительно функции распределения параметров флуоресценции БРОБ
/ (г) =
4 БоуЦа/ (г)/ п х(г)у2(/) + 40 / п,
(28)
где V (?) - средняя скорость движения центра масс мгновенного распределения фотонов вдоль их средней траектории [61, 62]. Модель описывается линейным интегральным уравнением Фредгольма 1-го рода
Г1 (г,, г,, Г) = Г / (г) Г/ (г,, г,, г, Г) й3 г .
(29)
При этом для функции чувствительности Ш/(г,, га, г, 0, ответственной за реконструкцию /(г), справедливо простое выражение
г
(Г,, Га, Г, Г) = Г фe (Г - Г,, Г')О/ (Га - Г, Г - Г'. (30)
0
Таким образом, решение обратной задачи БМЬТ можно получить в два шага. Первый - реконструкция / (г) для различных скоростей V (0, второй - разделение распределений ц/ (г) и х (г) посредством решения системы алгебраических уравнений. Метод реализован и верифицирован в [60] и [22] соответственно численным и физическим экспериментами.
2.2. Алгоритм моделирования функций чувствительности
Для расчета функций чувствительности в программе ТиГММС реализован так называемый вок-сельный подход, согласно которому исследуемый объект разбивается на множество кубических воксе-лей (ячеек) посредством наложения на него регулярной сетки. В процессе моделирования траекторий фотонов (пункт 1.3) вычисляются длины отрезков пересечения траекторий с ячейками. Именно эти длины играют ключевую роль. В результате удается рассчи-
тать веса (весовые коэффициенты) для каждой ,-й ячейки, положение которой определяется радиус-вектором г,, и таким образом определить дискретную функцию чувствительности
г
Ж (г, г,, г,, г) = { ф® (г, - г, г')в/ (г - г,, г - г')йг' .(31)
0
Из анализа выражения (31) следует, что для того, чтобы найти вес Ж/ (гЛ г^, г,, г), необходимо:
1) для каждого фотона возбуждения у-й истории рассчитать вес (гЛ г,) при его миграции из точки г в ячейку г,;
2) учесть изменение веса при поглощении фотона возбуждения и испускании флуоресцентного фотона;
3) для каждого ку -го флуоресцентного фотона рассчитать вес (гЛ г,) при его миграции из г,- в г^;
4) просуммировать результирующие веса по всем историям и флуоресцентным фотонам в каждой истории;
5) учесть, что вклад в сигнал вносят только флуоресцентные фотоны, удовлетворяющие условию 4,- < г.
Будем иметь
Wj (r, r,) = Wjfi exp
-Xi a (r, )ij (r,)
(32)
где Wjfl - начальный вес фотона возбуждения у-й истории, ца (г) - коэффициент поглощения в ячейке г,- на длине волны возбуждения, (г,) -длина отрезка пересечения ячейки г, траекторией фотона возбуждения у-й истории, ру- - число ячеек, которые пересекают фотон возбуждения у-й истории, мигрируя из точки г^ в ячейку г,. Вероятность поглощения фотона возбуждения и испускания флуоресцентного фотона в ячейке г, дается выражением (см., например, [41])
Ц j (r) = y{l - exp [-|aa/Lj (r, )]};
(33)
где Ьу (г,) - расстояние, которое прошел фотон возбуждения у-й истории в флуорофоре. Вес флуоресцентного фотона определим по формуле
wi (r, rrf) = wk ,o exp
Pi +qt]
- X lf (r )lk] (r,)
i=P] +1
(34)
где wуo - начальный вес Щ-го флуоресцентного фотона, ца (г,-) - коэффициент поглощения в ячейке г,-на длине волны флуоресценции, (г,) - длина отрезка пересечения ячейки г, траекторией ку-го флуоресцентного фотона, qkу - число ячеек, которые пересекает ку-й флуоресцентный фотон, мигрируя из ячейки г,, в точку г, Таким образом, суммируя веса по всем историям и флуоресцентным фотонам в каждой истории и учитывая временную зависимость, для дискретной функции чувствительности получим
W/ (r, rd, r,, t) = XX Wj.oT exp
j=1 k]
-Xlf (r, )lj (r,)
< {l - exp [-|a/Lj (r, )]| x exp <1(tk, < t).
Pi +«]
- X lf fr )lk] (!,■)
(35)
Если принять во внимание, что в нашем случае ца = ц£ = Ца/ = Ца, и разложить сомножитель 1 -ехр [- ца/Ьу (г,)] в ряд Тейлора, то формула (35) может быть упрощена до
W/ (r, rd, r,, t) =
N
= XXW].oY exp
]=1 k]
x|f (r,)L] (r, )1(tk] < t).
P] +«]
- X lf (r, )lk] (r,)
,=1
(36)
Именно формула (36) была запрограммирована и использовалась в программе TurbidMC для расчета функции чувствительности в режиме флуоресценции. Близкие формулы были получены также в [7, 41]. Заметим, что в случае выбора для расчета метода пробных частиц мы могли бы использовать ту же формулу (36). А вот в случае метода неаналогового моделирования потребовалось бы еще просуммировать (36) по всем ветвящимся траекториям.
3. Результаты моделирования и сравнение с экспериментальными данными
Основной целью разработки программы TurbidMC и моделирования FTPSF и функций чувствительности, как это указано во введении, была обработка результатов физического эксперимента по сканированию фантома с флуорофором и реконструкции пространственных распределений времени жизни флуоресценции. Подробное описание эксперимента и метода реконструкции FMLT-изображений не является задачей настоящей статьи. Это сделано в [22]. Однако для корректного сравнения модельных и экспериментальных данных имеет смысл дать краткое описание экспериментальной части проделанной работы.
3.1. Эксперимент по сканированию фантома с флуорофором
Эксперимент проводился в ФГУ «Федеральный исследовательский центр «Фундаментальные основы биотехнологии РАН». Экспериментальная установка была сконструирована на основе системы счета одиночных фотонов с корреляционной обработкой во времени (TCSPC) компании Becker & Hickl GmbH (Берлин, Германия): детектор счета одиночных фотонов PMC-100, TCSPC-плата регистрации SPC-150. Для возбуждения флуоресценции использовали импульсный лазер суперконтинум SC-450-6 компании Fianium (Саутгемптон, Англия) с акустооптическим модулятором длины волны испускаемого света. Для
регистрации экспериментальных FTR в геометрии обратного рассеяния использовался трехканальный оптоволоконный зонд, в котором четыре волокна закреплены линейно с межосевым расстоянием 1,1 мм. Каждое из 4 волокон имеет сердцевину диаметром 400 мкм и числовой апертурой 0,2. Первое волокно использовалось для ввода возбуждающего пучка света, последующие три - для регистрации FTR. В качестве фантома использовался параллелепипед размером 33x40x15 мм3 из однородного имитирующего биологические ткани материала INO Biomimic на основе полиуретана с добавлением рассеивающих частиц из диоксида титана. Фантом имеет цилиндрическое отверстие радиусом Rf = 1,5 мм, в которое помещается флуоресцирующий раствор Cy5 с концентрацией 5-10 -7 моль/литр. Флуорофор возбуждался на длине волны 640 нм, для регистрации использовался полосный фильтр HQ720/60 и блокирующий фильтр HQ650LP компании Chroma Technology Corp. (Вермонт, США). Зонд перемещался по поверхности фантома c шагом 0,5 мм с помощью микрометрической подвижки. Сначала сканирование происходило в направлении, перпендикулярном оси цилиндра с флуороформ. Затем после прохождения 19 положений зонд смещался в начало следующей строки и вновь двигался в том же направлении. Потом осуществлялся переход к следующей строке и так далее по зигзагообразной траектории. Всего было проска-нировано 19 строк. Таким образом, размер области сканирования составил 9x 9 мм2. Для каждого из 361 положения зонда регистрировались три FTR для трех различных расстояний между источником и детектором: 1,1; 2,2 и 3,3 мм. Кроме того, чтобы иметь возможность выполнить необходимую предобработку экспериментальных FTR (параграф 4) и сравнить с ними модельные FTPSF (пункт 3.2), был измерен инструментальный отклик с использованием технологии, изложенной в [5].
3.2. Результаты моделирования FTPSF
и сравнение с экспериментальными FTR
Для решения задачи моделирования TPSF и FTPSF, постановка которой дана в пункте 1.1, всего выполнен 51 расчет, по 17 для каждого значения расстояния Rsd: 3,3; 2,2 и 1.1 мм. При этом координата Xf оси флуорофора (рис. 1) принимала значения в пределах от -4 мм до 4 мм с шагом 0,5 мм. Моделирование выполнялось с использованием аналогового метода, при этом среднее количество моделируемых историй для одного расчета составило около 1,6-108 историй. Погрешность для интенсивности фотонов как возбуждающего излучения J1 ^, так и флуоресценции J]) не превышала 1,5 %. Среднее время одного расчета на многопроцессорном вычислительном комплексе РФЯЦ-ВНИИТФ в варианте параллельного счета на 16 потоках было равно 2 часам.
На рис. 4 и 5 представлены в виде графиков некоторые полезные для анализа результаты. На
рис. 4а представлена зависимость от х/. Из рис. 4а, например, следует, что, как и следовало ожидать при равных значениях оптических параметров объекта и флуорофора, место локализации флуорофора никак не влияет на интенсивность возбуждающего излучения. Рис. 4б, где дана зависимость /2) от х/, показывает, что максимум интенсивности флуоресценции соответствует минимальному смещению флуорофора относительно х/= 0 в сторону детектора. „ хЮ"4
и с я
о 4
5 я
sr
£ 3 н
о
о
я
и „
5 2
я
-T-
■—i— ■— ■—i— I I .—i—:
-
R =1,1 мм sa R =2,2 мм sa R =3,3 мм sa
- -
а)
х 10"
Р 7
§4
я
« „
5 3
я
о»
S2
R = 1,1 MM sa =2,2 мм sa
R R
sd '
/ / 1 N \ !
1 ' / ..... - : \\
1
, мм
б)
Рис. 4. Интенсивности возбуждающего излучения /1) и флуоресценции /2^ как функции координаты х/
Рис. 5 дает представление о значении времени задержки детектора t, на которое следует ориентироваться при синтезе измерительных данных для реконструкции и расчете функций чувствительности.
Чтобы сравнить модельные БТР8Б и зарегистрированные в эксперименте FTR, мы выполнили следующую последовательность действий. Во-первых, аппроксимировали смоделированные по программе TurbidMC ступенчатые БТР8Б (рис. 5) гладкими кривыми. Во-вторых, свернули результаты аппроксимации с инструментальным откликом, который был измерен в эксперименте. И, наконец, в-третьих, согласовали геометрию рис. 1 с геометрией сканирования фантома. Для аппроксимации и свертки мы ис-
пользовали стандартные операторы пакета МЛТЬЛБ inteгp(•) и сопу(^) соответственно. Для согласования геометрий мы выбрали некоторую среднюю строку сканирования (взяли 5-ю строку) и построили графики, похожие на графики рис. 4б, на которых каждому отчету данной строки ставилось в соответствие значение интеграла экспериментального БТЯ по времени. Сопоставление максимумов графиков, построенных по экспериментальным данным, и максимумов графиков рис. 4б позволило нам согласовать отсчеты моделирования х/, когда менял свое место локализации флуорофор, и отсчеты сканирования, при котором зонд перемещался по грани фантома перпендикулярно оси флуорофора [22]. В результате мы получили информацию о том, какой модельный импульс БТРБР какому экспериментальному БТЯ соответствует.
На рис. 6 представлены графики, на которых сравниваются свернутые с инструментальным откликом модельные БТРБР (серые кривые) и экспериментальные БТЯ (черные кривые). На каждом графике также показаны модельные БТРБР до свертки (штриховые кривые) и инструментальный отклик (пунктирные кривые), с которым свертка выполнялась. Для удобства сравнения все импульсы нормированы на
свои максимальные значения. Рис. 6 наглядно демонстрируют в целом неплохое согласование серых и черных кривых, полному совпадению которых друг с другом мешает присутствие шума на экспериментальных БТЯ. Таким образом, результаты, представленные в настоящем параграфе, позволяют заключить, что моделирование БТРБР с помощью программы Тиг-bidMC выполняется вполне корректно.
10""
® 10"
? 10-
Я 10"
10"
возбуждение флуоресценция
10
10'
10'
10'
10
время, пс
Рис. 5. Временные распределения ТР8Г (сплошная кривая) и РТР8Р (пунктирная кривая) для случая Я^ = 3,3 мм и х/= 0
8000
яппп
2000 Рис
4000 6000 время, пс
яппп
8000
яппп
яппп
8000
яппп
4000 6000 время, пс
4000 время, пс
6. Сравнение модельных ГТР8Р и экспериментальных РТЯ: слева направо - х/= -3; 0; 3 мм, сверху вниз - Язс1 = 3,3; 2,2; 1,1 мм
яппп
3.3. Результаты моделирования функций чувствительности
В численном эксперименте [60] расчет функций чувствительности выполнялся для среды с приведенным однородным распределением флуорофора. В результате эти функции выглядели одинаково для одних и тех же расстояний между источником и детектором и имели вид, близкий к симметричному. В принципе такой подход, как показали исследования [60], вполне эффективен, если реконструкции подлежит флуорофор, распределение которого представляет собой локальные включения малого размера. Однако в случае эксперимента [22] допущение о локальном характере распределения флуорофора не может считаться полностью оправданным, поскольку цилиндр с флуорофором занимает хоть и меньшую, но все же значимую часть фантома. Поэтому было естественно предположить, что функции чувствительности для пар источник-детектор вблизи цилиндра с флуорофором и на удалении от него должны как-то отличаться друг от друга. Таким образом, целью наших расчетов в геометрии рис. 1 было, в частности, исследование того, как меняется характер распределения функций чувствительности при изменении взаимоположения пары источник-детектор и места локализации флуорофора. В пункте 3.2 мы показали, что геометрия рис. 1 легко приводится к геометрии сканирования фантома [22]. Таким образом, мы знаем, какой модельный расчет функции чувствительности какой паре источник-детектор одной строки сканирования соответствует. Это означает, что наши пространственно-зависимые расчеты функций чувствительности мы можем использовать для реконструкции фантома по экспериментальным данным его сканирования (параграф 4).
Как и в случае моделирования БТР8Б (пункт 3.2), по программе ^ГММС мы выполнили 51 расчет функций чувствительности, варьируя значения и х/. Для времени задержки детектора изначально было выбрано значение / = 200 пс, которое всегда соответствует передним фронтам БТР8Б и лежит внутри диапазона 75 - 90 % от их максимальных значений.
На рис. 7 некоторые результаты моделирования представлены в виде 2Б-изображений. Эти 2Б-изображения являются сечениями рассчитанных 3Б-распределений функций чувствительности плоскостью ХО2, проходящей через центры источника S и приемника D, и несут наиболее полезную информацию о характере распределений в целом. Рис. 7 демонстрирует, что, когда х/ = 0 (центральный столбец изображений), распределения функции чувствительности выглядят близкими к осесимметричным (симметричным относительно плоскости УО2 в 3Б-случае). Когда же ось флуорофора перемещается, то распределения теряют симметричный вид, и их центр тяжести смещается в сторону места локализации
флуорофора. Эффект смещения распределений в сторону флуорофора можно объяснить тем, что при всех прочих равных условиях больший вклад в сигнал флуоресценции вносят те ячейки объекта, которые расположены ближе к флуорофору. Таким образом, мы видим, что характер пространственного распределения функции чувствительности, ответственной за реконструкцию БРББ, зависит от места локализации флуо-рофора в объекте. В [22] показано, что учет такой зависимости при реконструкции позволяет заметно улучшить качество флуоресцентных томограмм.
4. Пример реконструкции распределений параметров флуоресценции
Для реконструкции пространственных распределений параметров флуоресценции требовалось:
1) выполнить предобработку экспериментальных
2) восстановить БРББ /(г) для различных значений средней скорости миграции фотонов V (0;
3) разделить пространственные распределения коэффициента поглощения флуорофора Ца/ (г) и времени жизни флуоресценции х(г). Предобработка FTR состояла в компенсации шума
с использованием МЛТЬЛБ оператора sgolayfilt(•), деконволюции сглаженных импульсов с инструментальным откликом с помощью оператора deconvlucy(•) и калибровке с целью определения начала отсчетов времени задержки детектора Для калибровки мы использовали модельные БТР8Б.
Для решения обратной задачи относительно БРББ мы обращали систему линейных алгебраических уравнений, построенную в соответствии с линейной моделью (29). С целью формирования векторов измерительных данных и матриц чувствительности была разработана стратегия [22], согласно которой использовалось не одно, а несколько времен задержки, а именно t = 66; 100; 133; 200 пс. Разумеется, расчеты функций чувствительности выполнялись для всех необходимых сочетаний (Я, /). В пункте 1.1 отмечено, что при моделировании расстояние х/ принимало 17 значений. Поскольку реально в эксперименте в строке сканирования было на два отчета больше (19 отсчетов), мы выполнили экстраполяцию, согласно которой этим крайним отчетам были поставлены в соответствие соседние функции чувствительности (экстраполяция по «ближайшему соседу»). Далее для всех 19 строк сканирования мы брали одинаковые сочетания функций чувствительности. Запись в матрицу дискретных функций чувствительности производилась согласно разработанной в [22] стратегии. Число элементов вектора измерительных данных и соответственно строк матрицы чувствительности принимало значения, кратные 361: 1083 и 722. Реконструкции подлежала область размером 20*20*15 мм3 (рис. 8а), размер ячейки был выбран равным 0,1 мм, соответственно сетка имела размер 200*200*150. Та-
ким образом, мы имели сильно недоопределенную систему, для обращения которой требовался алгоритм с регуляризацией. В [22] в качестве такого алгоритма был выбран АЯТ-ТУ-Р^Т. сочетающий ал-
I3
-3 -2-10 1 X, мм
гебраическую реконструкцию (ART) [63], регуляризацию посредством минимизации нормы полной вариации (TV) [64, 65] и мягкую пороговую фильтрацию (FIST) [66].
О 0.039
I с
и
1
о'
II
-2 -1
0 1 2 А", мм
Рис. 7. Информативные 2Б-сечения 3Б-функций чувствительности: слева направо — х/=—3; 0; 3 мм,
сверху вниз — = 3,3; 2,2; 1,1 мм
а) X, мм о 20 V, мм (У) а, мм в) X, мм г)
Рис. 8. Численная модель фантома с флуорофором: (а) реконструируемая область фантома с выделенным фрагментом; (б) увеличенный фрагмент, используемый для визуализации томограмм; (в) и (г) 2- и X- сечения этого фрагмента соответственно. Воспроизведено из [22] в соответствии с лицензией СС-ВУ
Функция FPDF f(r) восстанавливалась для трех различных значений v (t). Поэтому при разделении распределений ц^ (r) и т (r) мы имели теперь уже переопределенную систему, которая решалась с помощью известного алгоритма наименьших квадратов с QR-факторизацией (QR-factorization least squares, LSQR) [67], реализованного в пакете MATLAB оператором lsqr(-).
На рис. 8а вместе с полной областью реконструкции показан также ее фрагмент, который мы используем для визуализации результатов. На рис. 8б этот фрагмент представлен в увеличенном виде. Также показаны горизонтальная ^-сечение) и вертикальная (Х-сечение) плоскости, которые используются для визуализации информативных сечений фрагментов томограмм. На рис. 8в и 8г квадратом и стрелками
соответственно показаны область сканирования и направления ввода света при сканировании. Видно, что область сканирования и, следовательно, зона чувствительности охватывают далеко не всю область реконструкции и даже не весь фрагмент, который мы выбрали для визуализации результатов.
На рис. 9 представлен лучший результат, который был получен при разделении параметров флуоресценции: на рис. 9а показано распределение коэффициента поглощения флуорофора ^ (г), а на рис. 96 -распределение времени жизни флуоресценции т (г). На каждом из рис. 9а и 96 слева дана 3Б-визуализация фрагмента томограммы. А остальные два изображения каждого рисунка показывают 2Б 2-и Х-сечения этого фрагмента. В нижней части справа каждого из рисунков показана палитра, шкала которой проградуирована в обратных миллиметрах в случае рис. 9а и в пикосекундах в случае рис. 96.
Из визуального сравнения рис. 8 и 9 следует, что получены далеко не идеальные изображения. Как и следовало ожидать, область вне зоны чувствительности не воспроизводится совсем или воспроизводится частично с искажениями. На изображении рис. 9а просматриваются артефакты вблизи плоскости сканирования. Но в целом оба распределения выглядят вполне адекватно. Распределение т (г) имеет достаточно четкие границы и, в общем, корректно отображает номинальное значение времени жизни (900 пс). Поскольку знание априорной информации об объекте позволило нам построить модель области реконструкции (рис. 8а), мы использовали такие характеристики, как коэффициент корреляции и показатель отклонения [22] для количественной оценки качества полученных изображений. Результаты расчетов характеристик сведены в табл. 1.
Рис. 9. Результаты разделения распределений параметров флуоресценции: коэффициент поглощения флуорофора Ц#(т); время жизни флуоресценции т(г). Воспроизведено из [22] в соответствии с лицензией СС-БУ
Табл. 1. Количественные характеристики качества изображений рис. 9
Количественная характеристика Номер рисунка Значение
Коэффициент корреляции 9а 0,7946
9б 0,7419
Показатель отклонения 9а 0,5462
9б 0,6751
Данные табл. 1 в целом подтверждают сделанный выше вывод о том, что получены не идеальные, но в целом адекватные томограммы. Тот факт, что коэффициент корреляции все-таки далек от 1, а показатель отклонения - от 0, по-видимому, можно объяснить тем, что на изображениях существуют области вне зоны чувствительности, где значения параметров флуоресценции вообще не воспроизводятся. Таким образом, приведенный пример реконструкции распределений параметров флуоресценции позволяет сделать вывод о том, что расчет функций чувствительности по программе TurbidMC выполняется вполне корректно.
Заключение
В статье описан алгоритм программы ТиГ^МС, реализующей пертурбационный метод Монте-Карло и предназначенной для моделирования временных функций рассеяния точки и функций чувствительности для задач время-разрешенной флуоресцентной молекулярной томографии. Спецификой алгоритма является то, что он ориентирован на работу с конкретным методом томографии времени жизни флуоресценции [22], в основе которого реконструкция обобщенной функции, содержащей распределения коэффициента поглощения флуорофора и времени жизни флуоресценции, которые далее подлежат разделению. Корректность работы программы проверена сравнением тестовых расчетов флуоресцентных временных функций рассеяния точки с данными эксперимента по сканированию фантома с флуорофором трехканальным зондом в мезоскопическом режиме обратного рассеяния. Показано, что в случае, когда флуорофор занимает часть области исследования, размерами которой нельзя пренебречь, существует
зависимость результатов расчетов функций чувствительности от места локализации флуорофора. Продемонстрировано, что, если априорная информация о месте локализации используется, удается получить вполне адекватные реконструкции распределений параметров флуоресценции с использованием метода, представленного в [22]. Заметим, что в противном случае во избежание потери качества томограмм в [22] рекомендован двухшаговый подход, согласно которому изначально восстанавливается распределение коэффициента поглощения флуорофора, а затем полученная апостериорная информация используется для реконструкции распределения времени жизни флуоресценции. Полагаем, что разработанная программа вполне готова для эффективной эксплуатации при решении задач экспериментальной время-разрешенной флуоресцентной молекулярной томографии. Можно добавить, что программа обладает до конца не исследованным потенциалом. Ожидается, что реализованный в программе метод пробных частиц в дальнейшем непременно проявит свою эффективность в случае малых размеров источников и детекторов, характерных, например, для бесконтактного режима облучения и регистрации сигнала, к которому мы планируем перейти в дальнейших экспериментах.
References
[1] Gao F, Zhao H-J, Tanikawa Y, Yamada Y. A linear, featured- data scheme for image reconstruction in timedomain fluorescence molecular tomography. Opt Express 2006; 14(16): 7109-7124. DOI: 10.1364/OE.14.007109.
[2] Kumar ATN, Raymond SB, Boverman G, Boas DA, Bacskai BJ. Time resolved fluorescence tomography of turbid media based on lifetime contrast. Opt Express 2006; 14(25): 12255-12270. DOI: 10.1364/OE.14.012225.
[3] Kumar ATN, Raymond SB, Dunn AK, Bacskai BJ, Boas DA. A time domain fluorescence tomography system for small animal imaging. IEEE Trans Med Imaging 2008; 27(8): 1152-1163. DOI: 10.1109/TMI.2008.918341.
[4] Nothdurft RE, Patwardhan SV, Akers W, Ye Y-P, Achilefu S, Culver JP. In vivo fluorescence lifetime tomography. J Biomed Opt 2009; 14(2): 024004. DOI: 10.1117/1.3086607.
[5] Gao F, Li J, Zhang L, Poulet P, Zhao H, Yamada Y. Simultaneous fluorescence yield and lifetime tomography from time-resolved transmittances of small-animal-sized phantom. Appl Opt 2010; 49(16): 3163-3172. DOI: 10.1364/AO.49.003163.
[6] Raymond SB, Boas DA, Bacskai BJ, Kumar ATN. Lifetime-based tomographic multiplexing. J Biomed Opt 2010; 15(4): 046011. DOI: 10.1117/1.3469797.
[7] Chen J, Venugopal V, Intes X. Monte Carlo based method for fluorescence tomographic imaging with lifetime multiplexing using time gates. Biomed Opt Express 2011; 2(4): 871-886. DOI: 10.1364/BOE.2.000871.
[8] Gao F, Li J, Zhang W, Yi X, Wang X, Zhang L, Zhou Z, Zhao H. A CT-analogous scheme for time-domain diffuse fluorescence tomography. J Xray Sci Technol 2012; 20(1): 91-105. DOI: 10.3233/XST-2012-0321.
[9] Rice WL, Kumar ATN. Preclinical whole body time domain fluorescence lifetime multiplexing of fluorescent proteins. J Biomed Opt 2014; 19(4): 046005. DOI: 10.1117/1.JBO.19.4.046005.
[10] Hou SS, Rice WL, Bacskai BJ, Kumar ATN. Tomographic lifetime imaging using combined early- and late-arriving photons. Opt Lett 2014; 39(5): 1165-1168. DOI: 10.1364/OL.39.001165.
[11] Rice WL, Shcherbakova DM, Verkhusha VV, Kumar ATN. In vivo tomographic imaging of deep-seated cancer using fluorescence lifetime contrast. Cancer Res 2015; 75(7): 12361243. DOI: 10.1158/0008-5472.CAN-14-3001.
[12] Cai C, Zhang L, Zhang J, Bai J, Luo J. Direct reconstruction method for time-domain fluorescence molecular lifetime tomography. Opt Lett 2015; 40(17): 4038-4041. DOI: 10.1364/OL.40.004038.
[13] Zhang L, Cai C, Lv Y, Luo J. Early-photon guided reconstruction method for time-domain fluorescence lifetime tomography. Chin Opt Lett 2016; 14(7): 071702. DOI: 10.3788/COL201614.071702.
[14] Cai C, Zhang L, Cai W, Zhang D, Lv Y, Luo J. Nonlinear greedy sparsity-constrained algorithm for direct reconstruction of fluorescence molecular lifetime tomography. Biomed Opt Express 2016; 7(4): 1210-1226. DOI: 10.1364/BOE.7.001210.
[15] Cai C, Cai W, Cheng J, Yang Y, Luo J. Self-guided reconstruction for time-domain fluorescence molecular lifetime tomography. J Biomed Opt 2016; 21(12): 126012. DOI: 10.1117/1.JBO.21.12.126012.
[16] Zhang P, Liu J, Hui H, An Y, Wang K, Yang X, Tian J. Linear scheme for the direct reconstruction of noncontact time-domain fluorescence molecular lifetime tomography. Appl Opt 2020; 59(26): 7961-7967. DOI: 10.1364/AO.398967.
[17] Cheng J, Zhang P, Cai C, Gao Y, Liu J, Hui H, Tian J, Luo J. Depth-recognizable time-domain fluorescence molecular tomography in reflective geometry. Biomed Opt Express 2021; 12(7): 3806-3818. DOI: 10.1364/BOE.430235.
[18] Becker W. Fluorescence lifetime imaging - techniques and applications. J Microsc 2012; 247(Part2): 119-136. DOI: 10.1111/j.1365-2818.2012.03618.x.
[19] Wang XF, Periasamy A, Herman B, Coleman DM. Fluorescence lifetime imaging microscopy (FLIM): instrumentation and applications. Crit Rev Anal Chem 1992; 23(5): 369-395. DOI: 10.1080/10408349208051651.
[20] Datta R, Heaster TM, Sharick JT, Gillette AA, Skala MC. Fluorescence lifetime imaging microscopy: fundamentals and advances in instrumentation, analysis, and applications. J Biomed Opt 2020; 25(7): 071203. DOI: 10.1117/1.JBO.25.7.071203.
[21] Dmitriev RI, Intes X, Barroso MM. Luminescence lifetime imaging of three-dimensional biological objects. J Cell Sci 2021; 134(9): jcs254763. DOI: 10.1242/jcs254763.
[22] Konovalov AB, Vlasov VV, Samarin SI, Soloviev ID, Sa-vitsky AP, Tuchin VV. Reconstruction of fluorophore absorption and fluorescence lifetime using early photon mesoscopic fluorescence molecular tomography: a phantom study. J Biomed Opt 2022; 27(12): 126001. DOI: 10.1117/1.JBO.27.12.126001.
[23] Abou-Elkacem L, Bjorn S, Doleschel D, Ntziachristos V, Schulz R, Hoffman RM, Kiessling F, Lederle W. High accuracy of mesoscopic epi-fluorescence tomography for non-invasive quantitative volume determination of fluorescent protein-expressing tumours in mice. Eur Radiol 2012; 22(9): 1955-1962. DOI: 10.1007/s00330-012-2462-x.
[24] Ozturk MS, Lee VK, Zhao L, Dai G, Intes X. Mesoscopic fluorescence molecular tomography of reporter genes in bioprinted thick tissue. J Biomed Opt 2013; 18(10): 100501. DOI: 10.1117/1.JBO.18.10.100501.
[25] Yang F, Ozturk MS, Zhao L, Cong W, Wang G, Intes X. High-resolution mesoscopic fluorescence molecular tomography based on compressive sensing. IEEE Trans Biomed Eng 2015; 62(1): 248-255. DOI: 10.1109/TBME.2014.2347284.
[26] Tang Q, Tsytsarev V, Frank A, Wu Y, Chen C-W, Er-zurumlu RS, Chen Y. In vivo mesoscopic voltage-sensitive dye imaging of brain activation. Sci Rep 2016; 6: 25269. DOI: 10.1038/srep25269.
[27] Azimipour M, Sheikhzadeh M, Baumgartner R, Cullen PK, Helmstetter FJ, Chang W-J, Pashaie R. Fluorescence laminar optical tomography for brain imaging: system implementation and performance evaluation. J Biomed Opt 2017; 22(1): 016003. DOI: 10.1117/1.JB0.22.1.016003.
[28] Ozturk MS, Montero MG, Wang L, Chaible LM, Jechlinger M, Prevedel R. Intravital mesoscopic fluorescence molecular tomography allows non-invasive in vivo monitoring and quantification of breast cancer growth dynamics. Commun Biol 2021; 4: 556. DOI: 10.1038/s42003-021-02063-8.
[29] Arridge SR, Schotland JC. Optical tomography: forward and inverse problems. Inverse Probl 2009; 25(12): 123010. DOI: 10.1088/0266-5611/25/12/123010.
[30] Kuz'min VL, Val'kov AYu, Zubkov LA. Photon diffusion in random media and anisotropy of scattering in the Henyey-Greenstein and Rayleigh-Gans models. J Exp Theor Phys 2019; 128(3): 396-406. DOI: 10.1134/S1063776119020109.
[31] Lu Y, Zhu B, Shen H, Rasmussen JC, Wang G, Sevick-Muraca EM. A parallel adaptive finite element simplified spherical harmonics approximation solver for frequency domain fluorescence molecular imaging. Phys Med Biol 2010; 55(16): 4625-4645. DOI: 10.1088/0031-9155/55/16/002.
[32] Kim HK, Lee JH, Hielscher AH. PDE-constrained fluorescence tomography with the frequency-domain equation of radiative transfer. IEEE J Sel Top Quantum Electron 2010; 16(4): 793-803. DOI: 10.1109/JSTQE.2009.2038112.
[33] Guo H, Hou Y, He X, Yu J, Cheng J, Pu X. Adaptive hp finite element method for fluorescence molecular tomography with simplified spherical harmonics approximation. J Innov Opt Health Sci 2014; 7(2): 1350057. DOI: 10.1142/S1793545813500570.
[34] He X, Guo H, Yu J, Zhang X, Hou Y. Effective and robust approach for fluorescence molecular tomography based on CoSaMP and SP3 model. J Innov Opt Health Sci 2016; 9(6): 1650024. DOI: 10.1142/S1793545816500243.
[35] Crilly RJ, Cheong W-F, Wilson B, Spears JR. Forward-adjoint fluorescence model: Monte Carlo integration and experimental validation. Appl Opt 1997; 36(25): 65136519. DOI: 10.1364/AO.36.006513.
[36] Finlay JC, Foster TH. Recovery of hemoglobin oxygen saturation and intrinsic fluorescence with a forward-adjoint model. Appl Opt 2005; 44(10): 1917-1933. DOI: 10.1364/AO.44.001917.
[37] Haykawa CK, Spanier J, Venugopalan V. Coupled forward-adjoint Monte Carlo simulations of radiative transport for the study of optical probe design in heterogeneous tissues. SIAM J Appl Math 2007; 68(1): 253-270. DOI: 10.1137/060653111.
[38] Chen J, Intes X. Time gated perturbation Monte Carlo for whole body functional imaging in small animals. Opt Express 2009; 17(22): 19566-19579. DOI: 10.1364/OE.17.019566.
[39] Chen J, Intes X. Comparison of Monte Carlo methods for fluorescence molecular tomography - computational efficiency. Med Phys 2011; 38(10): 5788-5798. DOI: 10.1118/1.3641827.
[40] Gardner AR, Haykawa CK, Venugopalan V. Coupled forward-adjoint Monte Carlo simulation of spatial-angular light fields to determine optical sensitivity in turbid media. J Biomed Opt 2014; 19(6): 065003. DOI: 10.1117/1.JBO.19.6.065003.
[41] Jiang X, Deng Y, Luo Z, Wang K, Lian L, Yang X, Meglinski I, Luo Q. Evaluation of path-history-based fluorescence Monte Carlo method for photon migration in heterogeneous media. Opt Express 2014; 22(26): 3194831965. DOI: 10.1364/OE.22.031948.
[42] Yao R, Intes X, Fang Q. Direct approach to compute Jaco-bians for diffuse optical tomography using perturbation Monte Carlo-based photon "replay." Biomed Opt Express 2018; 9(10): 4588-4603. DOI: 10.1364/BOE.9.004588.
[43] Wang L, Jacques SL, Zheng L. MCML - Monte Carlo modeling of light transport in multi-layered tissues. Comput Methods Programs Biomed 1995; 47(2): 131-146. DOI: 10.1016/0169-2607(95)01640-F.
[44] Welch AJ, Gardner C, Richards-Kortum R, Chan E, Criswell G, Pfefer J, Warren S. Propagation of fluorescent light. Lasers Surg Med 1997; 21(2): 166-178. DOI: 10.1002/(SICI)1096-9101(1997)21:2<166::AID-LSM8>3.0.CO;2-O
[45] Agostinelli S, et al. Geant4 - a simulation toolkit. Nucl Instrum Methods Phys Res A 2003; 506(3): 250-303. DOI: 10.1016/S0168-9002(03)01368-8.
[46] Doronin A, Meglinski I. Online object oriented Monte Carlo computational tool for the needs of biomedical optics. Biomed Opt Express 2011; 2(9): 2461-2469. DOI: 10.1364/BOE.2.002461.
[47] Ren S, Chen X, Wang H, Qu X, Wang G, Liang J, Tian J. Molecular Optical Simulation Environment (MOSE): A platform for the simulation of light propagation in turbid media. PLoS ONE 2013; 8(4): e61304. DOI: 10.1371/journal.pone.0061304.
[48] Leino AA, Pulkkinen A, Tarvainen T. ValoMC: a Monte Carlo software and MATLAB toolbox for simulating light transport in biological tissue. OSA Continuum 2019; 2(3): 957-972. DOI: 10.1364/OSAC.2.000957.
[49] Serov I, John T, Hoogenboom J. A new effective Monte Carlo midway coupling method in MCNP applied to a well logging problem. Appl Radiat Isot 1998; 49(12): 17371744. DOI: 10.1016/S0969-8043(98)00055-4.
[50] Serov I, John T, Hoogenboom J. A midway forward-adjoint coupling method for neutron and photon Monte Carlo transport. Nucl Sci Eng 1999; 133(9): 55-72. DOI: 10.13182/NSE99-A2072.
[51] Briesmeister J. MCNP - a general Monte Carlo N-particle transport code. Los Alamos National Laboratory Report LA-13709-M 2000.
[52] Kandiev YaZ, Malyshkin GN, Zatsepin OV. Monte Carlo code PRIZMA for calculation of particle transport problems. Proc Joint Int Conf on Supercomputing in Nuclear Applications and Monte Carlo [CD-ROM]; 2010.
[53] Lux I, Koblinger L. Monte-Carlo transport methods: Neutron and photon calculations. Boca Raton: CRC Press; 2000. ISBN: 0-8493-6074-9.
[54] Dorosev AS, Kostjuchenko VI, Samarin SI. Direct account of experimental data uncertainties in modeling of a 160 MeV proton beam interaction with a multilayer Faraday cup [In Russian]. Meditsinskaya Fizika 2014; 2(62): 24-31.
[55] Samarin SI. Certificate of Governmental Registration of Computer Program in FIPS. No. 2018666251; 2018.
[56] Born M, Wolf E. Principles of optics. 7th ed. Cambridge: Cambridge University Press; 1999. ISBN: 0-521-64222-1.
[57] Sobol IM. Numerical Monte Carlo methods [In Russian]. Moscow: "Nauka" Publisher; 1973.
[58] Henyey IG, Greenstein JI. Diffuse radiation in the galaxy. Astrophys J 1941; 93: 70-83. DOI: 10.1086/144246.
[59] Akkerman AF. Modeling of charge particle trajectories in matter [In Russian]. Moscow: "Energoatomizdat" Publisher; 1991. ISBN: 5-283-02924-7.
[60] Konovalov AB, Vlasov VV, Uglov AS. Early-photon reflectance fluorescence molecular tomography for small animal imaging: Mathematical model and numerical experiment. Int J Numer Method Biomed Eng 2021; 37(1): e3408. DOI: 10.1002/cnm.3408.
[61] Lyubimov VV, Kalintsev AG, Konovalov AB, Lyamtsev OV, Kravtsenyuk OV, Murzin AG, Golubkina OV, Mordvinov GB, Soms LN, Yavorskaya LM. Application of the photon average trajectories method to real-time reconstruction of tissue inhomogeneities in diffuse optical tomography of strongly scattering media. Phys Med Biol 2002; 47(12): 2109-2128. DOI: 10.1088/00319155/47/12/308.
[62] Konovalov AB, Vlasov VV, Lyubimov VV. Statistical characteristics of photon distributions in a semi-infinite
turbid medium: Analytical expressions and their application to optical tomography. Optik 2013; 124(23): 60006008. DOI: 10.1016/j.ijleo.2013.04.111.
[63] Gordon R, Bender R, Herman GT. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography. J Theor Biol 1970; 29(3): 471-482. DOI: 10.1016/0022-5193(70)90109-8.
[64] Yu H, Wang G. Compressed sensing based interior tomography. Phys Med Biol 2009; 54(9): 2791-2805. DOI: 10.1088/0031-9155/54/9/014.
[65] Vlasov VV, Konovalov AB, Kolchugin SV. Hybrid algorithm for few-views computed tomography of strongly absorbing media: algebraic reconstruction, TV-regularization, and adaptive segmentation. J Electron Imaging 2018; 27(4): 043006. DOI: 10.1117/1.JEI.27.4.043006.
[66] Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J Imaging Sci 2009; 2(1): 183-202. DOI: 10.1137/080716542.
[67] Paige CC, Sanders MA. LSQR: An algorithm for sparse linear equations and sparse least squares. ACM Trans Math Softw 1982; 8(1): 43-71. DOI: 10.1145/355984.355989.
Приложение А. О способе задания нормального распределения интенсивности источника
Мы исходим из того, что пространственно-угловое распределение источника является независимым по каждому из направлений: р (х, Эх,у, Эу) = Рх (х, Эх) • Ру (у, (рис. 2) и в принципе может определяться с учетом корреляции между парами параметров (х, Эх) и (у, Эу). Тогда в случае нормального распределения и по пространству, и по углу справедливы соотношения
P (х, Эх ) =
1
nJDX
exp
AxЭ2Х ~ 2ВххЭх + Cxx2 Dx
P (У, Эу ) =
1
n D
exp
Ay Эу - 2ВууЭу + Cyy2 ' Dy
(А1)
где Ax, у = 2а2
Cx, у = 2<
Вх, у = kx, у
C
у х, у
Dx,у = Ax,у ' Cx,у Вх,у ,
а а^у, аЭ э и кХуу - дисперсии и
коэффициенты корреляции соответственно. То есть, вообще говоря, для задания нормального распределения источника необходимо определить шесть параметров: Ах, Вх, Сх, Ау, Ву и Су. При отсутствии корреляции и с учетом того, что по х и по у мы берем равномерное распределение, их остается только два: Сх и Су или аэх и аэу. А с учетом симметрии - вообще один: аэх = аэу = а.
Таким образом, задача задания распределения интенсивности источника сводится к определению связи среднеквадратического отклонения а и апертурного угла Эар. Мы предположили, что апертурный угол соответствует значению распределения, равному половине от его максимума (рис. А1).
Тогда будем иметь
А
Рис. А1. К определению среднеквадратического отклонения а
Pmax = P(0) =
1
2na2
^m21 = f (Эар ) = 1
1
Эар
, э^р 1и2 =expI-
exp I---
2na2 l 2а2
^ а =
-\/2ln2
(А2)
То есть приходим к формулам (2) пункта 1.1.
Сведения об авторах
Самарин Сергей Иванович, окончил Челябинский государственный университет, кандидат физико-математических наук, руководитель лаборатории ФГУП «Российский федеральный ядерный центр - ВНИИ технической физики им. академика Е.И. Забабахина». Область научных интересов: Монте-Карло моделирование, ядерная и радиационная безопасность, радиационное повреждение кристаллов, программирование. E-mail: [email protected] .
Коновалов Александр Борисович, окончил Санкт-Петербургский государственный университет аэрокосмического приборостроения, кандидат физико-математических наук, ведущий научный сотрудник ФГУП «Российский федеральный ядерный центр - ВНИИ технической физики им. академика Е.И. Забабахина». Область научных интересов: математическое моделирование, обработка изображений, обратные и некорректные задачи, биомедицинская оптика. E-mail: [email protected] .
Власов Виталий Викторович, окончил Южноуральский государственный университет, руководитель лаборатории ФГУП «Российский федеральный ядерный центр - ВНИИ технической физики им. академика Е.И. Забабахина». Область научных интересов: обработка изображений, томография, распознавание образов, искусственные нейронные сети. E-mail: [email protected] .
Соловьев Илья Дмитриевич, окончил Московский государственный университет им. М.В. Ломоносова, кандидат химических наук, младший научный сотрудник ФГУ "Федеральный исследовательский центр «Фундаментальные основы биотехнологии РАН» Область научных интересов: флуоресцентная микроскопия и спектроскопия, молекулярный имиджинг. E-mail: [email protected] .
Савицкий Александр Павлович, доктор химических наук, профессор, главный научный сотрудник, заведующий лабораторией ФГУ «Федеральный исследовательский центр "Фундаментальные основы биотехнологии РАН». Область научных интересов: молекулярный имиджинг, биология клетки, флуоресцентная микроскопия и спектроскопия, молекулярные нанобиосенсоры, фотобиохимия и фотодинамическая терапия, технология мече-ния на основе лантанидов и фосфоресценции. E-mail: [email protected] .
Тучин Валерий Викторович, доктор физико-математических наук, профессор, член-корреспондент РАН, заведующий кафедрой оптики и биофотоники Саратовского государственного университета, работает также с несколькими другими институтами и университетами, заслуженный деятель науки Российской Федерации, почетный член международных научных обществ SPIE и OPTICA, лауреат ряда престижных премий этих обществ. Область научных интересов: биомедицинская оптика, лазерная медицина, оптическое просветление тканей, нанобиофотоника. E-mail: [email protected] .
ГРНТИ: 29.31.23
Поступила в редакцию 19 февраля 2023 г. Окончательный вариант - 3 апреля 2023 г.
Monte Carlo modeling of temporal point spread functions and sensitivity functions for mesoscopic time-resolved fluorescence molecular tomography
S.I. Samarin1, A.B. Konovalov1, V.V. Vlasov1, I.D. Solovyev2, A.P. Savitsky2, V.V. Tuchin 2,3 1FSUE "Russian Federal Nuclear Center - Zababakhin All-Russia Research Institute of Technical Physics", 456770, Russia, Chelyabinsk Region, Snezhinsk, Vasiliev Str. 13;
2 Bach Institute of Biochemistry, Research Center of Biotechnology of the Russian Academy of Science, 119071, Russia, Moscow, Leninsky Avenue 33, bld. 2;
3 Chernyshevsky Saratov State University, 410012, Russia, Saratov, Astrakhanskaya Str., 83
Abstract
The paper describes a TurbidMC code that implements a perturbative Monte Carlo method to model temporal point spread functions and sensitivity functions for time-resolved fluorescence molecular tomography (FMT). The code is aimed at working with a particular FMT method published earlier (Ref. [22]) which defines the specificity of sensitivity function calculation. The method solves the inverse problem first for a generalized fluorescence parameter distribution function and then calculates separate distributions for the fluorophore absorption coefficient and the fluorescence lifetime. The proper operation of the code was verified through a comparison between fluorescence temporal point spread functions from test calculations and data from experiments where a phantom with a fluorophore was scanned with a three-channel probe in the mesoscopic reflectance regime. An example is given on the reconstruction of fluorescence parameter distributions. It shows that the sensitivity functions are calculated correctly.
Keywords: TurbidMC code, Monte Carlo method, fluorescence molecular tomography, temporal point spread function, sensitivity function, fluorophore absorption coefficient, fluorescence lifetime.
Citation: Samarin SI, Konovalov AB, Vlasov VV, Solovyev ID, Savitsky AP, Tuchin VV. Monte Carlo modeling of temporal point spread functions and sensitivity functions for mesoscopic time-resolved fluorescence molecular tomography. Computer Optics 2023; 47(5): 673-690. DOI: 10.18287/2412-6179-CO-1295.
Authors' information
Sergei Ivanovich Samarin graduated from Chelyabinsk State University, candidate of Physics and Mathematics, laboratory head at FSUE "Russian Federal Nuclear Center - Zababakhin All-Russia Research Institute of Technical Physics". Areas of scientific interest: Monte Carlo simulations, nuclear and radiation safety, radiation damage of crystals, and programming. E-mail: [email protected] .
Alexander Borisovich Konovalov graduated from Saint-Petersburg State University of Aerospace Instrumentation, candidate of Physics and Mathematics, leading scientist at FSUE "Russian Federal Nuclear Center - Zababakhin All-Russia Research Institute of Technical Physics". Areas of scientific interest: mathematical modeling, image processing, inverse and ill-posed problems, biomedical optics. E-mail: [email protected] .
Vitaly Viktorovich Vlasov graduated from South-Ural State University, laboratory head at FSUE "Russian Federal Nuclear Center - Zababakhin All-Russia Research Institute of Technical Physics". Areas of scientific interest: image processing, tomography, pattern recognition, artificial neural networks. E-mail: [email protected] .
Ilya Dmitrievich Solovyev graduated from Lomonosov Moscow State University, candidate of Chemistry, junior researcher at Bach Institute of Biochemistry, Research Center of Biotechnology of the Russian Academy of Science. Areas of scientific interest: fluorescence microscopy and spectroscopy, molecular imaging. E-mail: [email protected] .
Alexander Pavlovich Savitsky is a doctor of Chemistry, professor, chief researcher, laboratory head at Bach Institute of Biochemistry, Research Center of Biotechnology of the Russian Academy of Science. Areas of scientific interest: molecular imaging, cell biology, fluorescence microscopy and spectroscopy, molecular nanobiosensors, photobio-chemistry and photodynamic therapy, and lanthanide and phosphorescent labeling technology. E-mail: [email protected] .
Valery Victorovich Tuchin is a doctor of Physics and Mathematics, professor, corresponding member of the RAS, Optics and Biophotonics department at Chernyshevsky Saratov State University, working with several other institutes and universities. He was awarded the Honored Science Worker of the Russia, Fellow and Laureate of SPIE h OPTICA International Societies. Areas of scientific interest: biomedical optics, laser medicine, tissue optical clearing, and nanobiophotonics. E-mail: [email protected] .
Received February 19, 2023. The final version - April 3, 2023.
