Conectividade cerebral funcional em repouso sem dor prevê a sensibilidade individual à dor

Considerações gerais

O desenho do estudo foi estabelecido com cuidadosa consideração das recomendações, requisitos e padrões recentes para biomarcadores de neuroimagem50 (neuromarcadores) e motivado pelas seguintes considerações

Maximizar o desempenho preditivo

Empregamos um pipeline de pré-processamento padronizado para assegurar a sensibilidade ótima do neuromarcador, pois o tamanho do efeito suficiente é um requisito básico de qualquer utilidade clínica50. Utilizamos alinhamento de imagem de alta precisão, incorporando anatomia individual na extração de dados de fMRI timeseries. Além disso, adotamos recomendações e protocolos recentes51 relativos à redução de artefatos e otimizamos nosso fluxo de trabalho para atender às necessidades especiais da análise baseada em conexão. Utilizámos a nossa biblioteca de software python de código aberto desenvolvida internamente, Pipelines Utilising a Modular Inventory (PUMI, https://github.com/spisakt/PUMI), que é baseada no nipype52, um projecto Python baseado na comunidade que fornece uma interface uniforme ao software de neuroimagem existente e, em parte, o código reutilizado do C-PAC53 e dos projectos de código aberto niworkflows54. Uma abordagem de modelagem preditiva (aprendizagem de máquina) foi utilizada para explorar os dados ricos fornecidos pelas redes cerebrais funcionais em estado de repouso e, potencialmente, tirar vantagem da hiperacuidade do fMRI55.

Avaliando o poder preditivo sob condições realistas

Usamos uma estratégia de validação externa pré-registada, que separou rigorosamente o treinamento modelo e a avaliação de desempenho. Para o treinamento com modelo, utilizamos exclusivamente dados do Estudo 1. Realizámos dois subestudos independentes (Estudos 2 e 3) em diferentes centros de investigação, com diferentes equipamentos e diferentes equipas de investigação para validação. Utilizámos um alinhamento liberal dos ambientes de investigação, permitindo uma razoável heterogeneidade nos procedimentos, equipamentos, sequências de imagem, linguagem de comunicação entre participantes-investigador e centros de estudo, introduzindo uma razoável heterogeneidade no procedimento de validação para garantir a generalização.

Certifique-se de que a previsão é guiada pelo sinal neural e é específica para a sensibilidade à dor

Para garantir que o marcador proposto de sensibilidade à dor seja de fato guiado por sinais neurais associados à sensibilidade à dor, avaliamos a correlação da pontuação prevista com várias variáveis de confundidor e validador pré-definidas (e pré-registadas).

Asseguir a acessibilidade dos resultados

Aplicamos um pré-registo abrangente e disponibilizamos o código fonte do método de código aberto e gratuito para a comunidade. Além disso, nós fornecemos um container docker independente de plataforma e fácil de usar, que oferece a oportunidade de usar nosso modelo preditivo como um produto de pesquisa50, para obter previsões de sensibilidade à dor fora da caixa de quaisquer conjuntos de dados de imagem apropriados.

Participantes

Um total de N = 116 voluntários jovens e saudáveis estavam envolvidos em três subestudos. A idade e sexo dos participantes é relatada na Tabela Complementar 1. Estudo 1 envolveu N1 = 39 participantes (a mesma amostra que na ref. 8). Foi realizado na Universidade Ruhr de Bochum (Alemanha) pela MZ e TSW e utilizado como amostra de formação para a previsão da sensibilidade à dor baseada na aprendizagem mecânica e, adicionalmente, serviu como base para a validação interna da previsão. Os estudos 2 e 3 (N2 = 48, N3 = 29) foram realizados no Hospital Universitário de Essen (Alemanha) pela FS e TS e na Universidade de Szeged (Hungria) pela BK e TK, respectivamente, e serviram como amostras para validação externa. Os critérios de inclusão e os critérios de exclusão foram em grande parte idênticos nos três centros e estão listados na Tabela 3. As políticas de recrutamento e reembolso variaram entre os centros; os participantes receberam 20 euros/h nos Estudos 1 e 2 e nenhum reembolso no Estudo 3.

Insplante metálico, piercing inamovível, pacificador, tatuagem na posição cabeça/ pescoço, gravidez ou claustrofobia conhecida foram considerados como contra-indicação para a medição da RM. Os participantes foram obrigados a se abster de consumir cafeína duas horas antes dos experimentos (exceto no Estudo 3) e de consumir álcool no dia do teste e no dia anterior.

O estudo foi conduzido de acordo com a Declaração de Helsinque, obedece a todos os regulamentos éticos relevantes para o trabalho com participantes humanos e foi aprovado pelos comitês de ética locais ou nacionais (Números de Registro: 4974-14, 18-8020-BO e 057617/2015/OTIG na Universidade Ruhr de Bochum, Hospital Universitário de Essen e ETT TUKEB Hungria, respectivamente). Todos os participantes deram consentimento livre e esclarecido por escrito antes dos testes.

Testes sensoriais por imagem e quantitativos (QST) foram realizados no mesmo dia no Estudo 1 e em média 2-3 dias de intervalo nos Estudos 2 e 3 (ver Tabela Complementar 1 para detalhes). As medidas de RM sempre precederam a sessão de QST.

Medidas-funcionais de RM

Medições de fMRI de alta resolução anatômica e de estado de repouso com os olhos abertos foram adquiridas de todos os participantes. Os parâmetros de varredura (incluindo equipamentos) variaram entre os centros e estão listados na Tabela 4. Durante as medições, os participantes foram instruídos a ficar quietos e relaxados, sem adormecer, e a evitar qualquer movimento. O almofadado de espuma, e nos Estudos 1 e 2, foram utilizadas almofadas pneumáticas para restringir os movimentos da cabeça. Todas as medidas anatômicas de RM foram triadas para achados incidentais.

Medidas-QST

Calor (HPT), frio (CPT) e limiares de dor mecânicos (MPT) foram adquiridos de acordo com o protocolo QST28. Os limiares de calor (WDT), frio (CDT) e nos Estudos 2 e 3, mecânicos (MDT) de detecção foram obtidos como medidas de controle adicionais. Todas as medidas sensoriais foram obtidas a partir do antebraço esquerdo palmar, proximal à crista do punho. Dentro da estrutura QST, os limiares térmicos são determinados através de um método de limites. Para este fim, temperaturas crescentes e decrescentes foram aplicadas à pele com um estimulador térmico MSA (Somedic, Hörby, Suécia) nos Estudos 1 e Estimuladores térmicos Pathway (Medoc Ltd., Ramat Yishai, Israel) nos Estudos 2 e 3. Em todos os estudos, foram utilizados termodes ATS em uma superfície da pele de 30 × 30mm, com uma temperatura de base de 32 °C. Os participantes foram instruídos a indicar o início da dor por pressão de botão. Para todos os limiares térmicos 6, em vez de 3 (como no protocolo original)28, foram realizadas repetições de estímulos para reduzir a variância entre sujeitos. Além disso, a primeira medida foi descartada da análise como um estímulo de teste. HPT e CPT foram calculados como a média aritmética das cinco temperaturas limiares restantes. Os MPTs e os MDTs foram determinados usando um método de escada. Cinco estímulos crescentes e cinco decrescentes de pinprick (MRC Systems, Heidelberg, Alemanha) foram aplicados ao antebraço esquerdo palmar de forma alternada, enquanto o participante foi instruído a classificar os estímulos como nocivos, ou não nocivos. O limiar de detecção mecânica foi avaliado de forma análoga com os estímulos do filamento de von Frey. MPT e MDT foram computados como a força média geométrica log-transformada determinada em cinco degraus de escada ascendentes e descendentes.

Medidas adicionais

Idade, sexo, altura, peso e, para as participantes do sexo feminino, a data do primeiro dia da última menstruação e o uso de contraceptivos, foi registrada antes de todas as medidas. Além disso, o auto-relato semanal de consumo de álcool e nível de escolaridade (escola primária, secundária, universitária) foi registrador para os Estudos 1 e 2. Antes do QST, os participantes preencheram o Questionário de Sensibilidade à Dor (PSQ)56, a Escala de Catastrofização da Dor (PCS)57, o Inventário de Ansiedade do Estado-Traço (STAI)58, a versão curta em alemão da Escala de Depressão (ADS-K, Center for Epidemiologic Studies)59 e, adicionalmente nos Estudos 2 e 3, o Índice de Qualidade do Sono de Pittsburgh (PSQI)60 e o questionário de estresse percebido (PSQ20)61. Nos Estudos 2 e 3, a pressão arterial foi medida tanto antes da ressonância magnética como das medidas de QST. Além disso, para a Amostra 1, os valores de T50 estavam disponíveis em um experimento paralelo realizado no dia anterior ao teste de fMRI. T50 representa a temperatura (em °C) necessária para induzir uma classificação de calor-dor de 50 (em uma escala que varia de 0, sem dor a 100 dor insuportável). Os valores de T50 foram obtidos a partir de uma interpolação não linear (polinômio de segunda ordem) das classificações obtidas em resposta a 15 estímulos tônicos de calor-dores (duração: 16 s) entre 42,5 °C e 48 °C, apresentados de forma pseudo-randomizada de pesquisa em grade.

Cálculo da sensibilidade à dor

A variável alvo para a previsão foi uma medida composta única de sensibilidade à dor individual, resumindo HPT, CPT e MPT como definido na ref. 8,

No Estudo 1, HPT, CPT e MPT foram transformados em Z (média centrada e padronizada) e HPT, bem como MPT foram invertidos (multiplicados por -1), de modo que valores Z mais elevados denotaram uma maior sensibilidade à dor. Em seguida, a média aritmética das variáveis transformadas em Z foi computada para cada participante e definida como pain-sensitivity score. Nos Estudos 2 e 3, o mesmo procedimento foi aplicado, exceto que a transformação em Z foi baseada no meio populacional e no desvio padrão do Estudo 1, para garantir que a mesma escala fosse utilizada em todos os estudos. Os valores extremos de QST foram definidos usando os percentis normativos de 95% relatados na ref. 28; os participantes que apresentaram valores extremos de HPT, CPT ou MPT em pelo menos duas das três modalidades foram excluídos. Esta triagem resultou na exclusão de 0, 3 e 2 participantes nas amostras 1, 2 e 3, respectivamente (Tabela Complementar 2).

fMRI pré-processamento

As conexões funcionais baseadas em fMRI são suscetíveis a artefatos de movimento no scanner62,63, o pré-processamento apropriado e a limpeza do sinal são fundamentais para o sucesso da previsão baseada na conectividade. Os dados de ressonância magnética funcional foram pré-processados de forma idêntica nos três estudos. O fluxo de trabalho aplicado, baseado em nipype, está descrito na Fig. 1. Ele utilizou software de neuroimagem de terceiros, código adaptado das ferramentas de software C-PAC53 e niworkflows54, e rotinas internas de pitão.

Retirada do cérebro tanto da imagem anatômica quanto da estrutural, assim como a fragmentação dos tecidos das imagens anatômicas foi realizada com FSL bet e fast64. As imagens anatômicas foram registradas linearmente e não linearmente de forma co-registrada com a técnica de registro baseado no limite de 1mm de resolução MNI152 padrão do cérebro com ANTs65 (ver https://gist.github.com/spisakt/0caa7ec4bc18d3ed736d3a4e49da7415 para o código fonte).

As imagens funcionais foram co-registradas para as imagens anatômicas com a técnica de registro baseado no limite de FSL flirt. Todas as transformações resultantes foram salvas para uso posterior. O pré-processamento de imagens funcionais aconteceu no espaço de imagens nativas, sem reamostragem. A correção de movimento baseada em realigmentação foi realizada com o mcflirt FSL. As seis estimativas de movimento da cabeça resultantes (3 rotações, 3 traduções), suas versões quadradas, seus derivados e os derivados quadrados (conhecidos como a expansão Friston-2466) foram calculados e salvos para correção de incômodos. Além disso, o movimento da cabeça foi resumido como séries temporais de deslocamento de quadro (FD), de acordo com o método de Power63, para ser usado na censura e exclusão de dados. Após a correção do movimento, os outliers (por exemplo, picos de movimento) nos dados das séries temporais foram atenuados usando o despike AFNI67. A união dos mapas de matéria branca erodidos e das máscaras ventriculares foram transformadas para o espaço funcional nativo e utilizadas para extrair o sinal de ruído para a correção anatômica do CompCor68,

Em um passo de regressão de incômodo, 6 parâmetros CompCor (os 6 primeiros componentes principais das séries temporais de ruído-região), os parâmetros de movimento Friston-24 e a tendência linear foram removidos dos dados das séries temporais com um modelo linear geral. Nos dados residuais, a filtragem do passo de banda temporal foi realizada com a banda 3DBandpass da AFNI para reter a banda de frequência de 0,008-0,08 Hz. Espera-se que o uso prévio do despique do AFNI atenue a similaridade dos artefatos de movimento residual nos quadros temporais vizinhos durante a filtragem do passo de banda69. Para atenuar ainda mais o impacto dos artefactos de movimento, foram descartados dos dados (conhecidos como “scrubbing the data”)70 os intervalos de tempo potencialmente contaminados pelo movimento, definidos por um limite conservador FD > 0,15 mm. Os participantes foram excluídos de análises posteriores se a média de FD excedesse 0,15 mm, ou quando mais de 30% dos quadros foram esfregados. Isto resultou na exclusão de 4, 8 e 7 participantes nas amostras 1, 2 e 3, respectivamente (Tabela Complementar 2). O controle de qualidade (checagem de cadastro, carpet-plots, ver, por exemplo, Figs Suplementares 2-4) foi realizado durante todo o fluxo de trabalho.

Análise de conectividade funcional

A versão de 122-parcelas do atlas funcional multi-resolução MIST71 e as máscaras de material cinzento obtidas da imagem anatômica foram transformadas para o espaço funcional nativo. Este atlas (construído com o método BASC, ou seja, análise bootstrap de clusters estáveis) teve recentemente um bom desempenho na modelagem preditiva baseada em conectividade72. As regiões do atlas do espaço nativo foram mascaradas com as máscaras de matéria cinzenta que foram obtidas a partir da imagem anatômica e transformadas anteriormente para o espaço funcional. Com esta técnica de atlas-individualização, o sinal regional final terá origem – com uma elevada probabilidade – em voxels cinza-matéria para cada sujeito (que verificamos cuidadosamente manualmente para todos os sujeitos), enquanto que com o método convencional, uma proporção variável de voxels cinza e branco-matéria é incluída para cada sujeito. Portanto, espera-se que a entrada de informações do processo de segmentação de tecidos diminua a variabilidade de assunto para assunto (veja a Fig. 5 Suplementar para exemplos). As séries temporais Voxel foram calculadas em média sobre estas regiões MIST individualizadas e, juntamente com o sinal cinza-matéria média, retidas para análise de conectividade baseada em gráficos.

As séries temporais regionais foram ordenadas em módulos funcionais de grande escala (definidos pelo atlas MIST de 7 partes) para fins de visualização (Fig. 1). A correlação parcial foi computada em todos os pares de regiões (e matéria cinzenta global), conforme implementado no módulo python nilearn73. Correlações parciais, ao invés de simples correlações, foram usadas para descartar a conectividade indireta74. Nossa abordagem de modelagem gráfica assegurou que o sinal global de matéria cinzenta seja tratado como um confuso durante o cálculo dos coeficientes de correlação parcial mas, ao mesmo tempo, também o considerou como um sinal de interesse, uma vez que pode representar processos relacionados à vigilância75. Os coeficientes de correlação parciais foram organizados em 123 por 123 (122 regiões + sinal global de matéria cinzenta) matrizes de conectividade simétrica. O triângulo superior dessas matrizes foi usado como espaço de características para modelagem preditiva baseada em aprendizagem mecânica.

Treinamento e validação de modelos preditivos

Dados de conectividade funcional de estado de repouso de cérebro-inteiro do estudo 1 (N1 = 35, após todas as exclusões, como na ref. 8, Tabela Complementar 2) foi usado como espaço de características de entrada (P = 7503 características por participante) para prever pontuações individuais de sensibilidade à dor, levando a uma grande configuração P-pequeno N.

Construímos um pipeline de aprendizagem da máquina (https://github.com/spisakt/RPN-signature/blob/master/PAINTeR/model.py) em scikit-learn76, que consiste em escalas robustas de características (remove a mediana e escalas com quantil de dados), pré-seleção de características77, seleção das melhores características K com relações mais fortes com a variável alvo e um modelo de regressão Elastic Net78 (um modelo linear com combinações de L1 e L2-normas como regularizador). O uso da rede elástica foi uma decisão tomada antes da análise. Nossa principal motivação para escolher a rede elástica foi que ela permite otimizar a esparsidade (regularização L1 vs. L2) como um hiperparâmetro, de forma que não tivemos que fazer nenhuma suposição a priori sobre a esparsidade da verdade do terreno discriminativo (ver a ref. 79 para a fundamentação). Em resumo, os hiperparâmetros livres do gasoduto de aprendizagem da máquina foram o número de características pré-seleccionadas (K), o rácio de L1/L2-regularização e o peso (alfa) da regularização. Os hiperparâmetros foram otimizados com um procedimento de pesquisa em grade e erro médio quadrático negativo como função de custo. Os valores para K variaram de 10 a 200 com incrementos de 5, e incluídos para a relação L1/L2 para alfa. A otimização dos hiperparâmetros foi realizada em uma validação cruzada de leave-one-participant-out (fase de validação interna). A validação cruzada incorporou o pipeline completo de aprendizagem da máquina para evitar a introdução de dependências entre o treinamento e as amostras de teste. Observe que o pré-processamento fMRI foi independente entre os sujeitos, não sendo, portanto, incluído na validação cruzada. Os hiperparâmetros ideais foram K = 25, L1/L2-ratio = 0,999 e alfa = 0,005,

Avalidação externa foi realizada aplicando a assinatura RPN nos dados fMRI dos Estudos 2 e 3 (N2 = 37, N3 = 19, após exclusões, Tabela Complementar 2), simplesmente aplicando a transformação de característica (escala) obtida na Amostra 1 e, em seguida, calculando o produto de pontos entre matrizes de conectividade individual e os pesos de característica não zero obtidos na Amostra 1. As previsões resultantes foram comparadas com os escores de sensibilidade à dor baseados em QST observados através do cálculo do erro médio absoluto (MAE), erro médio quadrático (MSE) e da variância explicada. Os valores de p baseados em permissão foram obtidos para as três medidas, usando a embalagem de pitão mlxtend. Além disso, o bootstrapping com cobertura condicional80 foi usado para fornecer p-values para pesos de conectividade preditivos para ajudar na interpretação. Construímos 10000 amostras de bootstrap (com reposição), com um tamanho igual ao da amostra original, consistindo de dados pareados do cérebro e dos resultados. O modelo preditivo com os hiperparâmetros ideais foi adaptado a cada amostra. Os valores de P não corrigidos foram calculados para cada conexão selecionada com base na proporção de pesos abaixo ou acima de zero, como por exemplo na ref. 30. Note que a interpretação destes valores de p e intervalos de confiança (Tabela Suplementar 4) permanece limitada, uma vez que estão condicionados ao procedimento de selecção de características.

Análise do Fundador

Para explorar potenciais variáveis de confusão, as predições de sensibilidade à dor (ou predições cruzadas no caso da Amostra 1) foram contrastadas com a média e a mediana da DF, a porcentagem de volumes esfregados, a pressão arterial sistólica e diastólica antes da RM e da medida da QST (como a pressão arterial foi anteriormente relatada81 para ser associada à sensibilidade à dor mecânica), o tempo decorrido entre a RM e o teste de QST (para testar a estabilidade temporal da previsão), idade, sexo, IMC, número de dias desde o primeiro dia da última menstruação, consumo de álcool (unidades/semana), nível de educação, estado e ansiedade traumática (STAI), escore de sintomas depressivos (ADS-K), sensibilidade à dor auto-relatada (PSQ) e dor catastrófica (PCS), estresse percebido (PSQ20), qualidade do sono (PSQI), e limiares de detecção de QST não nocivos (CDT, WDT e MDT, quando disponíveis). Além disso, no Estudo 1 as previsões foram comparadas aos valores T50 e aos níveis de GABA e Glutamato/Glutamina baseados em espectroscopia de RM nas regiões cerebrais que processam a dor (ver ref. 8 para detalhes). As associações foram testadas com modelos lineares baseados em permutação.

Visualização da rede preditiva

As conexões preditivas inter-regionais destacadas pelos coeficientes de regressão não zero da assinatura RPN foram exibidas como um gráfico de fitas usando o circlize do pacote R (Fig. 3). As máscaras de região cerebral correspondentes individualizadas foram transformadas de volta ao espaço padrão para criar um mapa de probabilidade regional específico do estudo (refletindo a exatidão do co-registro e a variabilidade individual na morfologia). Os mapas de probabilidade foram multiplicados pela soma dos coeficientes de regressão correspondentes para criar um mapa de força preditiva regional, que foi então visualizado com FSLeyes e MRIcroGL.

(https://www.mccauslandcenter.sc.edu/mricrogl) (Fig. 3). A análise do envolvimento da rede de estado de repouso em larga escala (como definido pelo atlas cerebral MIST71) foi realizada resumindo e transformando em Z os valores de voxel através das sete regiões de interesse. O gráfico Polar foi feito com o pacote R ggplot2.

Disponibilidade do software

As pontuações da assinatura RPN podem ser computadas com base nos conjuntos de dados funcionais estruturais e de estado de repouso pela ferramenta de software com o mesmo nome. A ferramenta de software RPN-signature consiste no pipeline de processamento de RM descrito e no modelo preditivo funcional baseado em connectome-based. Ela está disponível como código fonte em https://github.com/spisakt/RPN-signature. Como o software segue a Estrutura de Dados de Brain Imaging (BIDS)82 e a especificação BIDS-App, ele fornece uma interface de linha de comando padrão e se baseia na tecnologia Docker. A imagem da doca é depositada no Docker Hub: (https://cloud.docker.com/repository/docker/tspisak/rpn-signature) e não depende de nenhum software fora da imagem do recipiente. Isto, juntamente com o desenvolvimento baseado em integração contínua totalmente transparente e a etiquetagem automática e versionamento, aumenta a disponibilidade do software e suporta a reprodutibilidade dos resultados da assinatura RPN.

Resumo do relatório

Outras informações sobre o design da pesquisa estão disponíveis no Nature Research Reporting Summary ligado a este artigo.