Улучшенный метод Монте-Карло для моделирования распространения зондирующего излучения в задачах оптической диффузионной спектроскопии

Математическое моделирование. Оптимальное управление Вестник Нижегородского университета им. Н.И. Лобачевского, 2014, № 1 (1), с. 239-247

УДК 004.94

УЛУЧШЕННЫЙ МЕТОД МОНТЕ-КАРЛО ДЛЯ МОДЕЛИРОВАНИЯ РАСПРОСТРАНЕНИЯ ЗОНДИРУЮЩЕГО ИЗЛУЧЕНИЯ В ЗАДАЧАХ ОПТИЧЕСКОЙ ДИФФУЗИОННОЙ СПЕКТРОСКОПИИ

© 2014 г. А.В. Горшков, 1 М.Ю. Кириллин, 2 В.П. Гергель1

'Нижегородский госуниверситет им. Н.И. Лобачевского 2Институт прикладной физики РАН, Н. Новгород

anton.v. gorshkov@gmail. com

П-сжупила в редакцию 27.11.2013

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

Ключевые слова: моделирование распространения света, моделирование рассеяния, метод Монте-Карло, метод смещенной выборки.

Введение

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

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

Применение метода оптической диффузионной спектроскопии позволяет решать такие задачи, как диагностика и лечение раковых опухолей, в частности, рака груди [1], мониторинг активности зон коры головного мозга [2], планирование фотодинамической терапии [3], мониторинг состояния пациента при хирургическом вмешательстве [4], определение состояния кожных покровов [5] и др.

Общее описание задачи

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

Идея метода Монте-Карло в данной задаче состоит в случайной трассировке набора фотонов в биоткани. Фотоны объединяются в пакеты, каждый пакет обладает весом. Далее понятия «фотон» и «пакет фотонов» будут отождествляться. Пакет фотонов начинает движение от источника излучения. Далее на каждом шаге трассировки случайным образом определяется его направление, величина смещения и погло-

щенный вес. Моделирование пакета завершается либо при его поглощении средой (когда вес пакета становится меньше минимального), либо если он вылетает за границы исследуемого объекта (рис. 1).

Впервые применять метод Монте-Карло для моделирования распространения света в биологических тканях предложил Wilson [6]. Впоследствии этот подход развивался многими учеными для получения более точных результатов. Так, Prahl ввел в модель учет анизотропии и внутреннего отражения света [7]. Wang предложил алгоритм моделирования распространения света в многослойных гетерогенных средах с плоскопараллельной геометрией, а также разработал программную систему MCML, реализующую предложенный алгоритм [8]. Позже метод был успешно усовершенствован для моделирования функциональной диагностики головного мозга человека [9, 10] и решения ряда других прикладных задач.

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

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

В работе [11] для повышения эффективности стандартного метода Монте-Карло предлагается использовать метод существенной выборки [12]. А именно, предлагается выбирать новое направление движения фотона близким к на-

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

Основным недостатком описанного в работе [11] метода является его применимость только к однородному полубесконечному слою. Целью данной работы является обобщение метода на случай гетерогенных биологических тканей с произвольной геометрией.

Описание предлагаемого метода

Предлагаемый улучшенный метод Монте-Карло для моделирования распространения зондирующего излучения в сложных биологических тканях основан на методах, описанных в работах [8, 11].

Общая трассировка фотона в ткани осуществляется на основании правил, описанных в работе [8]. Основными отличиями предлагаемого метода являются:

• поддержка моделирования тканей с произвольными границами [13];

• адаптация метода к задачам оптической диффузионной спектроскопии [14].

Идея улучшения метода состоит в применении техники существенной выборки при выборе нового направления движения пакета фотонов [11].

Вводится коэффициент выживаемости

с 1+™ V

S (Vf , Va ) =

0.1 + 0.9-

fa

2

где V/ - новое направление движения фотона, va -вектор притяжения, є - коэффициент притяжения. Векторы V/ и va имеют единичную норму.

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

Коэффициент выживаемости показывает степень близости выбранного направления движения к вектору притяжения.

Алгоритм выбора нового направления движения фотона состоит в следующем. Генерируется направление движения фотона V/ стандартным способом. После этого генерируется зна-

чение равномерно распределенной на интервале [0, 1) случайной величины £. Если S(vf, уа) > £, тогда сгенерированное направление Vf принимается. Иначе генерируется новое направление, и процесс повторяется.

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

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

Новый вес вычисляется по формуле:

= -

12-/

S(V/ , ^а V2-Здесь м> - вес пакета фотонов, V, - направление движения фотона на предыдущем шаге моделирования (до выбора нового направления), /х) - стандартная функция распределения, используемая для выбора нового направления движения фотона [8].

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

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

На каждом шаге моделирования вектор притяжения вычисляется как

г - г

Г - г

— с.

3 1

1 Г — Г2 1 1 Г — Г1 1

Значения коэффициентов с и с2 зависят от используемой модели распространения излучения. В случае если рассматривается только рассеяние фотонов, эти коэффициенты равны 1 [11]. А если перенос учитывает еще и поглощение среды (данная модификация модели впервые предложена авторами), то

с1 = ехР

С2 = еХР

Здесь г - вектор, начало которого совпадает с источником излучения, а конец - с текущим

положением фотона, Р = ,

ца - коэффициент поглощения, ц, - коэффициент рассеяния и g - коэффициент анизотропии,

г = ,

Г2 = Гі — (4- + 2АК

где га - вектор, начало которого совпадает с источником излучения, а конец - с центром це-

левого детектора, 4- = 3D - транспортная длина, D = 1[ц, (1 - g) + ца ]- коэффициент диффузии,

А = &1Г - длина экстраполяции, е2 = (0, 0, 1).

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

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

Подбор параметров предлагаемого метода

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

Для выбора оптимального набора параметров были проведены соответствующие эксперименты. Решалась задача моделирования распространения зондирующего излучения в голове человека с геометрией слоев, полученных из данных МРТ (2 186 446 треугольников). Использовался источник излучения с длиной волны зондирования 830 нм. Состав слоев и их оптические характеристики описаны ниже (таблица 1).

Для проведения экспериментов использовался кластер из 10 узлов с двумя центральными процессорами 1Пе1 Хеоп L5630 (2.13 ГГц, 4 ядра) и 24 ГБ оперативной памяти на каждом.

Коэффициент притяжения

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

Коэффициент притяжения может принимать значения, меньшие либо равные единице. Для проведения экспериментов были рассмотрены значения данного параметра от 1 до 0.125. Результаты экспериментов приведены на рис. 2.

а = С2

3

Таблица 1

Оптические характеристики слоев головы человека _______при длине волны зондирования 830 нм _____________

Иа, мм-1 И*, мм-1 g п

Кожа головы 0.025 14.3 0.86 1.4

Жировой слой 0.1 10 0.9 1.4

Череп 0.02 25 0.94 1.55

Цереброспинальная жидкость 0.004 1 0.99 1.33

Серое вещество 0.02 25 0.96 1.4

Белое вещество 0.02 26.7 0.85 1.4

Рис. 2. Сравнение мощности излучения на детекторе для улучшенного метода Монте-Карло при различных коэффициентах притяжения

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

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

это, в свою очередь, приводит к существенному увеличению ошибки.

Относительно корректные результаты в данной задаче можно получить при значении исследуемого параметра с 0.25 и меньше. Однако рекомендуется брать коэффициент притяжения, равный 0.125, так как время работы в этом случае всего на 5% больше, а средняя по всем детекторам относительная ошибка (относительно стандартного метода Монте-Карло) в 2.5 раза меньше и составляет примерно 6.8%. При этом относительная ошибка для целевого детектора составляет 1.9% в обоих случаях.

Брать еще меньший коэффициент притяжения в данной задаче смысла не имеет, так как отличие результатов от решения, полученного стандартным методом Монте-Карло, не значительно. Дело в том, что при мониторинге сигнала на детекторе в случае изменения оптических характеристик целевого объема существенным будет считаться изменение сигнала в 2 и более раз. И ошибка в несколько процентов не окажет

Рис. 3. Сравнение мощности излучения на детекторе для улучшенного метода Монте-Карло с учетом и без учета поглощения при вычислении вектора притяжения

значительного влияния на общие результаты мониторинга.

Учет поглощения фотонов средой при вычислении вектора притяжения

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

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

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

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

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

Выбор оптических характеристик при вычислении вектора притяжения

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

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

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

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

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

как она не приводит к дополнительным вычислениям, а значит, более эффективна.

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

Выбор коэффициента при расчете длины экстраполяции

Коэффициент а, который участвует в формуле (7) при вычислении длины экстраполяции, есть константа, определенная по-разному в различных моделях приближений уравнения теории переноса излучения. В данной работе анализируются три различные константы: 2/3 -константа, применяемая в классической диффузионной теории, 0.71 - константа, полученная при численном решении проблемы Милна, и

0.74 - эмпирически подобранная константа, используемая в работе [11].

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

Проведенные эксперименты не выявили значительного различия результатов распределения мощности излучения на детекторах при ис-

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

Проверка корректности предлагаемого метода

Для проверки корректности предложенного в работе метода было выполнено сравнение результатов его работы с результатами стандартного метода Монте-Карло. Алгоритм стандартного метода Монте-Карло для моделирования распространения зондирующего излучения в гетерогенных биологических тканях с произвольными границами слоев описан в работе [14].

Решалась та же задача, что и при выборе оптимальных параметров метода.

Результаты моделирования приведены на рис. 6.

В качестве целевого был выбран детектор, расположенный на расстоянии 30 мм от источника излучения. В результате трассировки 107 фотонов улучшенным методом Монте-Карло удалось получить 8597 фотонов на целевом детекторе. Примерно такое же количество фотонов попадает на целевой детектор только при моделировании 109 фотонов стандартным методом (9302 фотона). При этом результаты моделирования сигнала на детекторах достаточно хорошо совпадают.

При моделировании 107 фотонов стандартным методом Монте-Карло хорошо заметен недостаток статистики, что проявляется в

Рис. 5. Сравнение мощности излучения на детекторе для улучшенного метода Монте-Карло при использовании различных констант в формуле расчета длины экстраполяции

Рис. 6. Сравнение мощности излучения на детекторе для стандартного и улучшенного метода Монте-Карло

значительном расхождении величин мощности излучения на последних 10 детекторах. Отметим, что число фотонов, попавших на целевой детектор, в этом случае равно 94.

Траектории движения фотонов от источника к целевому детектору (в двумерном сечении) представлены на рис. 7. Результаты моделирования 109 фотонов стандартным методом и 107 фотонов улучшенным методом похожи и позволяют визуально выделить тот объем, через который проходит больше всего фотонов. В то же время траектории, полученные в результате моделирования 107 фотонов стандартным методом Монте-Карло, практически никакой полезной информации не несут.

Отметим, что цветом здесь отмечена частота посещения фотоном данной подобласти. Чем

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

Таким образом, трассировки 107 фотонов улучшенным методом Монте-Карло достаточно для получения результатов, которые соответствуют результатам трассировки 109 фотонов стандартным методом. При этом время моделирования в первом случае составляет 1 час 12 минут, а во втором случае - 54 часа 35 минут.

Следует отметить, что при равном количестве трассируемых фотонов для моделирования улучшенным методом требуется несколько больше времени, чем стандартным. Это связано, во-первых, с тем, что в улучшенном методе присутствуют дополнительные расходы на

Рис. 7. Двумерные сечения траекторий движения фотонов от источника излучения к детектору для стандартного метода с 107

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

Замеры времени производились на кластере из 10 узлов, каждый узел содержал 2 центральных процессора 1Пе1 Хеоп L5630 (2.13 ГГц, 4 ядра) и 24 ГБ оперативной памяти.

Заключение

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

Работа выполнена при поддержке гранта Правительства РФ (контракт 11.G34.31.0066) и грантов РФФИ (проекты № 14-02-31549 и № 14-02-00741) при

организационной поддержке Лаборатории информационных технологий ННГУ.

Список литературы

1. Shah N., Cernssi A.E., Jakubowski D. et al. The role of diffuse optical spectroscopy in the clinical management of breast cancer // Dis. Markers. 2004. 19. P. 95-105.

2. Sorvoja H.S.S., Myllyla T.S., Kirillin M.Yu. et al. Non-invasive, MRI-compatible fibreoptic device for functional near-IR reflectometry of human brain // Quantum Electorinics. 2010. 40. P. 1067-1073.

3. Dougherty T.J., Gomer C.J., Henderson B.W. et al. Photodynamic therapy // J. Nat. Cancer Inst. 1998. 90. P. 889-905.

4. Taillefer M.C., Denault A.Y. Cerebral near-infrared spectroscopy in adult heart surgery: systematic review of its clinical efficacy // Can. J. Anesthesia. 2005. 52. P. 79-87.

5. Marghoob A.A., Swindle L.D., Moricz C.Z.M. et al. Instruments and new technologies for the in vivo diagnosis of melanoma // J. Am. Acad. Dermatol. 2003. 49. P. 777-797.

6. Wilson B.C., Adam G. A Monte Carlo model for the absorption and flux distributions of light in tissue // Med. Phys. 1983. 10. P. 824-830.

7. Prahl S.A., Keijzer M., Jacques S.L., Welch A.J. A Monte Carlo model of light propagation in tissue // Proc. SPIE. 1989. Is. 5. P. 102-111.

8. Wang L.V., Jacques S.L., Zheng L.Q. MCML -Monte Carlo modeling of light transport in multi-layered tissues // Comput. Meth. Prog. Biol. 1995. 47. P. 131-146.

9. Boas D.A., Culver J.P., Stott J.J., Dunn A.K. Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head // Opt. Express. 2002. 10. P. 159-170.

10. Chuang C.C., Chen C.M., Hsieh Y.S. et al. Brain structure and spatial sensitivity profile assessing by near-infrared spectroscopy modeling based on 3D MRI data // J. Biophotonics. 2013. б. P. 2б7-274.

11. Chen N. Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry // Appl. Opt. 2007. 4б. P. 1597-1б03.

12. Соболь И.М. Численные методы Монте-Карло. М: Наука, 1973. 312 с.

13. Горшков А.В., Коршунова А.Л. Оптимальный алгоритм поиска пересечений в задаче Монте-Карло моделирования распространения зондирующего излучения в головном мозге человека // Вестник Нижегородского университета им. Н.И. Лобачевского. 2012. № 5(2). С. 73-80.

14. Gorshkov A.V., Kirillin M.Yu. Monte Carlo simulation of brain sensing by optical diffuse spectroscopy // J. Computational Science. 2012. 3. P. 498-503.

BIASED MONTE CARLO METHOD FOR SIMULATION OF PROPAGATION OF SOUNDING RADIATION IN DIFFUSE OPTICAL SPECTROSCOPY

A. V. Gorshkov, M.Yu. Kirillin, V.P. Gergel

The paper describes an effective technique for Monte Carlo simulation of sounding radiation propagation in complex biological tissues as applied to the problems of diffuse optical spectroscopy. The standard Monte Carlo method is modified on the basis of the biased sampling technique. A description of the proposed method is given and its parameters are analyzed. Estimates of the method validity and efficiency are presented.

Keywords: simulation of light propagation, scattering simulation, Monte Carlo method, biased sampling method.

References

1. Shah N., Cerussi A.E., Jakubowski D. et al. The role of diffuse optical spectroscopy in the clinical management of breast cancer // Dis. Markers. 2004. 19. P. 95-105.

2. Sorvoja H.S.S., Myllyla T.S., Kirillin M.Yu. et al. Non-invasive, MRI-compatible fibreoptic device for functional near-IR reflectometry of human brain // Quantum Electorinics. 2010. 40. P. 1067-1073.

3. Dougherty T.J., Gomer C.J., Henderson B.W. et al. Photodynamic therapy // J. Nat. Cancer Inst. 1998. 90. P. 889-905.

4. Taillefer M.C., Denault A.Y. Cerebral near-infrared spectroscopy in adult heart surgery: systematic review of its clinical efficacy // Can. J. Anesthesia. 2005. 52. P. 79-87.

5. Marghoob A.A., Swindle L.D., Moricz C.Z.M. et al. Instruments and new technologies for the in vivo diagnosis of melanoma // J. Am. Acad. Dermatol. 2003. 49. P. 777-797.

6. Wilson B.C., Adam G. A Monte Carlo model for the absorption and flux distributions of light in tissue // Med. Phys. 1983. 10. P. 824-830.

7. Prahl S.A., Keijzer M., Jacques S.L., Welch A.J. A Monte Carlo model of light propagation in tissue // Proc. SPIE. 1989. Is. 5. P. 102-111.

8. Wang L.V., Jacques S.L., Zheng L.Q. MCML -Monte Carlo modeling of light transport in multi-layered tissues // Comput. Meth. Prog. Biol. 1995. 47. P. 131-146.

9. Boas D.A., Culver J.P., Stott J.J., Dunn A.K. Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head // Opt. Express. 2002. 10. P. 159-170.

10. Chuang C.C., Chen C.M., Hsieh Y.S. et al. Brain structure and spatial sensitivity profile assessing by near-infrared spectroscopy modeling based on 3D MRI data // J. Biophotonics. 2013. 6. P. 267-274.

11. Chen N. Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry // Appl. Opt. 2007. 46. P. 1597-1603.

12. Sobol' I.M. Chislennye metody Monte-Karlo. M: Nauka, 1973. 312 s.

13. Gorshkov A.V., Korshunova A.L. Optimal'nyj algoritm poiska peresechenij v zadache Monte-Karlo modelirovanija rasprostranenija zondirujushhego izluchenija v golovnom mozge cheloveka // Vestnik Nizhe-gorodskogo universiteta im. N.I. Lobachevskogo. 2012. № 5(2). S. 73-80.

14. Gorshkov A.V., Kirillin M.Yu. Monte Carlo simulation of brain sensing by optical diffuse spectroscopy // J. Computational Science. 2012. 3. P. 498-503.