Smertefri funktionel hjerneforbindelse i hviletilstand forudsiger individuel smertefølsomhed

Generelle overvejelser

Undersøgelsens design blev udarbejdet under nøje hensyntagen til de seneste anbefalinger, krav og standarder for neuroimaging biomarkører50 (neuromarkører) og motiveret af følgende tanker.

Maximere prædiktiv ydeevne

Vi anvendte en standardiseret forbehandlingsrørledning for at sikre optimal følsomhed af neuromarkøren, da tilstrækkelig effektstørrelse er et grundlæggende krav for enhver klinisk anvendelighed50. Vi anvendte billedjustering med høj præcision og inkorporerede individuel anatomi ved udtrækning af fMRI-tidsseriedata. Desuden vedtog vi nylige anbefalinger og protokoller51 vedrørende artefaktreduktion og optimerede vores arbejdsgang for at opfylde de særlige behov for connectombaseret analyse. Vi anvendte vores egenudviklede, open source python-softwarebibliotek Pipelines Utilising a Modular Inventory (PUMI, https://github.com/spisakt/PUMI), som er baseret på nipype52 , et fællesskabsbaseret Python-projekt, der giver en ensartet grænseflade til eksisterende neuroimaging-software, og som til dels genbrugte kode fra open source-projekterne C-PAC53 og niworkflows54 . Der blev anvendt en tilgang til prædiktiv modellering (maskinlæring) for at udnytte de rige data, der leveres af funktionelle hjernenetværk i hviletilstand, og potentielt drage fordel af fMRI-hyperacuity55.

Vurdering af prædiktiv kraft under realistiske forhold

Vi anvendte en forudregistreret, ekstern valideringsstrategi, der strengt adskilte modeltræning og præstationsvurdering. Til modeltræning anvendte vi udelukkende data fra undersøgelse 1. Vi gennemførte to uafhængige delundersøgelser (undersøgelse 2 og 3) i forskellige forskningscentre, med forskelligt udstyr og forskerpersonale til validering. Vi anvendte en liberal tilpasning af forskningsmiljøer, der gav mulighed for en rimelig heterogenitet i procedurer, udstyr, billedsekvenser og sprog for kommunikation mellem deltagerne og forskerne på tværs af undersøgelsescentrene, hvilket introducerede en rimelig heterogenitet i valideringsproceduren for at sikre generaliserbarhed.

Sikre at forudsigelse er drevet af neurale signaler og er specifik for smertefølsomhed

For at sikre, at den foreslåede markør for smertefølsomhed faktisk er drevet af neurale signaler, der er forbundet med smertefølsomhed, evaluerede vi korrelationen af den forudsagte score med forskellige foruddefinerede (og forudregistrerede) confounder- og validatorvariabler.

Sikre tilgængelighed af resultater

Vi anvendte en omfattende præregistrering og gjorde kildekoden til metoden open-source og frit tilgængelig for samfundet. Desuden leverer vi en platformsuafhængig, let anvendelig docker-container, som giver mulighed for at bruge vores prædiktive model som et forskningsprodukt50 for at opnå out-of-the box-smertefølsomhedsforudsigelser fra alle relevante billeddannelsesdatasæt.

Deltagere

I alt N = 116 sunde, unge frivillige blev involveret i tre delundersøgelser. Deltagernes alder og køn er angivet i Supplerende tabel 1. Undersøgelse 1 omfattede N1 = 39 deltagere (den samme prøve som i ref. 8). Den blev udført på Ruhr-universitetet Bochum (Tyskland) af MZ og TSW og blev brugt som træningsprøve til den maskinlæringsbaserede forudsigelse af smertefølsomhed og tjente desuden som grundlag for den interne validering af forudsigelsen. Undersøgelse 2 og 3 (N2 = 48, N3 = 29) blev udført på University Hospital Essen (Tyskland) af henholdsvis FS og TS og på University of Szeged (Ungarn) af BK og TK og tjente som stikprøver til ekstern validering. Inklusions- og eksklusionskriterierne var stort set identiske i alle tre centre og er anført i tabel 3. Rekrutterings- og refusionspolitikker varierede på tværs af centrene; deltagerne modtog 20 €/h i undersøgelse 1 og 2 og ingen refusion i undersøgelse 3.

Metalimplantat, ikke-aftagelig piercing, peacemaker, tatovering i hoved/hals-stilling, graviditet eller kendt klaustrofobi blev betragtet som kontraindikation for MR-måling. Deltagerne skulle afholde sig fra at indtage koffein to timer før eksperimenterne (undtagen i undersøgelse 3) og fra at indtage alkohol på testdagen og den foregående dag.

Undersøgelsen blev gennemført i overensstemmelse med Helsinki-erklæringen, overholder alle relevante etiske regler for arbejde med menneskelige deltagere og blev godkendt af de lokale eller nationale etiske komitéer (Registernumre: 4974-14, 18-8020-BO og 057617/2015/OTIG ved henholdsvis Ruhr University Bochum, University Hospital Essen og ETT TUKEB Ungarn). Alle deltagere gav skriftligt informeret samtykke før testning.

Billeddannelse og kvantitativ sensorisk testning (QST) blev udført samme dag i undersøgelse 1 og i gennemsnit med 2-3 dages mellemrum i undersøgelse 2 og 3 (se Supplerende tabel 1 for detaljer). MRI-måling gik altid forud for QST-sessionen.

Målinger-funktionel MRI

Højopløselige anatomiske og åbenøjede hvilestatus fMRI-målinger blev erhvervet fra alle deltagere. Scanningsparametre (herunder udstyr) varierede på tværs af centrene og er anført i tabel 4. Under målingerne blev deltagerne instrueret om at ligge stille og afslappet, uden at falde i søvn, og undgå enhver bevægelse. Der blev anvendt skumpuder og i undersøgelse 1 og 2 pneumatiske puder for at begrænse hovedbevægelser. Alle anatomiske MRI-målinger blev screenet for tilfældige fund.

Målinger-QST

Smertetærskel for varme (HPT), kulde (CPT) og mekanisk (MPT) blev erhvervet i overensstemmelse med QST-protokollen28. Varme (WDT), kulde (CDT) og i undersøgelse 2 og undersøgelse 3 blev der opnået mekaniske (MDT) detektionstærskelværdier som yderligere kontrolforanstaltninger. Alle sensoriske målinger blev foretaget fra den venstre underarm i håndfladen, proximalt i forhold til håndledskammen. Inden for QST-rammen bestemmes termiske tærskelværdier ved hjælp af en metode med grænser. Med henblik herpå blev der påført stigende og faldende temperaturer på huden med en MSA-termisk stimulator (Somedic, Hörby, Sverige) i undersøgelse 1 og Pathway-termiske stimulatorer (Medoc Ltd., Ramat Yishai, Israel) i undersøgelse 2 og 3. I alle undersøgelser blev der anvendt ATS-termoder på en hudoverflade på 30 × 30 mm med en baseline-temperatur på 32 °C. Deltagerne blev instrueret om at angive smertens indtræden ved at trykke på en knap. For alle termiske tærskler 6, snarere end 3 (som i den oprindelige protokol)28, blev der udført stimulus gentagelser for at reducere variansen mellem forsøgspersoner. Desuden blev den første måling udelukket fra analysen som en teststimulus. HPT og CPT blev beregnet som de aritmetiske gennemsnit af de fem resterende tærskeltemperaturer. MPT’er og MDT’er blev bestemt ved hjælp af en trappemetode. Fem stigende og fem faldende træk af pinprick-stimuli (MRC Systems, Heidelberg, Tyskland) blev påført den venstre underarm skiftevis på palmar venstre underarm, mens deltageren blev instrueret om at kategorisere stimuli som skadelige eller ikke-skadelige. Mekanisk detektionstærskel blev vurderet analogt med von Frey-filament-stimuleringer. MPT og MDT blev beregnet som den log-transformerede geometriske middelkraft bestemt i fem opadgående og nedadgående trappetræk-tærskel-løb.

Tilbagevendende foranstaltninger

Alder, køn, selvrapporteret højde, vægt og, for kvindelige deltagere, datoen for den første dag af den sidste menstruation og brugen af præventionsmidler, blev registreret forud for alle målinger. Desuden blev selvrapporteret ugentligt alkoholforbrug og uddannelsesniveau (grundskole, gymnasium, universitet) registreret i forbindelse med undersøgelse 1 og 2. Før QST udfyldte deltagerne Pain Sensitivity Questionnaire (PSQ)56, Pain Catastrophizing Scale (PCS)57, State-Trait Anxiety Inventory (STAI)58, den korte tyske version af Depression Scale (ADS-K, Center for Epidemiologic Studies)59 og, yderligere i Undersøgelse 2 og 3, Pittsburgh Sleep Quality Index (PSQI)60 og det opfattede stressspørgeskema (PSQ20)61. I undersøgelse 2 og 3 blev blodtrykket målt både før MRI- og QST-målingerne. Desuden var T50-værdierne for prøve 1 tilgængelige fra et parallelt eksperiment, der blev udført dagen før fMRI-testningen. T50 repræsenterer den temperatur (i °C), der er nødvendig for at fremkalde en varme-smerte-rating på 50 (på en skala fra 0, ingen smerte til 100 uudholdelig smerte). T50-værdier blev opnået fra en ikke-lineær (andenordenspolynomial) interpolation af vurderinger opnået som reaktion på 15 toniske varme-smerte-stimuli (varighed: 16 s) mellem 42.5 °C og 48 °C, præsenteret i pseudo-randomiseret gitter-søgning.

Beregning af smertefølsomhed

Målvariablen for forudsigelsen var en enkelt sammensat måling af individuel smertefølsomhed, der opsummerer HPT, CPT og MPT som defineret i ref. 8.

I undersøgelse 1 blev HPT, CPT og MPT Z-transformeret (middelcentreret og standardiseret), og HPT såvel som MPT blev inverteret (ganget med -1), således at højere Z-værdier betød højere smertefølsomhed. Derefter blev det aritmetiske gennemsnit af de Z-transformerede variabler beregnet for hver deltager og defineret som smertefølsomhedsscore. I undersøgelse 2 og 3 blev den samme procedure anvendt, bortset fra at Z-transformationen blev baseret på populationsmiddelværdien og standardafvigelsen i undersøgelse 1 for at sikre, at der blev anvendt den samme skala på tværs af undersøgelserne. Ekstreme QST-værdier blev defineret ved hjælp af de normative 95 %-percentiler, der er rapporteret i ref. 28; deltagere, der viste ekstreme HPT-, CPT- eller MPT-værdier i mindst to af de tre modaliteter, blev udelukket. Denne screening resulterede i udelukkelse af 0, 3 og 2 deltagere i henholdsvis prøve 1, 2 og 3 (Supplerende tabel 2).

fMRI-forbehandling

Som fMRI-baseret funktionel konnektivitet er modtagelig for bevægelsesartefakter i scanneren62,63, er passende forbehandling og signalrensning nøglen til en vellykket konnektivitetsbaseret forudsigelse. Funktionelle MRT-data i hviletilstand blev forbehandlet på samme måde i alle tre undersøgelser. Den anvendte, nipype-baserede arbejdsgang er afbildet på supplerende figur 1. Den anvendte neuroimaging-software fra tredjepart, kode tilpasset fra softwareværktøjerne C-PAC53 og niworkflows54 og interne python-rutiner.

Hjerneudtræk fra både de anatomiske og strukturelle billeder samt vævssegmentering fra de anatomiske billeder blev udført med FSL bet og fast64. Anatomiske billeder blev lineært og ikke-lineært samregistreret til den 1 mm opløselige MNI152 standard hjerne-skabelon hjerne med ANTs65 (se https://gist.github.com/spisakt/0caa7ec4bc18d3ed736d3a4e49da7415 for kildekode).

Funktionelle billeder blev samregistreret til de anatomiske billeder med den grænsebaserede registreringsteknik i FSL flirt. Alle resulterende transformationer blev gemt til videre brug. Forbehandlingen af funktionelle billeder skete i det oprindelige billedrum uden resampling. Realignmentbaseret bevægelseskorrektion blev udført med FSL mcflirt. De resulterende seks hovedbevægelsesestimater (3 rotationer, 3 translokationer), deres kvadrerede versioner, deres derivater og de kvadrerede derivater (kendt som Friston-24-udvidelsen66 ) blev beregnet og gemt med henblik på korrektion af uhensigtsmæssigheder. Desuden blev hovedbevægelsen opsummeret som tidsrækker for framewise displacement (FD) i henhold til Powers metode63 , som skulle anvendes til datacensurering og udelukkelse. Efter bevægelseskorrektion blev outliers (f.eks. bevægelsesspidser) i tidsseriedata dæmpet ved hjælp af AFNI despike67. Foreningen af de eroderede white-matter-kort og ventrikelmasker blev transformeret til det native funktionelle rum og brugt til at udtrække støjsignal til anatomisk CompCor-korrektion68.

I et nuisance-regressionstrin blev 6 CompCor-parametre (de 6 første hovedkomponenter i støjregionens tidsserier), Friston-24-bevægelsesparametrene og den lineære trend fjernet fra tidsseriedataene med en generel lineær model. På de resterende data blev der udført tidsmæssig båndpasfiltrering med AFNI’s 3DBandpass for at bevare frekvensbåndet 0,008-0,08 Hz. Den forudgående brug af AFNI’s despike forventes at dæmpe aliasing af resterende bevægelsesartefakter i de tilstødende tidsrammer under båndpasfiltreringen69. For yderligere at mindske virkningen af bevægelsesartefakter blev potentielt bevægelsesforurenede tidsrammer, defineret ved en konservativ tærskelværdi FD > 0,15 mm, udeladt fra dataene (såkaldt scrubbing af dataene)70. Deltagerne blev udelukket fra yderligere analyse, hvis den gennemsnitlige FD oversteg 0,15 mm, eller hvis mere end 30 % af rammene blev scrubbed. Dette resulterede i udelukkelse af henholdsvis 4, 8 og 7 deltagere i stikprøve 1, 2 og 3 (Supplerende tabel 2). Kvalitetskontrol (registreringskontrol, tæppe-plots, se f.eks. Supplerende figurer 2-4) blev udført i hele arbejdsgangen.

Funktionel konnektivitetsanalyse

Den 122-parcel version af MIST71 multiopløsnings funktionelle hjerneatlas og gråstofmasker opnået fra det anatomiske billede blev transformeret til det native funktionelle rum. Dette atlas (konstrueret ved hjælp af BASC-metoden, dvs. bootstrap-analyse af stabile klynger) blev for nylig vist at fungere godt i forbindelse med konnektivitetsbaseret prædiktiv modellering72. Atlasregioner i det native rum blev maskeret med de gråstofmasker, der tidligere blev opnået fra det anatomiske billede og transformeret til det funktionelle rum. Med denne atlas-individualiseringsteknik vil det endelige regionale signal stamme – med høj sandsynlighed – fra voxler med grå substans for hver forsøgsperson (hvilket vi omhyggeligt kontrollerede manuelt for alle forsøgspersoner), mens der med den konventionelle metode indgår et variabelt forhold mellem voxler med grå substans og voxler med hvid substans for hver forsøgsperson. Derfor forventes det, at indlæsning af oplysninger fra vævssegmenteringsprocessen vil mindske variabiliteten fra emne til emne (se supplerende fig. 5 for eksempler). Voxel-tidsserier blev gennemsnit over disse individualiserede MIST-regioner og sammen med det gennemsnitlige gråstofsignal tilbageholdt til grafbaseret konnektivitetsanalyse.

Regionale tidsserier blev ordnet i storskala funktionelle moduler (defineret af 7-parcel MIST-atlaset) til visualiseringsformål (Fig. 1). Partiel korrelation blev beregnet på tværs af alle par af regioner (og global grå substans), som implementeret i nilearn73 pythonmodulet. Delvise, snarere end simple korrelationer blev brugt til at udelukke indirekte konnektivitet74. Vores grafmodelleringstilgang sikrede, at det globale grå stof-signal håndteres som en forvirring under beregningen af de partielle korrelationskoefficienter, men samtidig blev det også betragtet som et signal af interesse, da det kan repræsentere vaktualitetsrelaterede processer75. Delvise korrelationskoefficienter blev organiseret til 123 gange 123 (122 regioner + globalt gråstofsignal) symmetriske forbindelsesmatricer. Den øverste trekant af disse matricer blev brugt som funktionsrum til maskinlæringsbaseret prædiktiv modellering.

Prædiktiv modeltræning og validering

Hele hjernen i hviletilstand funktionelle konnektivitetsdata fra undersøgelse 1 (N1 = 35, efter alle udelukkelser, som i ref. 8, Supplerende tabel 2) blev brugt som input-funktionalitetsrum (P = 7503 funktioner pr. Deltager) til at forudsige individuelle smertefølsomhedsscorer, hvilket førte til en stor P-små N-indstilling.

Vi konstruerede en pipeline til maskinlæring (https://github.com/spisakt/RPN-signature/blob/master/PAINTeR/model.py) i scikit-learn76, der består af robust funktionsskalering (fjerner medianen og skalerer med datakvantiler), forvalg af funktioner77, udvælgelse af de K bedste funktioner med stærkeste relationer til målvariablen og en Elastic Net-regressionsmodel78 (en lineær model med kombinerede L1- og L2-normer som regulariser). Brugen af elastisk net var en beslutning, der blev truffet forud for analysen. Vores hovedmotivation for at vælge elastisk net var, at det giver mulighed for at optimere sparsomhed (L1 vs. L2-regularisering) som en hyperparameter, så vi ikke behøvede at foretage nogen a-præceptive antagelser om sparsomheden af den diskriminerende grundsandhed (se ref. 79 for begrundelse). For at opsummere kan man sige, at frie hyperparametre i pipelinen til maskinlæring var antallet af på forhånd udvalgte features (K), forholdet mellem L1/L2-regularisering og regulariseringens vægt (alpha). Hyperparametrene blev optimeret med en gitter-søgningsprocedure og negativ gennemsnitlig kvadreret fejl som omkostningsfunktion. Værdierne for K varierede fra 10 til 200 med stigninger på 5, og for L1/L2-forholdet for alpha blev inkluderet. Hyperparameteroptimering blev udført i en leave-one-participant-out cross-validation (intern valideringsfase). Krydsvalideringen omfattede den komplette pipeline for maskinlæring for at undgå at indføre afhængigheder mellem trænings- og testprøverne. Bemærk, at fMRI-forbehandlingen var uafhængig mellem forsøgspersoner og derfor ikke indgik i krydsvalideringen. Optimale hyperparametre blev fundet at være K = 25, L1/L2-forhold = 0,999 og alpha = 0,005.

Ekstern validering blev udført ved at anvende RPN-signaturen på fMRI-dataene fra undersøgelse 2 og 3 (N2 = 37, N3 = 19, efter udelukkelser, Supplerende tabel 2), simpelthen ved at anvende den funktionstransformation (skalering), der er opnået på prøve 1, og derefter beregne prikproduktet mellem individuelle konnektivitetsmatricer og de ikke-nul funktionsvægte, der er opnået i prøve 1. De resulterende forudsigelser blev sammenlignet med de observerede QST-baserede smertefølsomhedsscorer ved at beregne den gennemsnitlige absolutte fejl (MAE), den gennemsnitlige kvadrerede fejl (MSE) og den forklarede varians. Permutationsbaserede p-værdier blev opnået for alle tre mål ved hjælp af pythonpakken mlxtend. Desuden blev bootstrapping med betinget dækning80 anvendt til at tilvejebringe p-værdier for prædiktive konnektivitetsvægte for at lette fortolkningen. Vi konstruerede 10000 bootstrapprøver (med udskiftning), med en størrelse svarende til den oprindelige prøve, bestående af parrede hjerne- og resultatdata. Den prædiktive model med de optimale hyperparametre blev tilpasset til hver prøve. Ukorrigerede P-værdier blev beregnet for hver udvalgt forbindelse baseret på andelen af vægte under eller over nul, som i f.eks. ref. 30. Bemærk, at fortolkningen af disse p-værdier og konfidensintervaller (Supplerende tabel 4) forbliver begrænset, da de er betinget af proceduren for udvælgelse af funktioner.

Forsageranalyse

For at undersøge potentielle forstyrrende variabler blev de forudsagte smertefølsomheds-scorer (eller krydsvaliderede forudsigelser i tilfælde af prøve 1) sat i kontrast til gennemsnitlig og median FD, procentdelen af skrubbede volumener, systolisk og diastolisk blodtryk før både MRI- og QST-måling (da blodtryk tidligere er blevet rapporteret81 som værende forbundet med følsomhed over for mekanisk smerte), tidsforsinkelsen mellem MRI- og QST-måling (for at teste for tidsmæssig stabilitet af forudsigelsen), alder, køn, BMI, antal dage siden den første dag af sidste menstruation, alkoholforbrug (enheder/uge), uddannelsesniveau, tilstandsangst og trækangst (STAI), score for depressive symptomer (ADS-K), selvrapporteret smertefølsomhed (PSQ) og smertekatastrofering (PCS), opfattet stress (PSQ20), søvnkvalitet (PSQI) og ikke-noxiske QST-detektionstærskler (CDT, WDT og MDT, hvis de foreligger). Desuden blev forudsigelser i Studie 1 sammenlignet med T50-værdier og MR-spektroskopi-baserede GABA- og Glutamat/Glutamin-niveauer i smertebearbejdende hjerneområder (se ref. 8 for detaljer). Associeringer blev testet med permutationsbaserede lineære modeller.

Visualisering af det prædiktive netværk

De prædiktive interregionale forbindelser fremhævet af de ikke-nul regressionskoefficienter af RPN-signaturen blev vist som et båndplot ved hjælp af R-pakken circlize (Fig. 3). Tilsvarende individualiserede hjerneområde-masker blev transformeret tilbage til standardrum for at skabe et undersøgelsesspecifikt regionalt sandsynlighedskort (der afspejler co-registreringsnøjagtighed og individuel variabilitet i morfologi). Sandsynlighedskort blev multipliceret med summen af de tilsvarende regressionskoefficienter for at skabe et regionalt prædiktionsstyrkekort, som derefter blev visualiseret med FSLeyes og MRIcroGL.

(https://www.mccauslandcenter.sc.edu/mricrogl) (Fig. 3). Analysen af storskala hviletilstandsnetværksinvolvering (som defineret af MIST71-hjerneatlaset) blev udført ved at sammenfatte og Z-transformere voxelværdierne på tværs af de syv regioner af interesse. Polar plot blev lavet med R-pakken ggplot2.

Software tilgængelighed

RPN-signaturscorerne kan beregnes på grundlag af strukturelle og resting-state funktionelle datasæt ved hjælp af softwareværktøjet med samme navn. RPN-signatur-softwareværktøjet består af den beskrevne MRI-behandlingspipeline og den funktionelle connectombaserede forudsigelsesmodel. Det er tilgængeligt som kildekode på https://github.com/spisakt/RPN-signature. Da softwaren følger Brain Imaging Data Structure (BIDS)82 og BIDS-App-specifikationen, giver den en standardkommandolinjeinterface og er afhængig af Docker-teknologien. Docker-aftrykket er deponeret på Docker Hub: (https://cloud.docker.com/repository/docker/tspisak/rpn-signature) og er ikke afhængig af nogen software uden for containeraftrykket. Dette, sammen med den fuldt gennemsigtige kontinuerlige integrationsbaserede udvikling og automatiseret tagging og versionering, forbedrer softwarens tilgængelighed og understøtter reproducerbarheden af RPN-signaturresultater.

Rapporteringsresumé

Der findes yderligere oplysninger om forskningsdesign i Nature Research Reporting Summary, der er knyttet til denne artikel.