Pain-free resting-state functional brain connectivity predicts individual pain sensitivity

General considerations

Studieutformningen fastställdes med noggrann hänsyn till de senaste rekommendationerna, kraven och standarderna för biomarkörer för neuroimaging50 (neuromarkörer) och motiverades av följande tankar.

Maximera den prediktiva prestandan

Vi använde oss av en standardiserad förbehandlingspipeline för att säkerställa optimal känslighet hos neuromarkören, eftersom tillräcklig effektstorlek är ett grundläggande krav för all klinisk användbarhet50. Vi använde oss av bildjustering med hög precision och införlivade individuell anatomi när vi extraherade fMRI-tidsseriedata. Dessutom antog vi nya rekommendationer och protokoll51 om minskning av artefakter och optimerade vårt arbetsflöde för att tillgodose de särskilda behoven för connectombaserad analys. Vi använde vårt egenutvecklade pythonprogramvarubibliotek Pipelines Utilising a Modular Inventory (PUMI, https://github.com/spisakt/PUMI) med öppen källkod, som bygger på nipype52 , ett gemenskapsbaserat Pythonprojekt som tillhandahåller ett enhetligt gränssnitt till befintlig programvara för neuroimaging och delvis återanvände kod från projekten med öppen källkod C-PAC53 och niworkflows54 . En metod för prediktiv modellering (maskininlärning) användes för att utnyttja de rika data som tillhandahålls av funktionella hjärnnätverk i vilostadiet och, potentiellt, dra nytta av fMRI-hyperacuitet55.

Bedömning av prediktiv förmåga under realistiska förhållanden

Vi använde oss av en förregistrerad, extern valideringsstrategi, som strikt skiljde på modellträning och prestandabedömning. För modellträning använde vi uteslutande data från studie 1. För validering genomförde vi två oberoende delstudier (studierna 2 och 3) på olika forskningscentra, med olika utrustning och olika forskningspersonal. Vi använde oss av en liberal anpassning av forskningsmiljöer, vilket gav utrymme för en rimlig heterogenitet i förfaranden, utrustning, bildsekvenser och språket i kommunikationen mellan deltagarna och forskarna i de olika studiecentren, vilket medförde en rimlig heterogenitet i valideringsförfarandet för att säkerställa generaliserbarheten.

Säkerställa att förutsägelsen drivs av neurala signaler och är specifik för smärtkänslighet

För att säkerställa att den föreslagna markören för smärtkänslighet verkligen drivs av neurala signaler som är associerade med smärtkänslighet, utvärderade vi korrelationen mellan den förutsagda poängsumman och olika fördefinierade (och förregistrerade) förväxlings- och valideringsvariabler.

Säkerställa resultatens tillgänglighet

Vi tillämpade en omfattande förregistrering och gjorde metodens källkod öppen källkod och fritt tillgänglig för samhället. Dessutom tillhandahåller vi en plattformsoberoende, lättanvänd docker-container, som ger möjlighet att använda vår prediktiva modell som en forskningsprodukt50, för att få out-of-the box smärtkänslighetsförutsägelser från alla lämpliga bilddatamängder.

Deltagare

Totalt N = 116 friska, unga frivilliga deltog i tre delstudier. Deltagarnas ålder och kön redovisas i kompletterande tabell 1. Studie 1 omfattade N1 = 39 deltagare (samma urval som i ref. 8). Den utfördes vid Ruhruniversitetet Bochum (Tyskland) av MZ och TSW och användes som träningsprov för den maskininlärningsbaserade förutsägelsen av smärtkänslighet och tjänade dessutom som grund för den interna valideringen av förutsägelsen. Studierna 2 och 3 (N2 = 48, N3 = 29) utfördes vid universitetssjukhuset i Essen (Tyskland) av FS och TS respektive vid universitetet i Szeged (Ungern) av BK och TK och användes som prov för extern validering. Inklusions- och exklusionskriterierna var i stort sett identiska i alla tre centra och anges i tabell 3. Rekryterings- och ersättningspolicy varierade mellan centrumen; deltagarna fick 20 €/h i studie 1 och 2 och ingen ersättning i studie 3.

Metallimplantat, icke avtagbar piercing, peacemaker, tatuering i huvud-/halsläge, graviditet eller känd klaustrofobi ansågs som kontraindikation för MR-mätning. Deltagarna var skyldiga att avstå från att konsumera koffein två timmar före experimenten (utom i studie 3) och från att konsumera alkohol på testdagen och föregående dag.

Studien genomfördes i enlighet med Helsingforsdeklarationen, uppfyller alla relevanta etiska regler för arbete med mänskliga deltagare och godkändes av de lokala eller nationella etikkommittéerna (registernummer: 4974-14, 18-8020-BO och 057617/2015/OTIG vid Ruhruniversitetet Bochum, Universitetssjukhuset Essen respektive ETT TUKEB Ungern). Alla deltagare gav skriftligt informerat samtycke före testning.

Bildtagning och kvantitativ sensorisk testning (QST) utfördes samma dag i studie 1 och med i genomsnitt 2-3 dagars mellanrum i studierna 2 och 3 (se kompletterande tabell 1 för detaljer). MRI-mätning föregick alltid QST-sessionen.

Mätningar-funktionell MRT

Anatomiska mätningar med hög upplösning och fMRI-mätningar i viloläge med öppna ögon förvärvades från alla deltagare. Skanningsparametrar (inklusive utrustning) varierade mellan olika centra och anges i tabell 4. Under mätningarna instruerades deltagarna att ligga stilla och avslappnade, utan att somna, och att undvika all rörelse. Skumkuddar, och i studierna 1 och 2, pneumatiska kuddar användes för att begränsa huvudets rörelser. Alla anatomiska MRI-mätningar screenades för tillfälliga fynd.

Mätningar-QST

Smärttrösklar för värme (HPT), kyla (CPT) och mekanisk smärta (MPT) förvärvades i enlighet med QST-protokollet28. Värme (WDT), kyla (CDT) och i studie 2 och 3 mekaniska (MDT) tröskelvärden för upptäckt av värme (WDT) erhölls som ytterligare kontrollmått. Alla sensoriska mätningar gjordes från den palmara vänstra underarmen, proximalt till handledskammen. Inom ramen för QST fastställs de termiska trösklarna med hjälp av en gränsvärdesmetod. I detta syfte applicerades ökande och minskande temperaturer på huden med en MSA-värmestimulator (Somedic, Hörby, Sverige) i studie 1 och Pathway-värmestimulatorer (Medoc Ltd., Ramat Yishai, Israel) i studierna 2 och 3. I alla studier användes ATS-termoder på en hudyta på 30 × 30 mm, med en utgångstemperatur på 32 °C. Deltagarna instruerades att ange smärtdebut genom att trycka på en knapp. För alla termiska tröskelvärden 6, i stället för 3 (som i det ursprungliga protokollet)28, utfördes stimulusupprepningar för att minska variansen mellan försökspersonerna. Dessutom uteslöts den första mätningen från analysen som ett teststimulus. HPT och CPT beräknades som aritmetiska medelvärden av de fem återstående tröskeltemperaturerna. MPT och MDT fastställdes med hjälp av en trappstegsmetod. Fem ökande och fem minskande tåg av pinprick-stimuli (MRC Systems, Heidelberg, Tyskland) applicerades på den palmara vänstra underarmen på ett omväxlande sätt, medan deltagaren instruerades att kategorisera stimuli som skadliga eller icke-skadliga. Tröskeln för mekanisk upptäckt bedömdes analogt med von Frey-filament-stimuleringar. MPT och MDT beräknades som den log-transformerade geometriska medelkraften som bestämdes i fem stigande och fallande trappstegs-tröskelkörningar.

Tillkommande åtgärder

Ålder, kön, självrapporterad längd, vikt och, för kvinnliga deltagare, datumet för den första dagen av den senaste menstruationen och användning av preventivmedel, registrerades före alla mätningar. Dessutom registrerades självrapporterad alkoholkonsumtion per vecka och utbildningsnivå (grundskola, gymnasium, universitet) för studierna 1 och 2. Före QST fyllde deltagarna i frågeformuläret för smärtkänslighet (Pain Sensitivity Questionnaire, PSQ)56, Pain Catastrophizing Scale (PCS)57, State-Trait Anxiety Inventory (STAI)58, den korta tyska versionen av Depression Scale (ADS-K, Center for Epidemiologic Studies)59 och, dessutom i studierna 2 och 3, Pittsburgh Sleep Quality Index (PSQI)60 och frågeformuläret för upplevd stress (PSQ20)61. I studierna 2 och 3 mättes blodtrycket både före MRT- och QST-mätningarna. För prov 1 fanns dessutom T50-värden tillgängliga från ett parallellt experiment som utfördes dagen före fMRI-testningen. T50 representerar den temperatur (i °C) som krävs för att framkalla en värmesmärta på 50 (på en skala från 0, ingen smärta till 100 outhärdlig smärta). T50-värdena erhölls genom en icke-linjär (andra ordningens polynom) interpolering av betyg som erhölls som svar på 15 toniska värmesmärtstimuli (varaktighet: 16 s) mellan 42,5 °C och 48 °C, som presenterades i ett pseudorandomiserat rutnätssökningsläge.

Beräkning av smärtkänslighet

Målvariabeln för förutsägelsen var ett enda sammansatt mått på individuell smärtkänslighet som sammanfattade HPT, CPT och MPT enligt definitionerna i ref. 8.

I studie 1 omvandlades HPT, CPT och MPT till Z-transformerade (medelvärdescentrerade och standardiserade) och HPT, liksom MPT, inverterades (multiplicerades med -1), så att högre Z-värden betecknade högre smärtkänslighet. Därefter beräknades det aritmetiska medelvärdet av de Z-transformerade variablerna för varje deltagare och definierades som poäng för smärtkänslighet. I studierna 2 och 3 tillämpades samma förfarande, förutom att Z-transformationen baserades på populationens medelvärde och standardavvikelse i studie 1, för att säkerställa att samma skala användes i alla studier. Extrema QST-värden definierades med hjälp av de normativa 95-procentiga percentilerna som redovisas i ref. 28. Deltagare som uppvisade extrema HPT-, CPT- eller MPT-värden i minst två av de tre modaliteterna exkluderades. Denna screening resulterade i att 0, 3 och 2 deltagare exkluderades i prov 1, 2 respektive 3 (kompletterande tabell 2).

fMRI-förbehandling

Då fMRI-baserad funktionell konnektivitet är känslig för rörelserelaterade artefakter i skannern62,63, är lämplig förbehandling och signalrensning nyckeln till framgångsrik konnektivitetsbaserad prediktion. Data från funktionell MRT i vilotillstånd förbehandlades på samma sätt i alla tre studierna. Det tillämpade, nipype-baserade arbetsflödet visas i kompletterande figur 1. Det använde sig av programvara för neuroimaging från tredje part, kod som anpassats från programvaruverktygen C-PAC53 och niworkflows54 samt interna pythonrutiner.

Hjärnutvinning från både de anatomiska och strukturella bilderna samt vävnadssegmentering från de anatomiska bilderna utfördes med FSL bet och fast64. Anatomiska bilder samregistrerades linjärt och icke-linjärt till den 1 mm upplösta MNI152 standard hjärnmallshjärnan med ANTs65 (se https://gist.github.com/spisakt/0caa7ec4bc18d3ed736d3a4e49da7415 för källkod).

Funktionella bilder samregistrerades till de anatomiska bilderna med den gränsbaserade registreringstekniken i FSL flirt. Alla resulterande transformationer sparades för vidare användning. Förbehandlingen av funktionella bilder skedde i det ursprungliga bildutrymmet, utan resampling. Realigneringsbaserad rörelsekorrigering utfördes med FSL mcflirt. De resulterande sex uppskattningarna av huvudets rörelse (3 rotationer, 3 translationer), deras kvadrerade versioner, deras derivat och de kvadrerade derivaten (kända som Friston-24-expansionen66 ) beräknades och sparades för störningskorrigering. Dessutom sammanfattades huvudrörelsen som tidsserier för ramvis förskjutning (FD) i enlighet med Powers metod63 för att användas vid censurering och uteslutning av data. Efter rörelsekorrigering dämpades outliers (t.ex. rörelsespikar) i tidsseriedata med hjälp av AFNI despike67. Föreningen av de eroderade white-matter-kartorna och ventrikelmaskerna transformerades till det ursprungliga funktionella utrymmet och användes för att extrahera brus-signal för anatomisk CompCor-korrigering68.

I ett störande regressionssteg avlägsnades 6 CompCor-parametrar (de 6 första huvudkomponenterna i tidsserierna för brus-regioner), Friston-24-rörlighetsparametrarna och den linjära trenden från tidsseriedata med en allmän linjär modell. På de kvarvarande uppgifterna utfördes temporal bandpassfiltrering med AFNI:s 3DBandpass för att behålla frekvensbandet 0,008-0,08 Hz. Den tidigare användningen av AFNI:s despike förväntas dämpa aliasing av kvarvarande rörelseartefakter i de angränsande tidsramarna under bandpassfiltreringen69. För att ytterligare dämpa effekterna av rörelseartefakter uteslöts potentiellt rörelsekontaminerade tidsramar, definierade med ett konservativt tröskelvärde FD > 0,15 mm, från datamaterialet (så kallad scrubbing av datamaterialet)70 . Deltagarna uteslöts från vidare analys om medelvärdet för FD översteg 0,15 mm eller om mer än 30 % av ramarna rensades bort. Detta resulterade i att 4, 8 och 7 deltagare i prov 1, 2 och 3 uteslöts (kompletterande tabell 2). Kvalitetskontroll (registreringskontroll, carpet-plots, se t.ex. kompletterande figurer 2-4) utfördes genom hela arbetsflödet.

Analys av funktionell konnektivitet

Den 122-parcellsversionen av den multiresolutionella funktionella hjärnatlasen MIST71 och gråämnesmasker som erhållits från den anatomiska bilden transformerades till det ursprungliga funktionella rummet. Denna atlas (konstruerad med BASC-metoden, dvs. bootstrap-analys av stabila kluster) har nyligen visat sig fungera väl vid konnektivitetsbaserad prediktiv modellering72. Atlasregionerna i det ursprungliga rummet maskerades med de gråa materialmasker som tidigare erhållits från den anatomiska bilden och omvandlats till det funktionella rummet. Med denna atlasindividualiseringsteknik kommer den slutliga regionala signalen att med hög sannolikhet härröra från grå-materievoxlar för varje försöksperson (vilket vi noggrant kontrollerade manuellt för alla försökspersoner), medan det med den konventionella metoden ingår ett varierande förhållande mellan grå- och vit-materievoxlar för varje försöksperson. Därför förväntas information från vävnadssegmenteringsprocessen minska variabiliteten från försöksperson till försöksperson (se kompletterande figur 5 för exempel). Voxeltidsserier medelvärdesbildades över dessa individualiserade MIST-regioner och, tillsammans med den genomsnittliga signalen för grå substans, bevarades för grafbaserad konnektivitetsanalys.

Regionala tidsserier ordnades i storskaliga funktionella moduler (definierade av MIST-atlasen med 7 paket) för visualiseringsändamål (fig. 1). Partiell korrelation beräknades över alla par av regioner (och global grå substans), vilket implementerades i pythonmodulen nilearn73. Partiella, snarare än enkla korrelationer användes för att utesluta indirekt konnektivitet74. Vår metod för grafmodellering säkerställde att den globala signalen för grå substans hanteras som en förvirring vid beräkningen av de partiella korrelationskoefficienterna, men samtidigt betraktades den också som en signal av intresse, eftersom den kan representera vaksamhetsrelaterade processer75. Partiella korrelationskoefficienter organiserades till 123 gånger 123 (122 regioner + global signal från grå substans) symmetriska konnektivitetsmatriser. Den övre triangeln i dessa matriser användes som funktionsutrymme för maskininlärningsbaserad prediktiv modellering.

Träning och validering av prediktiv modell

Data om funktionell konnektivitet i vila i hela hjärnan i studie 1 (N1 = 35, efter alla uteslutningar, som i ref. 8, kompletterande tabell 2) användes som inmatningsfunktionalitet (P = 7503 funktioner per deltagare) för att förutsäga individuella värderingar av smärtkänslighet, vilket ledde till en stor P-small N-inställning.

Vi konstruerade en pipeline för maskininlärning (https://github.com/spisakt/RPN-signature/blob/master/PAINTeR/model.py) i scikit-learn76, bestående av robust funktionsskalning (tar bort medianen och skalar med datakvantiteter), förval av funktioner77, val av de K bästa funktionerna med starkaste relationer till målvariabeln och en Elastic Net-regressionsmodell78 (en linjär modell med kombinerade L1- och L2-normer som regulariserare). Användningen av elastiskt nät var ett beslut som fattades före analysen. Vår huvudsakliga motivering för att välja elastiskt nät var att det gör det möjligt att optimera gleshet (L1 vs. L2-regularisering) som en hyperparameter, så att vi inte behövde göra några a-prioritära antaganden om glesheten hos den diskriminerande grundsanningens gleshet (se ref. 79 för motivering). Sammanfattningsvis kan man säga att fria hyperparametrar i pipelinen för maskininlärning var antalet förvalda funktioner (K), förhållandet mellan L1/L2-regulariseringen och vikten (alfa) för regulariseringen. Hyperparametrarna optimerades med ett rutnätssökförfarande och negativt medelkvadratfel som kostnadsfunktion. Värden för K varierade från 10 till 200 med ökningar på 5, och inkluderade för L1/L2-förhållandet för alfa. Hyperparameteroptimeringen utfördes i en korsvalidering utan deltagare (intern valideringsfas). Korsvalidering innefattade hela pipeline för maskininlärning för att undvika att införa beroenden mellan tränings- och testproverna. Observera att fMRI-förbehandlingen var oberoende mellan försökspersoner och därför inte ingick i korsvalidering. Optimala hyperparametrar befanns vara K = 25, L1/L2-förhållandet = 0,999 och alfa = 0,005.

Extern validering utfördes genom att tillämpa RPN-signaturen på fMRI-data från studierna 2 och 3 (N2 = 37, N3 = 19, efter uteslutningar, tilläggstabell 2), helt enkelt genom att tillämpa den funktionstransformation (skalning) som erhållits på prov 1 och sedan beräkna prickprodukten mellan individuella konnektivitetsmatriser och de icke-nollfunktionsvikter som erhållits i prov 1. De resulterande förutsägelserna jämfördes med de observerade QST-baserade värdena för smärtkänslighet genom att beräkna medelabsolut fel (MAE), medelkvadratfel (MSE) och förklarad varians. Permutationsbaserade p-värden erhölls för alla tre måtten med hjälp av pythonpaketet mlxtend. Dessutom användes bootstrapping med villkorlig täckning80 för att få fram p-värden för prediktiva konnektivitetsvikter för att underlätta tolkningen. Vi konstruerade 10 000 bootstrapprover (med ersättning), med en storlek som var lika stor som det ursprungliga provet, bestående av parade hjärn- och resultatdata. Den prediktiva modellen med de optimala hyperparametrarna anpassades till varje prov. Okorrigerade P-värden beräknades för varje valt samband baserat på andelen vikter under eller över noll, som i t.ex. ref. 30. Observera att tolkningen av dessa p-värden och konfidensintervall (kompletterande tabell 4) förblir begränsad eftersom de är betingade av förfarandet för funktionsurval.

Förstöraranalys

För att utforska potentiella förväxlingsvariabler kontrasterades de förutspådda smärtkänslighetssiffrorna (eller korsvaliderade förutsägelser när det gäller prov 1) mot medel- och medianvärdet av FD, procentandelen av skruvade volymer, systoliskt och diastoliskt blodtryck före både MRT- och QST-mätning (eftersom blodtrycket tidigare har rapporterats81 att vara förknippat med känslighet för mekanisk smärta), tidsfördröjning mellan MRT- och QST-mätning (för att testa om förutsägelsen är tidsmässigt stabil), ålder, kön, BMI, antal dagar sedan den första dagen av den senaste menstruationen, alkoholkonsumtion (enheter/vecka), utbildningsnivå, ångesttillstånd och ångestdrag (STAI), poäng för depressiva symtom (ADS-K), självrapporterad smärtkänslighet (PSQ) och smärtkatastrofering (PCS), upplevd stress (PSQ20), sömnkvalitet (PSQI) och tröskelvärden för detektering av icke-noxisk QST (CDT, WDT och MDT, om tillgängligt). Dessutom jämfördes förutsägelserna i studie 1 med T50-värden och MR-spektroskopibaserade GABA- och Glutamat/Glutamin-nivåer i smärtbearbetande hjärnregioner (se ref. 8 för detaljer). Associationer testades med permutationsbaserade linjära modeller.

Visualisering av det prediktiva nätverket

De prediktiva interregionala förbindelserna som lyfts fram av RPN-signaturens regressionskoefficienter som inte är noll visades som en bandplott med hjälp av R-paketet circlize (fig. 3). Motsvarande individualiserade hjärnregionsmasker transformerades tillbaka till standardrymd för att skapa en undersökningsspecifik regional sannolikhetskarta (som återspeglar samregistreringsnoggrannhet och individuell variabilitet i morfologi). Sannolikhetskartorna multiplicerades med summan av motsvarande regressionskoefficienter för att skapa en regional karta över prediktiv styrka, som sedan visualiserades med FSLeyes och MRIcroGL.

(https://www.mccauslandcenter.sc.edu/mricrogl) (fig. 3). Analysen av storskalig nätverksinvolvering i vilostadiet (enligt definitionen i hjärnatlasen MIST71) utfördes genom att sammanfatta och Z-transformera voxelvärdena över de sju regioner av intresse. Polar plot gjordes med R-paketet ggplot2.

Mjukvarutillgänglighet

RPN-signaturpoäng kan beräknas baserat på strukturella och funktionella datamängder i vilostadiet med hjälp av mjukvaruverktyget med samma namn. Mjukvaruverktyget RPN-signature består av den beskrivna MRI-bearbetningspipeline och den funktionella connectombaserade prediktiva modellen. Det finns tillgängligt som källkod på https://github.com/spisakt/RPN-signature. Eftersom programvaran följer specifikationen för Brain Imaging Data Structure (BIDS)82 och BIDS-App, tillhandahåller den ett standardkommandoradsgränssnitt och förlitar sig på Docker-teknik. Docker-avbildningen är deponerad på Docker Hub: (https://cloud.docker.com/repository/docker/tspisak/rpn-signature) och är inte beroende av någon programvara utanför containeravbildningen. Detta, tillsammans med den helt transparenta kontinuerliga integrationsbaserade utvecklingen och automatiserad taggning och versionering, förbättrar programvarans tillgänglighet och stöder reproducerbarheten av RPN-signaturens resultat.

Rapporteringssammanfattning

Förlängre information om forskningsdesign finns i Nature Research Reporting Summary som är länkad till denna artikel.