Система компьютерной реконструкции 3D-распределений нейронов

По данным доклада Всемирной организации здравоохранения (ВОЗ), неврологическими нарушениями страдает около одного миллиарда человек во всем мире. При этом большой интерес вызывают нейродегенеративные заболевания, связанные с гибелью нейронов. К ним относятся болезни Паркинсона и Альцгеймера, рассеянный склероз и др., которые десятилетиями развиваются в скрытой форме и нередко проявляются на пике профессиональной и социальной активности человека. Несмотря на то что лечить больных начинают сразу после появления симптомов, заболевания быстро прогрессируют, приводя к инвалидизации и летальному исходу.

Поэтому одной из важнейших задач нейронаук является разработка экспериментальных моделей социально-значимых заболеваний, в частности, связанных с гибелью нейронов, в том числе в определенной части головного мозга — черной субстанции (ЧС), что приводит к развитию болезни Паркинсона (БП). Такие модели предназначены для разработки новых методов и технологий диагностики и лечения подобных заболеваний. Характерной особенностью БП является то, что она развивается в течение длительного времени (нескольких лет) в скрытой форме — без проявления симптомов, т.е. в «доклинической» стадии. Только после потери «порогового» числа нейронов черной субстанции (»70%) начинают проявляться характерные симптомы — нарушение моторики, и заболевание переходит в клиническую стадию.

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

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

1) приготовление замороженных срезов;

2) иммуноцитохимическое выявление дофаминергических нейронов и их аксонов по содержащемуся в них ключевому ферменту синтеза дофамина — тирозин-гидроксилазе, и приготовление микроскопических препаратов;

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

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

В качестве источника экспериментальных данных для построения модели БП используются цифровые микроскопические изображения нейронов (МИН) и волокон срезов головного мозга экспериментальных животных (разрешение 0.0117 мкм2/пиксел2). Выжившие нейроны на фотографиях представляют собой темные округлые клетки со светлым ядром и выявляются на фронтальных серийных срезах компактной части черной субстанции (КЧС) толщиной 20 мкм (рис. 1). Они являются ключевым звеном регуляции моторного поведения. Прогрессирующая дегенерация нейронов приводит к развитию БП.

р1

Автоматизация анализа МИН и отростков нейронов позволяет снизить материальные затраты на порядок, а временные затраты на два порядка. Однако помимо задач распознавания нейронов на каждом изображении срезов ЧС встает задача упрощения анализа характера распределения нейронов по всему объему черной субстанции.

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

В процессе проектирования системы реконструкции ЗD-распределений нейронов на основе модели БП и серий изображений 2D-cpe30B мозга было выяснено, что ее функционирование должно осуществляться в виде цикла (конвейера) последовательных шагов. Такими шагами должны являться следующие:

1) выравнивание слоев — выделение на изображениях срезов согласованных между собой областей анализа и введение локальных (2D) координатных систем срезов, удобным образом связанных с выделенными областями, определение глобальной трехмерной (3D) системы координат и соотнесение с ней локальных координат всех срезов;

2) выделение нейронов — редактирование и согласование по срезам данных автоматического распознавания нейронов, определение их локальных координат, пересчет локальных 2D-координат к ЗD-глобальным с целью последующей 3D-визуализации самих нейронов;

3) кластеризация нейронов — построение двумерных непрерывных распределений плотности нейронов в срезах на основе нейробиологической модели БП и аппроксимации распределений моделями гауссовых смесей, известной как задача кластеризации;

4) выравнивание кластеров — определение типов кластеров и выравнивание их в соседних слоях по типу и затем по положению с целью получения плавных и соответствующих модели БП трехмерных объектов — 3D-поверхностей постоянного уровня распределений плотности нейронов;

5) 3D-реконструкция распределений плотности нейронов по объему ЧС на основе аппроксимации поверхностей постоянного уровня сплайнами для межслойной интерполяции данных.

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

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

Модельно распределения нейронов принято представлять в виде некоторых скоплений — кластеров, причем основными интересующими специалистов характеристиками кластеров являются количество содержащихся в них нейронов и пространственное расположение самих кластеров.

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

р2

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

1. Общая эллиптическая кривая задается достаточно простым математическим выражением (а и b — полуоси, φ — угол наклона):

621

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

2. В полярных координатах эллиптическая кривая задается формулой:

621.1

что позволяет почти равномерно распределить точки по кривой (для ψ=2п/k) и отобразить их на цифровом экране.

3. Квадратичная форма в экспоненте Гауссова распределения:

622

при фиксированном ее значении задает эллипс, что позволяет рассматривать последний как линию постоянного значения (уровня) распределения.

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

К счастью, на сегодняшний день известны хорошие алгоритмы решения приведенной выше формализованной задачи. Одним из таких широко известных алгоритмов, позволяющих эффективно работать с большими объемами данных, является ЕМ-алгоритм. Его название происходит от слов «expectation-maximization», что переводится как «ожидание-максимизация». Это связано с тем, что каждая итерация содержит два шага — вычисление математических ожиданий (expectation) и максимизацию (maximisation). Перейдем к изложению адаптированной к данному проекту версии ЕМ-алгоритма.

р3

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

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

Среди преимуществ ЕМ-алгоритма можно выделить следующие:

— Надежная статистическая основа.

— Линейное увеличение сложности при росте объема данных.

— Устойчивость к шумам и пропускам в данных.

— Возможность построения желаемого числа кластеров.

— Быстрая сходимость при удачной инициализации.

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

Как отмечалось выше, ЕМ-алгоритм предполагает, что кластеризуемые данные подчиняются линейной комбинации (смеси) нормальных (гауссовых) распределений, плотность вероятности которых имеет вид:

623

где [x], [y] — математическое ожидания, σx, σy — дисперсия, ρ — коэффициент корреляции.

Введем в рассмотрение функцию δ2(х,у) = (ξ2 — 2ρξη + η2)/(l — ρ2), которую будем рассматривать как обобщенное расстояние от (х, у) до ((x),(у)).

Алгоритм предполагает, что данные подчиняются смеси двумерных нормальных распределений. Модель, представляющая собой смесь гауссовых распределений, задается в виде:

624

где F(x,y I i) — нормальное распределение для i-гo кластера, ω — доля (вес) i-гo кластера в гауссовой смеси.

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

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

624.1

Целью алгоритма является оценка средних значений [x], [у], параметров σхi, σyi и ρi весов смеси ω. для функции распределения вероятности F(x, у), описанной выше, таких что они доставляют максимум функции правдоподобием L (принцип максимума правдоподобия — МП) при заданном наборе координат нейронов

624.2

Помимо введенных выше параметров введем еще так называемые «скрытые» параметры g = P(х.,у) — апостериорные вероятности принадлежности (x1,yi) (хi,уj) i-му кластеру. Рекуррентная максимизация функции правдоподобия L на основе ЕМ-алгоритма осуществляется по следующей принципиальной схеме:

шаг Е: по формуле Байеса уточняются (оцениваются) «скрытые» параметры gij:

624

шаг М: с использованием gij уточняются параметры ([xj],[yi]), σxi, σyi, ρi, ωi (на самом деле приведенные ниже оценки максимизируют L при заданных gij):

625

Алгоритм начинает работу с инициализации, т.е. некоторого начального числа кластеров и некоторого приближенного решения, которое может быть выбрано случайно или задано пользователем исходя из некоторых априорных сведений об исходных данных (например, на основе обобщенной модели Паркинсона). Наиболее общим способом инициализации является присвоение элементам ([xi],[yi]) математических ожиданий некоторых априорных значений, а начальные ω, σxi, σyi и ρi задаются как в алгоритме k-means: σxi = σyi = 1, ρi = 0.

Для проверки работы адаптированного алгоритма EM нами было разработано небольшое приложение на MatLab, СЕМ, внешний вид интерфейсной формы которого приведен на рис. 4.

р4

Результаты тестирования СЕМ (MatLab) для проверки адаптированного алгоритма EM на реальных срезах ЧС показали замечательные его характеристики: при заданном количестве кластеров и при достаточно грубых априорных расположениях их центров ([xi],[yi]) распределение нейронов по кластерам устанавливается за 4-5 итераций и параметры кластеров σxi, σyi и ρi хорошо соответствуют ожидаемым значениям.

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

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

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

C(u) = Σ Bj,p (u)Pj,

где С(u) — сплайн, параметризованный переменной и (обычно берут u ∈[0;1]), {Bj,p(u)} — базис В-сплайнов. Набор п контрольных точек {Pj} однозначно задает

форму и положение в пространстве сплайна С(u), который, приближая последовательность {Pj}, может и не проходить ни через одну из этих точек. Другими словами, С (и) — это кривая наперед заданной степени гладкости, максимально плотноприближающая пространственную ломаную с вершинами в точках {Рj}. Характерным свойством сплайнов является их локальность: изменение одной из контрольных точек {Pj} ввиду финитности каждого из В-сплайнов Bj,p(u) приведет к деформации С(u) лишь в ближайшей окрестности Pi.

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

Простая методика моделирования пространственных поверхностей сложной формы включает следующие шаги:

1) поверхность представляется набором принадлежащих ей характерных точек. Выбирается расположение координатных осей u и v. Характерные точки принимаются за опорные точки В-сплайновых примитивов;

2) образуется «окно», включающее 16 опорных точек (по 4 вдоль каждой координатной оси). Происходит вычисление координат текущих точек примитива;

3) по окончании развертывания очередного примитива окно перемещается на один ряд опорных точек по оси u или v. В результате перемещения в окно входят 12 прежних опорных точек и 4 новые точки. Для этого набора опорных точек по выражениям вычисляются координаты текущих точек;

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

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

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

Лофтинг-метод восстанавливает поверхность 3D-объекта неравномерными рациональными В-сплайнами, которые, в частности, хорошо аппроксимируют регулярные контуры объекта. Лофт-объекты строятся путем формирования оболочки по опорным сечениям, расставляемым вдоль некоторой заданной траектории. Оболочка как бы натягивается на сечения вдоль указанного пути, а в результате получается трехмерная модель. Данный метод моделирования прекрасно подходит для тех моделей, форма которых может быть охарактеризована некоторым набором поперечных сечений, как в нашем случае. В основе любого подобного объекта всегда лежат траектория (путь) и одно или более сечений. Путь задает основную линию лофт-объекта и может иметь форму прямой, окружности, произвольной кривой и т.п., а сечения определяют его форму и тоже могут быть самыми разнообразными (в разумных пределах). При использовании нескольких сечений они размещаются вдоль пути по указанному пользователем принципу, а в случае одного сечения данная форма размещается на обоих концах пути. Оба типа структурных элементов — путь и сечения — в общем случае представлены обычными сплайнами, т.е. наборами опорных вершин. Число опорных вершин поперечных сечений может быть произвольным, но это число должно быть во всех задействованных в лофт-объекте сечениях одинаково. Кроме того, если сечение представлено составными формами, то данные формы должны иметь одинаковый порядок вложения. Это означает, что если первое сечение содержит внутри одного контура два других контура, то и все последующие сечения должны быть сформированы по тому же принципу. А если в копии данной составной формы переместить два внутренних сплайна вне исходного, то ее уже нельзя будет указать в качестве второго сечения. Правда, при желании в ряде случаев данное ограничение можно обойти, превратив обычные замкнутые сплайны в разомкнутые. Форма лофт-объекта определяется не только за счет пути и определенного набора сечений — не менее важны положение внутренних сечений вдоль пути и согласование первых вершин каждой формы поперечного сечения. От размещения сечений вдоль пути зависит то, как и в какой момент будет смоделирован переход от одного сечения к другому, а согласование вершин позволяет избежать перекручивания моделей при переходе от сечения к сечению или, наоборот, при необходимости искусственно скручивать объекты. Качество лофтинг-метода сильно зависит от качества контуров, особенно в случае ручной обводки.

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

— Плотные каркасы показывают большее число деталей, чем разреженные.

— Плотные каркасы деформируются более однородно, чем разреженные.

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

— Разреженные каркасы используют меньше памяти и отображаются быстрее.

— С разреженными каркасами, как правило, проще и быстрее работать, нежели с плотными.

Лучшее решение достигается путем создания как можно более разреженных каркасов, удовлетворяющих требованиям проекта по визуализации распределений. Впрочем, в программное обеспечение, предназначенное для визуализации распределений нейронов, при разработке закладывается возможность ручного выбора «гладкости» 3D-объектов — плотности каркаса поверхности.

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

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

— на форме представлен следующий набор вкладок: выравнивание слоев, выделение нейронов, кластеризация, выравнивание кластеров, 3D-представление (рис. 5). На каждой вкладке обеспечиваются средства ручного и автоматического редактирования для соответствующего шага процесса.

р5

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

Модуль «выравнивание слоев» предназначен для выравнивания изображений соседних срезов (текущего относительно предыдущего) и выделения эллиптичес

кой области анализа (для текущего среза), с которой связывается система локальных координат данного среза (рис. 6а). Модуль «выделение нейронов» (рис. 6б) предназначен для редактирования координат нейронов, помеченных вручную или автоматически распознанных с помощью алгоритмов распознавания, разработанных соисполнителями по проекту ВЦ им. А.А. Дородницына РАН.

р6

Вкладка «кластеризация» предназначена для автоматического выделения групп нейронов на срезах с использованием адаптированного ЕМ-алгоритма (рис. 7). С помощью модуля «выравнивание кластеров» согласовываются кластеры на соседних срезах по типам. Однако перед согласованием необходимо каждому из кластеров приписать его уникальный на срезе тип. Эта процедура может быть осуществлена как вручную, так и автоматически по критерию ближнего по евклидову расстоянию кластера предыдущего слоя (если кластеры предыдущего слоя уже типизированы).

Последний шаг процесса непосредственно « 3D -реконструкция» (рис. 8), на котором на основе плоских, параллельных контуров, как на формах лофт-объектов, строятся 3D-поверхности пространственных распределений. На этой вкладке пользователь осуществляет визуальный анализ рассчитанных распределений, представленных 3D-объектами. Средства DirectX-3D управления камерой предоставляют возможность быстро и просто реализовать такие функции анализа и визуализации 3D-модели, как изменение точки обзора, повороты, приближение к объекту или удаление от него и т.д.

р7

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

По результатам работы над проектом к настоящему моменту:

— разработана методология реконструкции трехмерного распределения нейронов по данным их распознавания в изображениях серийных срезов головного мозга;

— адаптирован к задаче выделения контуров распределения нейронов в срезах алгоритм кластеризации, основанный на известном ЕМ-методе (оценивание-максимизация);

— разработана полиномиальная межслойная аппроксимация параметров двумерных контуров кластеров с целью уменьшения количества срезов для трехмерной реконструкции плотности нейронов;

— разработано экспериментальное программное обеспечение ВТРН;

— в ПО включена возможность использования выходных данных программных комплексов автоматического распознавания нейронов по изображениям срезов мозга.

Подпишитесь на свежую email рассылку сайта!

Читайте также