La connettività cerebrale funzionale a riposo senza dolore predice la sensibilità individuale al dolore

Considerazioni generali

Il disegno dello studio è stato stabilito con un’attenta considerazione delle recenti raccomandazioni, requisiti e standard per i biomarcatori di neuroimmagini50 (neuromarcatori) e motivato dai seguenti pensieri.

Massimizzare le prestazioni predittive

Abbiamo impiegato una pipeline di pre-elaborazione standardizzata per garantire una sensibilità ottimale del neuromarcatore, poiché una dimensione dell’effetto sufficiente è un requisito fondamentale di qualsiasi utilità clinica50. Abbiamo usato l’allineamento dell’immagine ad alta precisione, incorporando l’anatomia individuale durante l’estrazione dei dati fMRI timeseries. Inoltre, abbiamo adottato recenti raccomandazioni e protocolli51 per quanto riguarda la riduzione degli artefatti e ottimizzato il nostro flusso di lavoro per soddisfare le particolari esigenze di analisi connectome-based. Abbiamo usato il nostro in-house sviluppato, open-source libreria software python Pipelines Utilising a Modular Inventory (PUMI, https://github.com/spisakt/PUMI), che si basa su nipype52, un progetto Python community-based che fornisce un’interfaccia uniforme al software esistente neuroimaging e, in parte, riutilizzato il codice dal C-PAC53 e il niworkflows54 progetti open-source. È stato utilizzato un approccio di modellazione predittiva (apprendimento automatico) per sfruttare i ricchi dati forniti dalle reti cerebrali funzionali a riposo e, potenzialmente, trarre vantaggio dall’iperacuità fMRI55.

Valutazione del potere predittivo in condizioni realistiche

Abbiamo utilizzato una strategia di validazione esterna pre-registrata, che separava rigorosamente l’addestramento del modello dalla valutazione delle prestazioni. Per la formazione del modello, abbiamo utilizzato esclusivamente i dati dello studio 1. Abbiamo condotto due sottostudi indipendenti (Studi 2 e 3) in diversi centri di ricerca, con diverse attrezzature e diverso personale di ricerca per la convalida. Abbiamo usato un allineamento liberale delle impostazioni di ricerca, consentendo una ragionevole eterogeneità nelle procedure, nelle attrezzature, nelle sequenze di imaging, nella lingua di comunicazione partecipante-ricercatore tra i centri di studio, introducendo una ragionevole eterogeneità nella procedura di convalida per garantire la generalizzabilità.

Assicurare che la predizione sia guidata dal segnale neurale e sia specifica per la sensibilità al dolore

Per assicurare che il marker proposto per la sensibilità al dolore sia effettivamente guidato da segnali neurali associati alla sensibilità al dolore, abbiamo valutato la correlazione del punteggio previsto con vari confonditori e variabili di validazione predefiniti (e preregistrati).

Assicurare l’accessibilità dei risultati

Abbiamo applicato una pre-registrazione completa e reso il codice sorgente del metodo open-source e liberamente disponibile per la comunità. Inoltre, forniamo un contenitore docker indipendente dalla piattaforma e facile da usare, che offre l’opportunità di utilizzare il nostro modello predittivo come prodotto di ricerca50, per ottenere previsioni di sensibilità al dolore out-of-the box da qualsiasi set di dati di imaging appropriato.

Partecipanti

Un totale di N = 116 volontari giovani e sani sono stati coinvolti in tre sotto-studi. L’età e il sesso dei partecipanti sono riportati nella tabella supplementare 1. Lo studio 1 ha coinvolto N1 = 39 partecipanti (lo stesso campione del rif. 8). È stato eseguito presso la Ruhr University Bochum (Germania) da MZ e TSW e utilizzato come campione di formazione per la previsione basata sull’apprendimento automatico della sensibilità al dolore e, inoltre, è servito come base per la validazione interna della previsione. Gli studi 2 e 3 (N2 = 48, N3 = 29) sono stati eseguiti presso l’ospedale universitario di Essen (Germania) da FS e TS e presso l’Università di Szeged (Ungheria) da BK e TK, rispettivamente, e sono serviti come campioni per la validazione esterna. I criteri di inclusione e di esclusione erano in gran parte identici in tutti e tre i centri e sono elencati nella tabella 3. Le politiche di reclutamento e di rimborso variavano tra i centri; i partecipanti hanno ricevuto 20 €/h negli studi 1 e 2 e nessun rimborso nello studio 3.

Impianto metallico, piercing inamovibile, peacemaker, tatuaggio in posizione testa/collo, gravidanza o claustrofobia nota sono stati considerati come controindicazione alla misurazione della risonanza magnetica. Ai partecipanti è stato richiesto di astenersi dal consumare caffeina due ore prima degli esperimenti (tranne nello studio 3) e dal consumare alcol il giorno del test e il giorno precedente.

Lo studio è stato condotto in conformità con la Dichiarazione di Helsinki, è conforme a tutte le norme etiche pertinenti per il lavoro con i partecipanti umani ed è stato approvato dai comitati etici locali o nazionali (numeri di registro: 4974-14, 18-8020-BO e 057617/2015/OTIG rispettivamente presso la Ruhr University Bochum, l’University Hospital Essen e l’ETT TUKEB Ungheria). Tutti i partecipanti hanno dato il consenso informato scritto prima del test.

Imaging e test sensoriale quantitativo (QST) sono stati eseguiti lo stesso giorno nello studio 1 e in media 2-3 giorni di distanza negli studi 2 e 3 (vedi tabella 1 supplementare per i dettagli). La misurazione della risonanza magnetica ha sempre preceduto la sessione di QST.

Misure di risonanza magnetica funzionale

Sono state acquisite misure di fMRI anatomica ad alta risoluzione e in stato di riposo ad occhi aperti da tutti i partecipanti. I parametri di scansione (comprese le attrezzature) variavano tra i centri e sono elencati nella tabella 4. Durante le misurazioni, i partecipanti sono stati istruiti a stare fermi e rilassati, senza addormentarsi, ed evitare qualsiasi movimento. Per limitare i movimenti della testa sono state utilizzate imbottiture di schiuma e, negli studi 1 e 2, cuscini pneumatici. Tutte le misurazioni anatomiche MRI sono state esaminate per i risultati accidentali.

Misure-QST

Le soglie del dolore caldo (HPT), freddo (CPT) e meccanico (MPT) sono state acquisite secondo il protocollo QST28. Calore (WDT), freddo (CDT) e nello studio 2 e nello studio 3, le soglie di rilevamento meccanico (MDT) sono state ottenute come misure di controllo aggiuntive. Tutte le misurazioni sensoriali sono state ottenute dall’avambraccio sinistro palmare, prossimale alla cresta del polso. Nell’ambito del QST, le soglie termiche sono determinate utilizzando un metodo di limiti. A tal fine, temperature crescenti e decrescenti sono state applicate alla pelle con uno stimolatore termico MSA (Somedic, Hörby, Svezia) nello studio 1 e con gli stimolatori termici Pathway (Medoc Ltd., Ramat Yishai, Israele) negli studi 2 e 3. In tutti gli studi, i termodi ATS sono stati utilizzati su una superficie cutanea di 30 × 30 mm, con una temperatura di base di 32 °C. I partecipanti sono stati istruiti a indicare l’inizio del dolore premendo un pulsante. Per tutte le soglie termiche 6, piuttosto che 3 (come nel protocollo originale)28, sono state eseguite ripetizioni dello stimolo per ridurre la varianza tra soggetti. Inoltre, la prima misurazione è stata scartata dall’analisi come stimolo di prova. HPT e CPT sono stati calcolati come le medie aritmetiche delle cinque temperature di soglia rimanenti. MPT e MDT sono stati determinati utilizzando un metodo a scale. Cinque treni crescenti e cinque decrescenti di stimoli pinprick (MRC Systems, Heidelberg, Germania) sono stati applicati all’avambraccio palmare sinistro in modo alternato, mentre il partecipante è stato istruito a classificare gli stimoli come nocivi o non nocivi. La soglia di rilevamento meccanico è stata valutata in modo analogo con le stimolazioni dei filamenti di von Frey. MPT e MDT sono stati calcolati come la forza media geometrica log-trasformata determinata in cinque scale ascendenti e discendenti di soglia.

Misure aggiuntive

Età, sesso, altezza auto dichiarata, peso e, per le partecipanti donne, la data del primo giorno delle ultime mestruazioni e l’uso di contraccettivi, è stato registrato prima di tutte le misurazioni. Inoltre, il consumo settimanale di alcol e il livello di istruzione (scuola primaria, scuola secondaria, università) sono stati registrati per gli studi 1 e 2. Prima del QST, i partecipanti hanno compilato il Pain Sensitivity Questionnaire (PSQ)56, la Pain Catastrophizing Scale (PCS)57, lo State-Trait Anxiety Inventory (STAI)58, la versione breve tedesca della Depression Scale (ADS-K, Center for Epidemiologic Studies)59 e, in aggiunta negli Studi 2 e 3, il Pittsburgh Sleep Quality Index (PSQI)60 e il questionario dello stress percepito (PSQ20)61. Negli studi 2 e 3, la pressione sanguigna è stata misurata sia prima della risonanza magnetica che delle misurazioni QST. Inoltre, per il campione 1, i valori T50 erano disponibili da un esperimento parallelo eseguito il giorno prima del test fMRI. T50 rappresenta la temperatura (in °C) necessaria per indurre una valutazione di calore-dolore di 50 (su una scala che va da 0, nessun dolore a 100 dolore insopportabile). I valori T50 sono stati ottenuti da un’interpolazione non lineare (polinomiale del secondo ordine) delle valutazioni ottenute in risposta a 15 stimoli tonici di calore-dolore (durata: 16 s) tra 42,5 °C e 48 °C, presentati in modo pseudo-randomizzato di ricerca a griglia.

Calcolo della sensibilità al dolore

La variabile obiettivo per la previsione era una singola misura composita della sensibilità al dolore individuale che riassume HPT, CPT e MPT come definito in rif. 8.

Nello studio 1, HPT, CPT e MPT sono stati trasformati a Z (media centrata e standardizzata) e HPT, così come MPT sono stati invertiti (moltiplicati per -1), in modo che i valori Z più alti denotassero una maggiore sensibilità al dolore. Poi, la media aritmetica delle variabili trasformate in Z è stata calcolata per ogni partecipante e definita come punteggio di sensibilità al dolore. Negli studi 2 e 3, è stata applicata la stessa procedura, tranne che la trasformazione Z è stata basata sulla media della popolazione e la deviazione standard dello studio 1, per garantire che la stessa scala fosse usata in tutti gli studi. I valori estremi di QST sono stati definiti utilizzando i percentili normativi del 95% riportati in rif. 28; i partecipanti che mostravano valori estremi di HPT, CPT o MPT in almeno due delle tre modalità sono stati esclusi. Questo screening ha portato ad escludere 0, 3 e 2 partecipanti nei campioni 1, 2 e 3, rispettivamente (Tabella supplementare 2).

preprocessing fMRI

Come la connettività funzionale basata su fMRI è suscettibile di artefatti di movimento in-scanner62,63, preprocessing appropriato e pulizia del segnale è la chiave per la previsione di successo basato sulla connettività. I dati di risonanza magnetica funzionale a riposo sono stati preprocessati in modo identico in tutti e tre gli studi. Il flusso di lavoro applicato, nipype-based è raffigurato nella Fig. 1 supplementare. Ha utilizzato software di neuroimaging di terze parti, il codice adattato dagli strumenti software C-PAC53 e niworkflows54, e in-house routine python.

Estrazione del cervello da entrambe le immagini anatomiche e strutturali, così come, tessuto-segmentazione dalle immagini anatomiche è stata eseguita con FSL bet e veloce64. Immagini anatomiche sono stati linearmente e non linearmente co-registrati alla risoluzione 1mm standard MNI152 cervello modello cervello con ANTs65 (vedi https://gist.github.com/spisakt/0caa7ec4bc18d3ed736d3a4e49da7415 per il codice sorgente).

Immagini funzionali sono stati co-registrati alle immagini anatomiche con la tecnica di registrazione boundary-based di FSL flirt. Tutte le trasformazioni risultanti sono state salvate per un uso successivo. La pre-elaborazione delle immagini funzionali è avvenuta nello spazio nativo dell’immagine, senza ricampionamento. La correzione del movimento basata sul riallineamento è stata eseguita con FSL mcflirt. Le sei stime di movimento della testa risultanti (3 rotazioni, 3 traslazioni), le loro versioni al quadrato, le loro derivate e le derivate al quadrato (note come Friston-24-expansion66) sono state calcolate e salvate per la correzione dei disturbi. Inoltre, il movimento della testa è stato riassunto come serie temporali di spostamento framewise (FD), secondo il metodo di Power63, da utilizzare nella censura e nell’esclusione dei dati. Dopo la correzione del movimento, outlier (ad esempio picchi di movimento) nei dati di serie temporali sono stati attenuati utilizzando AFNI despike67. L’unione delle mappe erose bianco-materia e maschere ventricolo sono stati trasformati nello spazio funzionale nativo e utilizzato per l’estrazione di rumore-segnale per la correzione CompCor anatomico68.

In un passo di regressione fastidioso, 6 parametri CompCor (le 6 prime componenti principali della timeserie di rumore-regione), i parametri di movimento Friston-24 e la tendenza lineare sono stati rimossi dai dati timeserie con un modello lineare generale. Sui dati residui, il filtraggio passa-banda temporale è stato eseguito con 3DBandpass di AFNI per mantenere la banda di frequenza 0.008-0.08 Hz. L’uso precedente di AFNI despike è previsto per attenuare aliasing di artefatti di movimento residuo nel vicino time-frames durante il filtraggio bandpass69. Per attenuare ulteriormente l’impatto degli artefatti di movimento, potenzialmente motion-contaminated time-frames, definito da un conservatore FD > 0,15 mm soglia, sono stati eliminati dai dati (noto come scrubbing i dati)70. I partecipanti sono stati esclusi da ulteriori analisi se l’FD medio superava 0,15 mm, o quando più del 30% dei fotogrammi veniva eliminato. Questo ha portato all’esclusione di 4, 8 e 7 partecipanti nei campioni 1, 2 e 3, rispettivamente (tabella 2 supplementare). Il controllo di qualità (controllo di registrazione, carpet-plot, vedi ad esempio le figure supplementari 2-4) è stato eseguito durante tutto il flusso di lavoro.

Analisi della connettività funzionale

La versione 122-parcel dell’atlante cerebrale funzionale multirisoluzione MIST71 e le maschere di materia grigia ottenute dall’immagine anatomica sono state trasformate nello spazio funzionale nativo. Questo atlante (costruito con il metodo BASC, cioè l’analisi bootstrap di cluster stabili) è stato recentemente dimostrato di eseguire bene nella modellazione predittiva basata sulla connettività72. Native-spazio regioni atlante sono stati mascherati con le maschere di materia grigia che sono stati ottenuti dall’immagine anatomica e trasformato in spazio funzionale in precedenza. Con questa tecnica di atlante-individualizzazione, il segnale regionale finale avrà origine – con un’alta probabilità – da voxel di materia grigia per ogni soggetto (che abbiamo accuratamente controllato manualmente per tutti i soggetti), mentre con il metodo convenzionale, un rapporto variabile di grigio e bianco-materia voxel sono inclusi per ogni soggetto. Pertanto, l’inserimento di informazioni dal processo di segmentazione dei tessuti dovrebbe diminuire la variabilità da soggetto a soggetto (vedi Fig. 5 supplementare per esempi). Le serie temporali dei voxel sono state mediate su queste regioni MIST individualizzate e, insieme al segnale medio di materia grigia, conservate per l’analisi di connettività basata sul grafico.

Le serie temporali regionali sono state ordinate in moduli funzionali su larga scala (definiti dall’atlante MIST a 7 parcelle) per scopi di visualizzazione (Fig. 1). La correlazione parziale è stata calcolata attraverso tutte le coppie di regioni (e materia grigia globale), come implementato nel modulo python nilearn73. Parziale, piuttosto che semplici correlazioni sono stati utilizzati per escludere la connettività indiretta74. Il nostro approccio graph-modelling assicurato che il segnale di materia grigia globale è gestito come un confondimento durante il calcolo dei coefficienti di correlazione parziale, ma, allo stesso tempo, anche considerato come un segnale di interesse, in quanto può rappresentare processi correlati vigilanza75. I coefficienti di correlazione parziale sono stati organizzati in matrici di connettività simmetriche 123 per 123 (122 regioni + segnale globale di materia grigia). Il triangolo superiore di queste matrici è stato utilizzato come spazio delle caratteristiche per la modellazione predittiva basata sull’apprendimento automatico.

Formazione e convalida del modello predittivo

I dati di connettività funzionale a riposo dell’intero cervello dello studio 1 (N1 = 35, dopo tutte le esclusioni, come in rif. 8, Tabella supplementare 2) sono stati utilizzati come spazio delle caratteristiche di input (P = 7503 caratteristiche per partecipante) per prevedere i punteggi individuali di sensibilità al dolore, portando a una grande impostazione P-piccola N.

Abbiamo costruito una pipeline di machine-learning (https://github.com/spisakt/RPN-signature/blob/master/PAINTeR/model.py) in scikit-learn76, che consiste in un robusto feature scaling (rimuove la mediana e scala con i quantili dei dati), nella preselezione delle caratteristiche77, nella selezione delle K migliori caratteristiche con le relazioni più forti con la variabile target e in un modello di regressione Elastic Net78 (un modello lineare con L1 combinato e L2-norms come regolarizzatore). L’uso della rete elastica è stata una decisione presa prima dell’analisi. La nostra motivazione principale per scegliere Elastic Net è stata che permette di ottimizzare la sparsità (regolarizzazione L1 vs. L2) come iperparametro, in modo da non dover fare alcuna ipotesi a-prioritaria sulla sparsità della ground truth discriminativa (vedi rif. 79 per il razionale). Per riassumere, gli iperparametri liberi della pipeline di apprendimento automatico erano il numero di caratteristiche preselezionate (K), il rapporto della regolarizzazione L1/L2 e il peso (alfa) della regolarizzazione. Gli iperparametri sono stati ottimizzati con una procedura di ricerca a griglia e un errore quadratico medio negativo come funzione di costo. I valori per K andavano da 10 a 200 con incrementi di 5, e inclusi per il rapporto L1/L2 per alfa. L’ottimizzazione degli iperparametri è stata eseguita in una convalida incrociata leave-one-participant-out (fase di convalida interna). La convalida incrociata ha incorporato l’intera pipeline di apprendimento automatico per evitare di introdurre dipendenze tra i campioni di formazione e di prova. Si noti che la pre-elaborazione fMRI era indipendente tra i soggetti, quindi non inclusa nella convalida incrociata. Iperparametri ottimali sono stati trovati per essere K = 25, L1/L2-ratio = 0.999 e alfa = 0.005.

Convalida esterna è stata eseguita applicando la RPN-signature sui dati fMRI degli studi 2 e 3 (N2 = 37, N3 = 19, dopo le esclusioni, Tabella supplementare 2), semplicemente applicando la trasformazione caratteristica (scaling) ottenuto sul campione 1 e quindi calcolando il prodotto di punto tra le matrici di connettività individuale e le caratteristiche non zero pesi ottenuti nel campione 1. Le previsioni risultanti sono state confrontate con i punteggi osservati di sensibilità al dolore basati sul QST calcolando l’errore assoluto medio (MAE), l’errore quadratico medio (MSE) e la varianza spiegata. I valori p-valori basati sulla permutazione sono stati ottenuti per tutte e tre le misure, utilizzando il pacchetto python mlxtend. Inoltre, il bootstrapping con copertura condizionale80 è stato utilizzato per fornire valori p per i pesi di connettività predittiva per aiutare l’interpretazione. Abbiamo costruito 10000 campioni di bootstrap (con sostituzione), con una dimensione pari al campione originale, costituito da coppie di dati del cervello e di risultato. Il modello predittivo con gli iperparametri ottimali è stato montato su ogni campione. P-valori non corretti sono stati calcolati per ogni connessione selezionata in base alla percentuale di pesi sotto o sopra lo zero, come ad esempio in rif. 30. Si noti che l’interpretazione di questi p-valori e intervalli di confidenza (Tabella supplementare 4) rimane limitata in quanto sono condizionati alla procedura di selezione delle caratteristiche.

Analisi dei fattori di confondimento

Per esplorare potenziali variabili di confondimento, i punteggi di sensibilità al dolore previsti (o le previsioni convalidate in modo incrociato nel caso del campione 1) sono stati contrapposti alla FD media e mediana, alla percentuale di volumi strofinati, alla pressione sanguigna sistolica e diastolica prima della risonanza magnetica e della misurazione QST (poiché la pressione sanguigna è stata precedentemente riportata81 come associata alla sensibilità al dolore meccanico), il ritardo temporale tra la risonanza magnetica e il test QST (per verificare la stabilità temporale della predizione), l’età, il sesso, l’IMC, il numero di giorni trascorsi dal primo giorno delle ultime mestruazioni, il consumo di alcol (unità/settimana), il livello di istruzione, l’ansia di stato e di tratto (STAI), punteggio dei sintomi depressivi (ADS-K), sensibilità al dolore auto-riferita (PSQ) e catastrofizzazione del dolore (PCS), stress percepito (PSQ20), qualità del sonno (PSQI), e soglie di rilevamento QST non nocive (CDT, WDT e MDT, se disponibili). Inoltre, nello studio 1 le previsioni sono state confrontate con i valori T50 e i livelli di GABA e Glutammato/Glutammina basati sulla spettroscopia MR nelle regioni cerebrali che elaborano il dolore (vedi rif. 8 per i dettagli). Le associazioni sono state testate con modelli lineari basati sulla permutazione.

Visualizzazione della rete predittiva

Le connessioni interregionali predittive evidenziate dai coefficienti di regressione non nulli della RPN-signature sono state visualizzate come un grafico a nastro utilizzando il R-pacchetto circlize (Fig. 3). Corrispondenti maschere individualizzate regione del cervello sono stati trasformati di nuovo allo spazio standard per creare uno studio-specifico mappa di probabilità regionale (che riflette la precisione co-registrazione e la variabilità individuale in morfologia). Le mappe di probabilità sono state moltiplicate per la somma dei coefficienti di regressione corrispondenti per creare una mappa regionale forza predittiva, che è stato poi visualizzato con FSLeyes e MRIcroGL.

(https://www.mccauslandcenter.sc.edu/mricrogl) (Fig. 3). L’analisi del coinvolgimento della rete su larga scala allo stato di riposo (come definito dall’atlante cerebrale MIST71) è stata eseguita riassumendo e trasformando in Z i valori dei voxel nelle sette regioni di interesse. La trama polare è stata fatta con il pacchetto R ggplot2.

Disponibilità del software

I punteggi RPN-signature possono essere calcolati sulla base di dataset strutturali e funzionali a riposo dallo strumento software con lo stesso nome. Lo strumento software RPN-signature consiste nella pipeline di elaborazione MRI descritta e nel modello predittivo basato sul connettoma funzionale. È disponibile come codice sorgente a https://github.com/spisakt/RPN-signature. Poiché il software segue la Brain Imaging Data Structure (BIDS)82 e la specifica BIDS-App, fornisce un’interfaccia standard a riga di comando e si basa sulla tecnologia Docker. L’immagine docker è depositata su Docker Hub: (https://cloud.docker.com/repository/docker/tspisak/rpn-signature) e non dipende da alcun software al di fuori dell’immagine del contenitore. Questo, insieme allo sviluppo basato sull’integrazione continua completamente trasparente e al tagging e al versioning automatizzati, migliora la disponibilità del software e supporta la riproducibilità dei risultati delle firme RPN.

Reporting summary

Più informazioni sul design della ricerca sono disponibili nel Nature Research Reporting Summary collegato a questo articolo.