Conjuntos de datos
En este estudio utilizamos un enfoque meta-analítico de los estudios de resonancia magnética funcional descritos por Yarkoni et al.18 (http://Neurosynth.org). Descargamos la base de datos Neurosynth que contenía 3107 mapas funcionales reversos sin umbral y los detalles de 11.406 fuentes bibliográficas a partir del 25 de septiembre de 2017.
Los datos del conectoma estructural se derivaron del conjunto de datos de imágenes ponderadas por difusión de 163 participantes adquiridos a 7 Tesla por el equipo del Proyecto Conectoma Humano62 (http://www.humanconnectome.org/study/hcp-young-adult/) (Consorcio WU-Minn; Investigadores principales: David Van Essen y Kamil Ugurbil; 1U54MH091657). Este proyecto fue financiado por los 16 institutos y centros de los NIH que apoyan el Blueprint for Neuroscience Research de los NIH, y por el McDonnell Center for Systems Neuroscience de la Universidad de Washington.
Preprocesamiento de los datos de Neurosynth
Dos investigadores (V.R.K. y M.T.S.) actuaron como jueces, seleccionando los términos que, en su opinión, estaban relacionados con procesos cognitivos específicos. El procedimiento de selección consistió en dos etapas. En la primera, los jueces hicieron su selección de forma independiente. Se excluyeron sistemáticamente los términos anatómicos cerebrales (por ejemplo, «red de saliencia»), psiquiátricos (por ejemplo, «esquizofrenia») y patológicos (por ejemplo, «alzheimer»). Los dos jueces se pusieron de acuerdo en 422 términos como relacionados con los procesos cognitivos, así como en 2309 términos no relacionados que debían ser descartados (88% de reproducibilidad). Para el resto de términos, los jueces tomaron su decisión conjuntamente. Al final, se seleccionaron 590 términos cognitivos para el estudio.
En el presente análisis, corregimos las diferencias anatómicas entre los hemisferios izquierdo y derecho para centrarnos en las asimetrías funcionales. Dado que los mapas funcionales de Neurosynth se proporcionan en el espacio de plantilla MNI estándar de 2 mm, que no es simétrico, co-registramos de forma no lineal la plantilla MNI a una plantilla simétrica MNI, disponible en http://www.bic.mni.mcgill.ca/ServicesAtlases/ICBM152NLin2009, utilizando la tubería de normalización difeomórfica simétrica de Greedy (GreedySyN) distribuida con las herramientas de normalización avanzadas (ANTs, http://stnava.github.io/ANTs/)63. La plantilla simétrica se redujo a un tamaño de vóxel de 2 mm para que coincidiera con las dimensiones de vóxel de la plantilla estándar. La transformación estimada entre los espacios MNI no simétricos y simétricos se aplicó entonces a todos los mapas funcionales.
Se realizaron los siguientes pasos para obtener índices de lateralización para cada mapa funcional tras su corregistro con la plantilla simétrica. En primer lugar, dividimos los mapas funcionales en las partes del hemisferio izquierdo y derecho y suavizamos los mapas resultantes utilizando un filtro gaussiano de 6 mm FWHM. A continuación, invertimos los mapas del hemisferio izquierdo y los restamos de los mapas del hemisferio derecho no invertidos para obtener mapas de índices de lateralidad (LI) (véase24 para un enfoque similar). Los valores positivos y negativos en estos mapas significarían una mayor evidencia meta-analítica para, respectivamente, la lateralización derecha e izquierda de la función asociada a un término.
Preprocesamiento de los datos del conectoma estructural
Los parámetros de escaneo se han descrito previamente en Vu et al.62. En resumen, cada imagen ponderada por difusión consistió en un total de 132 cortes casi axiales adquiridos con un factor de aceleración de 3 (ref. 64), resolución isotrópica (1,05 mm3) y cobertura de toda la cabeza con un TE de 71,2 ms y con un TR de 7000 ms. En cada posición de corte, se adquirieron imágenes ponderadas en difusión con 65 gradientes distribuidos uniformemente en múltiples casillas del espacio Q65 y 6 imágenes sin gradiente de difusión aplicado. Esta adquisición se repitió cuatro veces con un valor b de 1000 y 2000 s mm-2 en pares con direcciones de codificación de fase de izquierda a derecha y de derecha a izquierda. Se aplicó a los datos la línea de preprocesamiento HCP por defecto (v3.19.0)66,67. En resumen, el campo de fuera de resonancia inducido por la susceptibilidad se estimó a partir de pares de imágenes con gradiente de difusión aplicado con distorsiones que iban en direcciones opuestas68 y se corrigió para todo el conjunto de datos ponderados por difusión utilizando TOPUP69. Posteriormente, se corrigieron el movimiento y la distorsión geométrica utilizando la herramienta EDDY implementada en FSL.
La caja de herramientas ExploreDTI para Matlab (http://www.exploredti.com70,71) se ha utilizado para extraer estimaciones de la fracción de agua axonal28. A continuación, descartamos los volúmenes con un valor b de 1000 s mm-2 y posteriormente se realizó una tractografía determinista de todo el cerebro en el espacio DWI nativo utilizando el software StarTrack (https://www.mr-startrack.com). Se aplicó un algoritmo de Richardson-Lucy amortiguado para las deconvoluciones esféricas72. Se adoptó una respuesta de fibra fija correspondiente a un factor de forma de α = 1,5 × 10-3 mm2 s-1, junto con el parámetro de amortiguación geométrica de 8. Se ejecutaron 200 iteraciones del algoritmo. El umbral absoluto se definió como tres veces la distribución esférica de la orientación de las fibras (FOD) de un vóxel isotrópico de materia gris y el umbral relativo como el 8% de la amplitud máxima de la FOD73. Se utilizó un algoritmo de Euler modificado74 para realizar la tractografía de líneas de corriente de todo el cerebro, con un umbral de ángulo de 35°, un tamaño de paso de 0,5 mm y una longitud de línea de corriente mínima de 15 mm.
Corregimos los datos del conectoma estructural al espacio estándar MNI de 2 mm utilizando los siguientes pasos: primero, la tractografía de líneas de corriente de todo el cerebro se convirtió en volúmenes de densidad de líneas de corriente donde las intensidades correspondían al número de líneas de corriente que cruzaban cada vóxel. En segundo lugar, se generó una plantilla específica para el estudio de los volúmenes de densidad de las líneas de corriente utilizando la línea de normalización difeomórfica simétrica de Greedy (GreedySyN) distribuida con ANTs. Esto proporcionó una plantilla media de los volúmenes de densidad de las líneas de corriente para todos los sujetos. La plantilla se co-registró con una plantilla estándar de 2 mm MNI152 utilizando la herramienta flirt implementada en FSL. Este paso produjo una plantilla de densidad de líneas de corriente en el espacio MNI152. En tercer lugar, los volúmenes de densidad de líneas de corriente individuales se registraron en la plantilla de densidad de líneas de corriente en la plantilla del espacio MNI152 y se aplicó la misma transformación a la tractografía de líneas de corriente individuales de todo el cerebro utilizando la herramienta trackmath distribuida con el paquete de software Tract Querier75, y a los mapas de fracción de agua axonal, utilizando ANTs GreedySyn. Este paso produjo una tractografía de líneas de flujo de todo el cerebro y mapas de fracción de agua axonal en el espacio estándar MNI152.
Determinación de regiones funcionalmente lateralizadas
En estos análisis, completados en dos pasos, pensamos en identificar las regiones con lateralización funcional significativa. Véase la figura suplementaria 5. En el primer paso, abordamos la redundancia preservando la riqueza de los datos de Neurosynth. Por ejemplo, muchos términos seleccionados estaban relacionados como formas singulares y plurales de la misma palabra (por ejemplo, «forma visual» y «formas visuales») y, por tanto, es probable que sus mapas sean muy similares. Para ello, redujimos la dimensionalidad de los datos mediante un análisis de componentes principales (PC) con rotación varimax implementado en SPSS (SPSS, Chicago, IL) con los mapas de LI como entradas76,77,78. Tras un análisis de componentes principales estándar, que implica la eigendecomposición de la matriz de covarianza, 171 componentes ortogonales extraídos con valores propios superiores a la gran media se sometieron al procedimiento de rotación varimax utilizando el criterio de normalización de Kaiser79, con un máximo de 1000 iteraciones para la convergencia. De este modo, se explicó el 72,6% de la varianza de los datos. La distribución de las cargas a lo largo de los componentes principales con rotación varimax suele ser asimétrica y sólo unos pocos elementos reciben grandes cargas. Por lo tanto, con el fin de discutir los resultados, los componentes fueron etiquetados de acuerdo con el término(s) con las mayores cargas (Tabla Suplementaria 3).
En el segundo paso, se empleó la modelización lineal general para identificar los vóxeles con una lateralización significativa asociada a un componente particular. En este análisis, los componentes principales se utilizaron como un conjunto de predictores para ajustar los mapas de LI y obtener mapas beta, es decir, mapas espaciales de componentes. Se realizó la prueba de permutación para identificar las regiones significativamente lateralizadas. Dado que la rotación varimax puede imponer algunas correlaciones entre las columnas de la matriz de componentes principales, realizamos permutaciones en las filas de la matriz sin rotar, aplicando posteriormente la rotación de componentes y calculando un mapa aleatorio en cada permutación del mismo modo que se hizo para los componentes principales reales. Este procedimiento nos permitió imitar la estructura correlacional de los datos no permutados y proporcionar una prueba de significación más sólida. Para tener en cuenta las comparaciones múltiples, se utilizó el enfoque de la estadística máxima, según el cual los valores del mapa espacial para los componentes principales reales se compararon con el valor máximo (positivo o negativo) de todo un mapa aleatorio en cada permutación. Se realizaron cinco mil permutaciones. Se consideró que los vóxeles mostraban una lateralización significativa si satisfacían simultáneamente dos criterios (1) los valores de su mapa espacial eran, en el 97,5% de los casos, superiores o inferiores, respectivamente, a los valores máximos positivos y negativos obtenidos a través de las permutaciones (es decir, p < 0,05, de dos colas y con corrección de FWE); (2) formaban un grupo de al menos 20 vóxeles. El segundo criterio se utilizó para excluir los efectos pequeños y posiblemente espurios observados en un pequeño número de vóxeles.
Incorporación multivariante de los mapas de lateralización
Con el fin de caracterizar una estructura de baja dimensión de la lateralización funcional del cerebro, se realizó una incrustación espectral de los mapas de LI utilizando la eigendecomposición del laplaciano normalizado del gráfico de la matriz de similitud80. El método buscaba descubrir características geométricas en las similitudes entre los mapas de lateralización convirtiendo estas similitudes en distancias entre los mapas de lateralización en el espacio incrustado (cuanto mayor sea la similitud entre los perfiles de lateralización, menor será la distancia). En este caso, nos concentramos únicamente en las varianzas que fueron contabilizadas por los 171 componentes analizados en el presente estudio. Para ello, los mapas de LI se «desnotizaron», en el sentido de que se reconstruyeron como el producto matricial de 171 componentes y sus mapas espaciales. Cada elemento de la matriz de similitud se calculó como un producto de puntos tomado para un par de mapas de LI «desnoizados» a través de todos los vóxeles (es decir, un elemento de la matriz de similitud era una suma de productos de los valores de los vóxeles para un par de mapas). Los valores negativos se redujeron a cero para permitir la estimabilidad. Las dimensiones de incrustación se ordenaron según sus valores propios, de menor a mayor. Se eliminó la primera dimensión no informativa asociada a un valor propio cero. En el análisis se trató de determinar si existe una estructura en una representación de baja dimensión de los datos, concretamente la triangularidad estructural de los datos, y si existe, en cuántas dimensiones se conserva esta estructura (para el gráfico de valores propios, véase la figura suplementaria 6). La estructura triangular se cuantificó como un ratio t, es decir, un ratio entre el área del casco convexo que engloba todos los puntos en el espacio incrustado y un triángulo envolvente de área mínima27. Estos valores se compararon con los cocientes t de los mapas LI aleatorios. Estos mapas aleatorios se obtuvieron generando 2000 conjuntos de 590 mapas aleatorios mediante la permutación del orden de los vóxeles. Para cada conjunto, se calcularon mapas LI aleatorios para cada par y luego se sometieron a un análisis varimax con el número de componentes principales = 171. El procedimiento de incrustación fue idéntico al aplicado a los mapas LI no aleatorios. La amplitud dimensional de la organización triangular se evaluó comprobando si el cociente t de los mapas de LI no aleatorios era mayor que los cocientes t de los mapas de LI aleatorios en cada subespacio bidimensional de incrustación (p < 0,05, corregido por Bonferroni). La etiqueta de los ejes se definió ad-hoc según uno o varios términos situados en los vértices del triángulo. Los mapas de arquetipos se aproximaron mediante un enfoque de regresión múltiple. En primer lugar, se realizó una regresión de los valores de cada vóxel en los mapas de LI «denotados» sobre las coordenadas de los mapas correspondientes en las primeras 171 dimensiones del espacio incrustado (es decir, coincidiendo con el número de componentes utilizados para la «denotada»). Esto proporcionó una contribución estimada de cada dimensión incrustada al índice de lateralización. A continuación, obtuvimos los mapas de arquetipos evaluando los coeficientes de regresión para las dimensiones en las que se observaba la estructura triangular en las ubicaciones estimadas de los arquetipos (es decir, en los vértices del «simplex» – triangular multidimensional).
Determinación de las regiones no lateralizadas
En los siguientes análisis contrastamos los perfiles de conectividad de las regiones lateralizadas con las regiones que no muestran una lateralización significativa pero que, sin embargo, muestran una implicación significativa al menos en una función. Estas últimas se identificaron repitiendo los análisis descritos en la sección «Determinación de las regiones funcionalmente lateralizadas» con los mapas funcionales originales de Neurosynth como entradas. Véase la figura suplementaria 7. De este modo se obtuvieron 69 componentes, que representan el 70,6% de la varianza. Para una mayor comparabilidad, el análisis se realizó en el espacio simétrico y para los hemisferios izquierdo y derecho por separado. Se consideró que los vóxeles no tenían una lateralización significativa si cumplían los siguientes criterios (1) superaban el umbral de significación para al menos un componente y un hemisferio; (2) no se solapaban con vóxeles lateralizados; y (3) eran homólogos de los vóxeles que cumplían los criterios (1) y (2) en el hemisferio opuesto. En el resto del texto se utilizó el término abreviado «regiones no lateralizadas» para denominar los vóxeles sin lateralización significativa. Esto proporciona un contraste conservador para las regiones lateralizadas porque, en virtud del enfoque estadístico frecuentista, las regiones no lateralizadas también incluirían vóxeles que demuestran una lateralización considerable pero que no cumplen los criterios estadísticos de significación utilizados en el estudio. El número de vóxeles no lateralizados fue 3,6 veces mayor que el número de vóxeles lateralizados.
Medidas de la fuerza de la conectividad para las relaciones estructura-función
Se utilizaron los siguientes pasos para las relaciones estructura-función. Primero, combinamos los mapas espaciales de los vóxeles significativamente lateralizados, independientemente de la polaridad izquierda y derecha de la lateralización. En segundo lugar, transformamos el mapa combinado de nuevo en el espacio MNI regular para un análisis conjunto con la información de difusión utilizando una inversa de las deformaciones de plantilla MNI no simétricas a MNI simétricas estimadas anteriormente. Por último, proyectamos el mapa combinado sobre el límite de la materia blanca de la plantilla MNI no simétrica en cada hemisferio y, posteriormente, seleccionamos la tractografía desde estos vóxeles hasta el cuerpo calloso. Los mismos procedimientos se aplicaron a los mapas de las regiones no lateralizadas.
Se analizaron dos medidas para la fuerza de la conectividad estructural interhemisférica. La primera medida, microestructural, se refería a la fracción de agua axonal, promediada entre los participantes de la muestra HCP, en los vóxeles del cuerpo calloso que fueron alcanzados por las líneas de corriente de las regiones lateralizadas (o no lateralizadas) seleccionadas. La segunda medida macroestructural de la conectividad se definió en términos de replicabilidad de la conexión30 entre los vóxeles del cerebro y el cuerpo calloso, es decir, como la proporción de participantes en los que existe una conexión entre los vóxeles del cerebro y el cuerpo calloso respecto al tamaño total de la muestra HCP. Nos referiremos a esta medida como «probabilidad de conexión» para abreviar.
Comparación de la conectividad entre regiones lateralizadas y no lateralizadas
La comparación de la conectividad entre regiones lateralizadas y no lateralizadas se realizó mediante el muestreo de subconjuntos de vóxeles (sin reemplazo) de los pools de vóxeles corticales lateralizados y no lateralizados. Una muestra de cada grupo era igual al 5% del número total de vóxeles en ese grupo (es decir, asegurando que la frecuencia espacial dentro del grupo de las muestras extraídas era igual entre los grupos). Para cada subconjunto calculamos un valor medio para la probabilidad de conexión y una media ponderada para la fracción de agua axonal callosa, donde un peso para un vóxel se dio como una replicabilidad de conexión entre este vóxel y cualquier vóxel en un subconjunto muestreado. Un valor negativo indicaría una conectividad más débil de los vóxeles lateralizados. Las distribuciones de la diferencia en las medidas de conectividad entre las regiones corticales lateralizadas y no lateralizadas se obtuvieron repitiendo el procedimiento 1000 veces y para cada hemisferio por separado.
Análisis de la dominancia hemisférica
El grado de dominancia funcional hemisférica se evaluó en radianes como arctangente de la relación entre las fuerzas de activación en dos hemisferios. Se restó Pi/4 de este valor para asegurar que la magnitud absoluta de este valor aumenta si la activación de la tarea es unilateral y disminuye si ambos hemisferios muestran niveles comparables de actividad en la tarea. Dado que es posible un solapamiento espacial parcial entre las regiones lateralizadas asociadas a los diferentes componentes, en los análisis elegimos los valores de dominancia asociados a los componentes que arrojaban la mayor puntuación z en un vóxel concreto. Para obtener una estimación robusta de la relación entre la dominancia hemisférica y la fuerza de la conectividad inter-hemisférica, los vóxeles fueron divididos por las probabilidades de conexión de tal manera que la anchura del vóxel más pequeño era de un tamaño igual a 1/163 y aumentaba con la probabilidad de conexión (dada por la función logspace en Matlab). Este procedimiento se utilizó para compensar parcialmente el hecho de que sólo un número muy limitado de vóxeles tenía una alta probabilidad de conexión con el cuerpo calloso, mientras que la mayoría se caracterizaba por valores pequeños. También se estimó la actividad media del vóxel entre los hemisferios izquierdo y derecho (es decir, (actividad del hemisferio izquierdo + derecho)/2) y se utilizó como covariable de no interés en los análisis que examinaban la relación entre la dominancia hemisférica y otras medidas.