La conectividad cerebral funcional en estado de reposo sin dolor predice la sensibilidad individual al dolor

Consideraciones generales

El diseño del estudio se estableció teniendo muy en cuenta las recientes recomendaciones, requisitos y normas para los biomarcadores de neuroimagen50 (neuromarcadores) y motivado por las siguientes reflexiones.

Mejorar el rendimiento predictivo

Empleamos una línea de preprocesamiento estandarizada para asegurar una sensibilidad óptima del neuromarcador, ya que un tamaño de efecto suficiente es un requisito básico de cualquier utilidad clínica50. Utilizamos una alineación de imágenes de alta precisión, incorporando la anatomía individual al extraer los datos de las series temporales de fMRI. Además, adoptamos recomendaciones y protocolos recientes51 sobre la reducción de artefactos y optimizamos nuestro flujo de trabajo para satisfacer las necesidades especiales del análisis basado en el conectoma. Utilizamos nuestra propia biblioteca de software python de código abierto Pipelines Utilising a Modular Inventory (PUMI, https://github.com/spisakt/PUMI), que se basa en nipype52, un proyecto Python basado en la comunidad que proporciona una interfaz uniforme para el software de neuroimagen existente y, en parte, reutiliza el código de los proyectos de código abierto C-PAC53 y niworkflows54. Se utilizó un enfoque de modelado predictivo (aprendizaje automático) para explotar los ricos datos proporcionados por las redes cerebrales funcionales en estado de reposo y, potencialmente, aprovechar la hiperactividad de la fMRI55.

Evaluación del poder predictivo en condiciones realistas

Utilizamos una estrategia de validación externa prerregistrada, que separaba estrictamente el entrenamiento del modelo y la evaluación del rendimiento. Para el entrenamiento del modelo, utilizamos exclusivamente los datos del Estudio 1. Llevamos a cabo dos subestudios independientes (Estudios 2 y 3) en diferentes centros de investigación, con diferentes equipos y diferente personal de investigación para la validación. Se utilizó una alineación liberal de los entornos de investigación, permitiendo una razonable heterogeneidad en los procedimientos, el equipo, las secuencias de imágenes, el lenguaje de la comunicación participante-investigador a través de los centros de estudio, introduciendo una razonable heterogeneidad en el procedimiento de validación para asegurar la generalizabilidad.

Asegurar que la predicción está impulsada por la señal neural y es específica para la sensibilidad al dolor

Para asegurar que el marcador de sensibilidad al dolor propuesto está realmente impulsado por las señales neurales asociadas a la sensibilidad al dolor, se evaluó la correlación de la puntuación predicha con diversas variables confusoras y validadoras predefinidas (y registradas previamente).

Aseguramos la accesibilidad de los resultados

Aplicamos un prerregistro exhaustivo e hicimos que el código fuente del método fuera de código abierto y estuviera disponible libremente para la comunidad. Además, proporcionamos un contenedor docker independiente de la plataforma y fácil de usar, que ofrece la oportunidad de utilizar nuestro modelo predictivo como producto de investigación50, para obtener predicciones de sensibilidad al dolor listas para usar a partir de cualquier conjunto de datos de imágenes apropiado.

Participantes

Un total de N = 116 voluntarios sanos y jóvenes participaron en tres subestudios. La edad y el sexo de los participantes se indican en la Tabla Suplementaria 1. En el estudio 1 participaron N1 = 39 personas (la misma muestra que en la referencia 8). Se realizó en la Universidad del Ruhr de Bochum (Alemania) por MZ y TSW y se utilizó como muestra de entrenamiento para la predicción de la sensibilidad al dolor basada en el aprendizaje automático y, además, sirvió de base para la validación interna de la predicción. Los estudios 2 y 3 (N2 = 48, N3 = 29) fueron realizados en el Hospital Universitario de Essen (Alemania) por FS y TS y en la Universidad de Szeged (Hungría) por BK y TK, respectivamente, y sirvieron como muestras para la validación externa. Los criterios de inclusión y exclusión fueron en gran medida idénticos en los tres centros y se enumeran en la Tabla 3. Las políticas de reclutamiento y reembolso variaron entre los centros; los participantes recibieron 20 euros/hora en los Estudios 1 y 2 y ningún reembolso en el Estudio 3.

Se consideraron contraindicaciones para la medición de la RM los implantes metálicos, los piercings inamovibles, los marcapasos, los tatuajes en posición de cabeza/cuello, el embarazo o la claustrofobia conocida. Se exigió a los participantes que se abstuvieran de consumir cafeína dos horas antes de los experimentos (excepto en el Estudio 3) y de consumir alcohol el día de la prueba y el día anterior.

El estudio se llevó a cabo de acuerdo con la Declaración de Helsinki, cumple con todas las normas éticas pertinentes para el trabajo con participantes humanos y fue aprobado por los comités de ética locales o nacionales (Números de registro: 4974-14, 18-8020-BO y 057617/2015/OTIG en la Universidad del Ruhr de Bochum, el Hospital Universitario de Essen y el ETT TUKEB de Hungría, respectivamente). Todos los participantes dieron su consentimiento informado por escrito antes de las pruebas.

Las imágenes y las pruebas sensoriales cuantitativas (QST) se realizaron el mismo día en el estudio 1 y con una media de 2-3 días de diferencia en los estudios 2 y 3 (véase la tabla suplementaria 1 para más detalles). La medición de la IRM siempre precedió a la sesión de QST.

Medidas de la IRM funcional

Se adquirieron mediciones de la IRMf de alta resolución anatómica y de estado de reposo con los ojos abiertos de todos los participantes. Los parámetros de escaneo (incluyendo el equipo) variaron entre los centros y se enumeran en la Tabla 4. Durante las mediciones, se indicó a los participantes que permanecieran quietos y relajados, sin dormirse, y que evitaran cualquier movimiento. Se utilizaron almohadillas de espuma y, en los Estudios 1 y 2, almohadas neumáticas para restringir los movimientos de la cabeza. Todas las mediciones anatómicas de la RM se examinaron en busca de hallazgos incidentales.

Medidas-QST

Los umbrales de dolor de calor (HPT), frío (CPT) y mecánico (MPT) se adquirieron de acuerdo con el protocolo QST28. Los umbrales de detección de calor (WDT), frío (CDT) y, en el Estudio 2 y el Estudio 3, mecánico (MDT) se obtuvieron como medidas de control adicionales. Todas las mediciones sensoriales se obtuvieron del antebrazo izquierdo palmar, proximal a la cresta de la muñeca. En el marco del QST, los umbrales térmicos se determinan mediante un método de límites. Para ello, se aplicaron temperaturas crecientes y decrecientes a la piel con un estimulador térmico MSA (Somedic, Hörby, Suecia) en el Estudio 1 y estimuladores térmicos Pathway (Medoc Ltd., Ramat Yishai, Israel) en los Estudios 2 y 3. En todos los estudios, se utilizaron termodos ATS en una superficie cutánea de 30 × 30 mm, con una temperatura de referencia de 32 °C. Se indicó a los participantes que indicaran la aparición del dolor pulsando un botón. Para todos los umbrales térmicos se realizaron 6, en lugar de 3 (como en el protocolo original)28, repeticiones del estímulo para reducir la varianza entre sujetos. Además, la primera medición se descartó del análisis como estímulo de prueba. La HPT y la CPT se calcularon como las medias aritméticas de las cinco temperaturas umbral restantes. Las MPT y MDT se determinaron mediante un método escalonado. Se aplicaron cinco trenes crecientes y cinco decrecientes de estímulos de pinchazos (MRC Systems, Heidelberg, Alemania) en la palma del antebrazo izquierdo de forma alternada, mientras que se instruía al participante para que clasificara los estímulos como nocivos o nocivos. El umbral de detección mecánica se evaluó de forma análoga a los estímulos de filamentos de von Frey. El MPT y el MDT se calcularon como la fuerza media geométrica transformada logarítmicamente determinada en cinco recorridos de umbralización ascendente y descendente.

Medidas adicionales

Antes de todas las mediciones se registraron la edad, el sexo, la altura y el peso autodeclarados y, en el caso de las mujeres, la fecha del primer día de la última menstruación y el uso de anticonceptivos. Además, se registró el consumo semanal de alcohol autodeclarado y el nivel de educación (escuela primaria, escuela secundaria, universidad) para los Estudios 1 y 2. Antes del QST, los participantes rellenaron el Cuestionario de Sensibilidad al Dolor (PSQ)56, la Escala de Catastrofización del Dolor (PCS)57, el Inventario de Ansiedad Estado-Rasgo (STAI)58, la versión corta en alemán de la Escala de Depresión (ADS-K, Centro de Estudios Epidemiológicos)59 y, además, en los Estudios 2 y 3, el Índice de Calidad del Sueño de Pittsburgh (PSQI)60 y el cuestionario de estrés percibido (PSQ20)61. En los Estudios 2 y 3, se midió la presión arterial tanto antes de la RM como de las mediciones del QST. Además, para la Muestra 1, se disponía de los valores T50 de un experimento paralelo realizado el día anterior a la prueba de RMNf. La T50 representa la temperatura (en °C) necesaria para inducir una calificación de dolor por calor de 50 (en una escala que va de 0, ningún dolor, a 100, dolor insoportable). Los valores de T50 se obtuvieron a partir de una interpolación no lineal (polinomio de segundo orden) de las puntuaciones obtenidas en respuesta a 15 estímulos tónicos de dolor por calor (duración: 16 s) entre 42,5 °C y 48 °C, presentados en forma de rejilla de búsqueda pseudoaleatoria.

Cálculo de la sensibilidad al dolor

La variable objetivo para la predicción fue una única medida compuesta de la sensibilidad individual al dolor que resumía la HPT, la CPT y la MPT, tal como se define en la ref. 8.

En el Estudio 1, la HPT, la CPT y la MPT se transformaron en Z (media centrada y estandarizada) y la HPT, así como la MPT, se invirtieron (multiplicadas por -1), de modo que los valores Z más altos denotaron una mayor sensibilidad al dolor. A continuación, se calculó la media aritmética de las variables transformadas en Z para cada participante y se definió como puntuación de sensibilidad al dolor. En los Estudios 2 y 3, se aplicó el mismo procedimiento, excepto que la transformación Z se basó en la media de la población y la desviación estándar del Estudio 1, para garantizar que se utilizara la misma escala en todos los estudios. Los valores extremos de QST se definieron utilizando los percentiles normativos del 95% reportados en la ref. 28; los participantes que mostraron valores extremos de HPT, CPT o MPT en al menos dos de las tres modalidades fueron excluidos. Esta selección dio como resultado la exclusión de 0, 3 y 2 participantes en las Muestras 1, 2 y 3, respectivamente (Tabla Suplementaria 2).

Preprocesamiento de la RMNf

Como la conectividad funcional basada en la RMNf es susceptible a los artefactos de movimiento en el escáner62,63, el preprocesamiento y la limpieza de la señal adecuados son clave para una predicción exitosa basada en la conectividad. Los datos de IRM funcional en estado de reposo se preprocesaron de forma idéntica en los tres estudios. El flujo de trabajo aplicado, basado en el nipype, se representa en la Fig. 1 suplementaria. Se utilizó software de neuroimagen de terceros, código adaptado de las herramientas de software C-PAC53 y niworkflows54, y rutinas propias de python.

La extracción del cerebro de las imágenes anatómicas y estructurales, así como la segmentación del tejido de las imágenes anatómicas se realizó con FSL bet y fast64. Las imágenes anatómicas se co-registraron lineal y no linealmente a la plantilla cerebral estándar MNI152 de 1mm de resolución con ANTs65 (ver https://gist.github.com/spisakt/0caa7ec4bc18d3ed736d3a4e49da7415 para el código fuente).

Las imágenes funcionales se co-registraron a las imágenes anatómicas con la técnica de registro basada en los límites de FSL flirt. Todas las transformaciones resultantes se guardaron para su uso posterior. El preprocesamiento de las imágenes funcionales se realizó en el espacio de imagen nativo, sin remuestreo. La corrección del movimiento basada en la realineación se realizó con FSL mcflirt. Las seis estimaciones de movimiento de la cabeza resultantes (3 rotaciones, 3 traslaciones), sus versiones al cuadrado, sus derivados y los derivados al cuadrado (conocidos como la expansión Friston-2466) se calcularon y guardaron para la corrección de molestias. Además, el movimiento de la cabeza se resumió como series temporales de desplazamiento en el marco (FD), de acuerdo con el método de Power63, para ser utilizado en la censura y exclusión de datos. Tras la corrección del movimiento, los valores atípicos (por ejemplo, los picos de movimiento) en los datos de las series temporales se atenuaron utilizando AFNI despike67. La unión de los mapas de materia blanca erosionados y las máscaras de los ventrículos se transformaron al espacio funcional nativo y se utilizaron para extraer la señal de ruido para la corrección anatómica CompCor68.

En un paso de regresión de molestia, se eliminaron 6 parámetros CompCor (los 6 primeros componentes principales de las series temporales de ruido), los parámetros de movimiento Friston-24 y la tendencia lineal de los datos de las series temporales con un modelo lineal general. En los datos residuales, se realizó un filtrado temporal de paso de banda con 3DBandpass de AFNI para conservar la banda de frecuencia de 0,008-0,08 Hz. Se espera que el uso previo del despike de AFNI atenúe el aliasing de los artefactos de movimiento residuales en los fotogramas de tiempo vecinos durante el filtrado de paso de banda69. Para atenuar aún más el impacto de los artefactos de movimiento, se eliminaron de los datos los fotogramas de tiempo potencialmente contaminados por el movimiento, definidos por un umbral conservador de FD > 0,15 mm (lo que se conoce como depuración de los datos)70 . Los participantes fueron excluidos de los análisis posteriores si la media de la FD superaba los 0,15 mm, o si se eliminaba más del 30% de los fotogramas. Esto dio lugar a la exclusión de 4, 8 y 7 participantes en las muestras 1, 2 y 3, respectivamente (tabla suplementaria 2). Se realizó un control de calidad (comprobación de registro, diagramas de alfombra, véase, por ejemplo, las figuras suplementarias 2-4) durante todo el flujo de trabajo.

Análisis de conectividad funcional

La versión de 122 parcelas del atlas cerebral funcional multirresolución MIST71 y las máscaras de materia gris obtenidas de la imagen anatómica se transformaron al espacio funcional nativo. Este atlas (construido con el método BASC, es decir, el análisis bootstrap de clusters estables) ha demostrado recientemente un buen rendimiento en la modelización predictiva basada en la conectividad72. Las regiones del atlas en espacio nativo se enmascararon con las máscaras de materia gris que se obtuvieron de la imagen anatómica y se transformaron al espacio funcional previamente. Con esta técnica de individualización del atlas, la señal regional final se originará -con una alta probabilidad- a partir de vóxeles de materia gris para cada sujeto (que comprobamos cuidadosamente de forma manual para todos los sujetos), mientras que con el método convencional, se incluye una proporción variable de vóxeles de materia gris y blanca para cada sujeto. Por lo tanto, se espera que la introducción de información del proceso de segmentación del tejido disminuya la variabilidad de un sujeto a otro (ver ejemplos en la Fig. 5). Las series de tiempo de los vóxeles se promediaron en estas regiones MIST individualizadas y, junto con la señal media de materia gris, se retuvieron para el análisis de conectividad basado en gráficos.

Las series de tiempo regionales se ordenaron en módulos funcionales a gran escala (definidos por el atlas MIST de 7 parcelas) para su visualización (Fig. 1). Se calculó la correlación parcial entre todos los pares de regiones (y la materia gris global), tal y como se implementó en el módulo python nilearn73. Se utilizaron correlaciones parciales, en lugar de simples, para descartar la conectividad indirecta74. Nuestro enfoque de modelado de gráficos garantiza que la señal de materia gris global se maneja como un factor de confusión durante el cálculo de los coeficientes de correlación parcial pero, al mismo tiempo, también se considera una señal de interés, ya que puede representar procesos relacionados con la vigilancia75. Los coeficientes de correlación parcial se organizaron en matrices de conectividad simétrica de 123 por 123 (122 regiones + señal global de materia gris). El triángulo superior de estas matrices se utilizó como espacio de características para el modelado predictivo basado en el aprendizaje automático.

Entrenamiento y validación del modelo predictivo

Los datos de conectividad funcional de todo el cerebro en estado de reposo del estudio 1 (N1 = 35, después de todas las exclusiones, como en la ref. 8, Tabla Suplementaria 2) se utilizaron como espacio de características de entrada (P = 7503 características por participante) para predecir las puntuaciones individuales de sensibilidad al dolor, lo que llevó a una configuración de P grande y N pequeña.

Construimos una línea de aprendizaje automático (https://github.com/spisakt/RPN-signature/blob/master/PAINTeR/model.py) en scikit-learn76, que consistía en un escalado robusto de las características (elimina la mediana y escala con los cuantiles de los datos), la preselección de las características77, la selección de las K mejores características con las relaciones más fuertes con la variable objetivo y un modelo de regresión de red elástica78 (un modelo lineal con las formas L1 y L2 combinadas como regularizador). El uso de la red elástica fue una decisión tomada antes del análisis. Nuestra principal motivación para elegir la red elástica fue que permite optimizar la dispersión (regularización L1 frente a L2) como hiperparámetro, de modo que no tuvimos que hacer ninguna suposición a priori sobre la dispersión de la verdad básica discriminatoria (véase la ref. 79 para la justificación). En resumen, los hiperparámetros libres de la línea de aprendizaje automático fueron el número de características preseleccionadas (K), la relación de la regularización L1/L2 y el peso (alfa) de la regularización. Los hiperparámetros se optimizaron con un procedimiento de búsqueda en cuadrícula y un error medio cuadrático negativo como función de coste. Los valores de K oscilaron entre 10 y 200 con incrementos de 5, e incluyeron la relación L1/L2 para alfa. La optimización de los hiperparámetros se llevó a cabo con una validación cruzada de un participante (fase de validación interna). La validación cruzada incorporó la tubería completa de aprendizaje automático para evitar la introducción de dependencias entre las muestras de entrenamiento y de prueba. Hay que tener en cuenta que el preprocesamiento fMRI fue independiente entre los sujetos, por lo que no se incluyó en la validación cruzada. Los hiperparámetros óptimos fueron K = 25, L1/L2-ratio = 0.999 y alpha = 0.005.

La validación externa se realizó aplicando la firma RPN en los datos fMRI de los Estudios 2 y 3 (N2 = 37, N3 = 19, después de las exclusiones, Tabla Suplementaria 2), simplemente aplicando la transformación de características (escalado) obtenida en la Muestra 1 y luego calculando el producto punto entre las matrices de conectividad individuales y los pesos de características no nulos obtenidos en la Muestra 1. Las predicciones resultantes se compararon con las puntuaciones de sensibilidad al dolor basadas en QST calculando el error medio absoluto (MAE), el error medio cuadrático (MSE) y la varianza explicada. Se obtuvieron valores p basados en permutaciones para las tres medidas, utilizando el paquete mlxtend python. Además, se utilizó el bootstrap con cobertura condicional80 para obtener los valores p de los pesos de conectividad predictiva para ayudar a la interpretación. Construimos 10000 muestras bootstrap (con reemplazo), con un tamaño igual al de la muestra original, consistente en datos emparejados de cerebro y resultado. El modelo predictivo con los hiperparámetros óptimos se ajustó a cada muestra. Los valores P no corregidos se calcularon para cada conexión seleccionada basándose en la proporción de pesos por debajo o por encima de cero, como en, por ejemplo, la ref. 30. Nótese que la interpretación de estos valores p e intervalos de confianza (Tabla Suplementaria 4) sigue siendo limitada, ya que están condicionados al procedimiento de selección de características.

Análisis de factores de confusión

Para explorar posibles variables de confusión, las puntuaciones de sensibilidad al dolor predichas (o las predicciones validadas de forma cruzada en el caso de la Muestra 1) se contrastaron con la media y la mediana de la FD, el porcentaje de volúmenes fregados y la presión arterial sistólica y diastólica antes de la medición de la IRM y la QST (ya que se informó anteriormente81 de que la presión arterial estaba asociada a la sensibilidad al dolor mecánico), el tiempo de retraso entre la RM y la prueba QST (para comprobar la estabilidad temporal de la predicción), la edad, el sexo, el IMC, el número de días transcurridos desde el primer día de la última menstruación, el consumo de alcohol (unidades/semana), el nivel de educación, la ansiedad estado y rasgo (STAI) puntuación de los síntomas depresivos (ADS-K), sensibilidad al dolor autodeclarada (PSQ) y catastrofismo del dolor (PCS), estrés percibido (PSQ20), calidad del sueño (PSQI) y umbrales de detección de QST no nocivos (CDT, WDT y MDT, cuando estaban disponibles). Además, en el Estudio 1 las predicciones se compararon con los valores T50 y los niveles de GABA y Glutamato/Glutamina basados en la espectroscopia de RMN en las regiones cerebrales que procesan el dolor (véase la ref. 8 para más detalles). Las asociaciones se probaron con modelos lineales basados en permutaciones.

Visualización de la red predictiva

Las conexiones interregionales predictivas resaltadas por los coeficientes de regresión no nulos de la firma RPN se mostraron como un gráfico de cinta utilizando el paquete R circlize (Fig. 3). Las correspondientes máscaras de regiones cerebrales individualizadas se transformaron de nuevo al espacio estándar para crear un mapa de probabilidad regional específico para el estudio (que refleja la precisión del corregistro y la variabilidad individual en la morfología). Los mapas de probabilidad se multiplicaron por la suma de los coeficientes de regresión correspondientes para crear un mapa de fuerza predictiva regional, que luego se visualizó con FSLeyes y MRIcroGL.

(https://www.mccauslandcenter.sc.edu/mricrogl) (Fig. 3). El análisis de la implicación de la red en estado de reposo a gran escala (según la definición del atlas cerebral MIST71) se realizó resumiendo y transformando en Z los valores de los vóxeles en las siete regiones de interés. El gráfico polar se realizó con el paquete R ggplot2.

Disponibilidad del software

Las puntuaciones de la firma RPN pueden calcularse basándose en conjuntos de datos estructurales y funcionales en estado de reposo mediante la herramienta de software con el mismo nombre. La herramienta de software RPN-signature consiste en la tubería de procesamiento de IRM descrita y el modelo predictivo funcional basado en el conectoma. Está disponible como código fuente en https://github.com/spisakt/RPN-signature. Como el software sigue la estructura de datos de imágenes cerebrales (BIDS)82 y la especificación BIDS-App, proporciona una interfaz de línea de comandos estándar y se basa en la tecnología Docker. La imagen Docker está depositada en Docker Hub: (https://cloud.docker.com/repository/docker/tspisak/rpn-signature) y no depende de ningún software fuera de la imagen del contenedor. Esto, junto con el desarrollo totalmente transparente basado en la integración continua y el etiquetado y versionado automatizados, mejora la disponibilidad del software y apoya la reproducibilidad de los resultados de la firma RPN.

Resumen del informe

Más información sobre el diseño de la investigación está disponible en el resumen del informe de Nature Research vinculado a este artículo.