В первой статье этой серии мы рассмотрели общую организацию наборов данных МРТ и фМРТ. Во второй статье мы двинулись дальше и исследовали, какие части мозга были активны во время сканирования фМРТ, выполнив корреляционный анализ между данными и идеализированным профилем ответа. Метод сработал довольно хорошо. Мы наблюдали активность в слуховой коре, как и ожидалось, во время слуховой стимуляции, но карта выглядела немного зашумленной, и мы хотели посмотреть, может ли общая линейная модель (GLM) дать нам лучшие результаты.
Что мы сделали до сих пор
Но прежде чем мы продолжим, давайте вспомним, как был собран рассматриваемый нами набор данных. Данные поступают от субъекта, лежащего в сканере МРТ, в то время как слуховой стимул («двухсложные слова») предъявляется периодически. Во время слуховой стимуляции мы ожидаем сигнал, зависящий от b люфта- o xygen- l уровня d (ЖИРНЫЙ), что является косвенным показателем активности мозга, чтобы увеличить.
Ниже вы можете увидеть матрицу дизайна, которую мы создали в предыдущей статье, которая отражает наше предположение о реакции мозга на слуховой стимул. Он состоит из двух компонентов, или регрессоров, как мы их будем называть. Первый регрессор - это константа, которая всегда имеет значение 1, в то время как ожидаемая реакция изменяется от 0 до 1 каждый раз, когда появляется слуховой стимул.
Для нашего корреляционного анализа мы смотрели только на «ожидаемый ответ» и игнорировали постоянную часть. Однако теперь, когда мы хотим использовать GLM, чтобы увидеть, какие части мозга были активными, нам также нужна постоянная часть. Но прежде чем мы действительно подойдем к модели, давайте посмотрим, что такое GLM и как мы можем использовать ее с нашими данными.
Общая линейная модель
Обычно GLM представляет собой анализ множественной регрессии, который пытается объяснить нашу зависимую переменную, сигнал BOLD, посредством линейной комбинации независимых опорных функций или регрессоров, как мы их здесь будем называть. В нашем случае мы уже определили два регрессора в нашей матрице дизайна, постоянную часть и ожидаемую реакцию. GLM определяется следующим образом:
Где i - i -ое наблюдение (объем во времени), а j - j -й регрессор ( постоянный или ожидаемый ответ). В то время как ε представляет собой ошибку. В конечном итоге это означает следующее: ЖИРНЫЙ сигнал y любого заданного воксела - это сумма нашей постоянной части, умноженной на β0, и ожидаемого отклика, умноженного на β1, плюс некоторый ошибка ε. Теперь наша задача - найти веса β0 и β1. Для этого полезно еще раз подумать о том, с какими данными мы имеем дело. Каждый срез имеет вокселы 64x64, и каждый временной курс вокселей содержит 96 томов. Это означает, что один воксель - это вектор с 96 измерениями. А наша матрица дизайна - это матрица размером 96x2. Таким образом, мы можем переписать приведенное выше уравнение в матричном формате:
Где X обозначает нашу матрицу плана 96x2, а b - двумерный вектор весов. Это означает, что наша ошибка e - это разница между Xb и y:
И эта ошибка e - это то, что нам нужно минимизировать, чтобы подогнать GLM к временной шкале вокселов. Один из эффективных способов сделать это - минимизировать сумму квадратов ошибок, как показано ниже:
Мы можем довольно легко найти веса β с помощью некоторых простых умножений матриц. Если вас интересует более подробная информация о том, как добраться до приведенного ниже уравнения, вы можете взглянуть на это и эту ссылку.
Хорошо, теперь давайте поместим это в код, который мы можем запустить с нашими данными. Приведенная ниже функция принимает в качестве входных данных два аргумента, матрицу дизайна и данные фМРТ, и вычисляет параметры GLM.
GLM как функция Python
>>> def do_GLM(X, y):
>>> def do_GLM(X, y):
# Make sure the design matrix has the right orientation
>>> if X.shape[1] > X.shape[0]:
>>> X = X.transpose()
# Calculate the dot product of the transposed design matrix
# and the design matrix and invert the resulting matrix.
>>> tmp = np.linalg.inv(X.transpose().dot(X))
# Now calculate the dot product of the above result and the
# transposed design matrix
>>> tmp = tmp.dot(X.transpose())
# Pre-allocate variables
>>> beta = np.zeros((y.shape[0], X.shape[1]))
>>> e = np.zeros(y.shape)
>>> model = np.zeros(y.shape)
>>> r = np.zeros(y.shape[0])
# Find the beta values, the error and the correlation coefficients
# between the model and the data for each voxel in the data.
>>> for i in range(y.shape[0]):
>>> beta[i] = tmp.dot(y[i,:].transpose())
>>> model[i] = X.dot(beta[i])
>>> e[i] = (y[i,:] - model[i])
>>> r[i] = np.sqrt(model[i].var()/y[i,:].var())
>>> return beta, model, e, r
Как видно из приведенного выше кода, мы вычисляем еще один параметр, о котором мы еще не говорили. Этот параметр r представляет собой коэффициент корреляции нашей подобранной модели и данных, что означает, что квадрат r представляет собой процент отклонения, объясняемого моделью. Мы можем использовать этот параметр, чтобы установить порог для нашей карты, как показано в следующем коде, который запускает нашу функцию «do_GLM» для данных fMRI и отображает результаты.
# Run the GLM
>>> beta, model, e, r = do_GLM(design_matrix, data)
# Reshape the correlation vector r to create a map
>>> r = r.reshape(x_size,y_size)
>>> map = r.copy()
>>> map[map<0.3] = np.nan
# Plot the result
>>> fig, ax = plt.subplots(1,3,figsize=(18, 6))
>>> ax[0].imshow(mean_data, cmap='gray')
>>> ax[0].set_title('1st EPI image', fontsize=25)
>>> ax[0].set_xticks([])
>>> ax[0].set_yticks([])
>>> ax[1].imshow(r, cmap='afmhot')
>>> ax[1].set_title('un-thresholded map', fontsize=25)
>>> ax[1].set_xticks([])
>>> ax[1].set_yticks([])
>>> ax[2].imshow(mean_data, cmap='gray')
>>> ax[2].imshow(map, cmap='afmhot')
>>> ax[2].set_title('thresholded map (overlay)', fontsize=25)
>>> ax[2].set_xticks([])
>>> ax[2].set_yticks([])
>>> plt.show()
Добавление третьего регрессора
На карте выше показаны два кластера, похожие на те, которые мы видели с картой корреляции, которую мы рассчитали в предыдущей статье, но она не сильно улучшилась. Однако преимущество GLM заключается в том, что теперь мы можем изменить матрицу проектирования, добавив больше регрессоров. Чтобы продемонстрировать это, давайте используем средний отклик вокселей из пороговой карты выше в качестве третьего регрессора. Как отмечалось ранее, амплитуда BOLD-сигнала может быть большой в абсолютных числах, но также может варьироваться для отдельных вокселей. Поэтому было бы хорошо стандартизировать все временные интервалы перед их усреднением. Это делается путем вычисления z-показателя:
Где x - отсчет времени, µ - среднее значение сигнала, а σ - стандартное отклонение. Функция, представленная ниже, выполняет расчет за нас.
>>> def z_score(data): >>> mean = data.mean(axis=1, keepdims=True) >>> std = data.std(axis=1, keepdims=True) >>> return (data-mean)/std
Теперь мы можем оценить наши сигналы по оси Z, масштабировать их, как мы делали в предыдущем посте, и добавить среднее значение в нашу матрицу дизайна.
>>> avg = z_score(data_ordered[~np.isnan(map),:]) >>> avg = scale(np.transpose(avg).mean(axis=1)) >>> design_matrix = np.array((constant, predicted_response, avg))
Результирующая матрица дизайна изображена ниже и содержит средний отклик всех вокселей с коэффициентом корреляции выше нашего порога (›0,3).
Если мы теперь снова запустим приведенный выше код с нашей новой матрицей дизайна и построим результирующие карты, мы увидим следующее.
Карта по-прежнему выглядит сопоставимой с нашими предыдущими картами, но размер кластеров увеличился. Также мы видим больше одиночных вокселей выше порога.
Какая карта лучшая?
Теперь это «лучший» результат, чем у наших предыдущих карт? Не обязательно. Очевидно, что наш первый подход к «ожидаемому ответу» - это чрезмерное упрощение. Он не принимает во внимание задержки ответа BOLD-сигнала и другие физиологические процессы, которые здесь играют решающую роль. Так что лучше использовать средний жирный сигнал вокселей выше порога? Опять же, здесь нет правильного или неправильного ответа. Конечно, этот временной курс намного лучше отражает поведение ЖИВОГО сигнала, чем модель простого товарного вагона. Однако подобные подходы также более чувствительны к артефактам. Например, если объект перемещался во время сканирования, это отразится на среднем течении времени, поэтому мы получим карты, отражающие движение, а не реакцию на стимул. Как и в любом другом случае, не существует универсального решения, и вы действительно должны спросить себя, что и почему я делаю то, что делаю.
Снижение шума - пространственное сглаживание
Теперь, прежде чем мы подойдем к концу этой статьи, давайте посмотрим, можем ли мы улучшить уровень шума наших данных. Один из способов сделать это - пространственное сглаживание, которое может быть выполнено путем свертки изображений EPI с гауссовым ядром. Здесь мы не будем вдаваться в подробности того, как построить гауссово ядро и как выполнить двумерную свертку. Если вас интересуют эти вещи, вам следует взглянуть на Блокнот Jupyter для этой статьи или проверить эту ссылку для гауссовского ядра и эту для свертки.
На рисунке выше показан результат сглаживания изображений EPI с помощью ядра Гаусса. В основном он размывает изображение, усредняя соседние воксели с разным весом. Как видите, мы теряем детали, но в то же время уменьшаем шум на изображениях EPI. Если мы сделаем это с каждым изображением EPI в нашем наборе данных, мы получим сглаженный временной курс, на котором мы можем снова запустить GLM и получить сглаженную карту активации ниже.
Как видите, карта выше выглядит чище и менее шумной, чем предыдущие карты. Также увеличился размер нашей активированной области. Но, как мы обсуждали ранее, нам нужно быть осторожными с интерпретацией этой карты. Активированная область на этом графике может быть больше, чем есть на самом деле. Также теряем пространственные детали. С другой стороны, это явно помогает нам с шумом в данных.
Заключение
В этой серии из трех статей мы рассмотрели общую организацию данных МРТ и фМРТ. Мы перешли от визуализации статических изображений МРТ к анализу динамики 4-мерных наборов данных фМРТ с помощью корреляционных карт и общей линейной модели. Наконец, мы уменьшили шум в данных с помощью пространственного сглаживания и увидели кластеры активности в слуховой коре головного мозга во время слуховой стимуляции.
Эти три статьи должны были дать вам представление о том, как организованы данные изображений мозга, полученные в результате экспериментов с фМРТ, и как к ним можно подойти, чтобы получить из них содержательную информацию. Естественно, эта тема сложнее, чем можно выразить в этом коротком формате, но я надеюсь, что вам понравилось читать и вы узнали кое-что.
Если вам нужен полный код этого проекта, вы можете найти его здесь. И, конечно, не стесняйтесь подписываться на меня в Twitter или подключаться через LinkedIn.
Полный код этого проекта можно найти на Github.
