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

Algemene overwegingen

De studieopzet is tot stand gekomen met zorgvuldige overweging van recente aanbevelingen, eisen en standaarden voor neuroimaging biomarkers50 (neuromarkers) en ingegeven door de volgende gedachten.

Maximaliseer voorspellende prestaties

We gebruikten een gestandaardiseerde voorbewerkingspijplijn om optimale gevoeligheid van de neuromarker te garanderen, aangezien voldoende effectgrootte een basisvereiste is voor elke klinische bruikbaarheid50. We gebruikten hoge-precisie beelduitlijning, waarbij individuele anatomie werd opgenomen bij het extraheren van fMRI-tijdenreeksgegevens. Bovendien hebben we aangenomen recente aanbevelingen en protocollen51 met betrekking tot artefact reductie en geoptimaliseerd onze workflow aan de speciale behoeften van connectoom-gebaseerde analyse te voldoen. We gebruikten onze in-house ontwikkelde, open-source python software bibliotheek Pipelines Utilising a Modular Inventory (PUMI, https://github.com/spisakt/PUMI), die is gebaseerd op nipype52, een community-based Python-project dat een uniforme interface biedt voor bestaande neuroimaging software en, gedeeltelijk, code hergebruikt van de C-PAC53 en de niworkflows54 open-source projecten. Een voorspellende modellering (machine learning) aanpak werd gebruikt om de rijke gegevens van resting-state functionele hersennetwerken te exploiteren en, potentieel, te profiteren van fMRI hyperacuiteit55.

Beoordeling van voorspellend vermogen onder realistische omstandigheden

We gebruikten een vooraf geregistreerde, externe validatie strategie, die een strikte scheiding maakte tussen model training en prestatie beoordeling. Voor de modeltraining gebruikten we uitsluitend gegevens uit studie 1. Voor de validatie voerden wij twee onafhankelijke deelstudies uit (Studie 2 en 3) in verschillende onderzoekscentra, met verschillende apparatuur en verschillend onderzoekspersoneel. We gebruikten een liberale afstemming van onderzoeksinstellingen, waardoor een redelijke heterogeniteit in procedures, apparatuur, beeldvorming sequenties, taal van de deelnemer-onderzoeker communicatie in studie-centra, waardoor een redelijke heterogeniteit in de validatieprocedure om de generaliseerbaarheid te waarborgen.

Vergewis u ervan dat de voorspelling wordt aangedreven door neurale signalen en specifiek is voor pijngevoeligheid

Om er zeker van te zijn dat de voorgestelde marker van pijngevoeligheid inderdaad wordt aangedreven door neurale signalen die verband houden met pijngevoeligheid, evalueerden we de correlatie van de voorspelde score met verschillende vooraf gedefinieerde (en vooraf geregistreerde) confounder en validator variabelen.

Zorg voor toegankelijkheid van de resultaten

We pasten een uitgebreide preregistratie toe en maakten de broncode van de methode open-source en vrij beschikbaar voor de gemeenschap. Bovendien bieden we een platformonafhankelijke, eenvoudig te gebruiken docker-container, die de mogelijkheid biedt om ons voorspellende model te gebruiken als een onderzoeksproduct50, om out-of-the-box pijngevoeligheidsvoorspellingen te verkrijgen van elke geschikte beeldvormingsdataset.

Deelnemers

Er waren in totaal N = 116 gezonde, jonge vrijwilligers betrokken bij drie deelstudies. Leeftijd en geslacht van de deelnemers wordt gerapporteerd in de aanvullende tabel 1. Studie 1 betrof N1 = 39 deelnemers (dezelfde steekproef als in ref. 8). Deze studie werd uitgevoerd aan de Ruhr Universiteit Bochum (Duitsland) door MZ en TSW en werd gebruikt als het trainingsmonster voor de machine-learning gebaseerde voorspelling van pijngevoeligheid en diende bovendien als basis voor de interne validatie van de voorspelling. Studies 2 en 3 (N2 = 48, N3 = 29) werden uitgevoerd in het Universitair Ziekenhuis Essen (Duitsland) door FS en TS en aan de Universiteit van Szeged (Hongarije) door BK en TK, respectievelijk, en dienden als steekproeven voor de externe validatie. De inclusie- en exclusiecriteria waren in alle drie de centra grotendeels identiek en staan vermeld in tabel 3. Het wervings- en vergoedingsbeleid verschilde tussen de centra; deelnemers ontvingen 20 €/h in studie 1 en 2 en geen vergoeding in studie 3.

Metaalimplantaat, niet-verwijderbare piercing, peacemaker, tatoeage in hoofd-/nekpositie, zwangerschap of bekende claustrofobie werden beschouwd als contra-indicatie voor MR-meting. De deelnemers werden verzocht zich te onthouden van de consumptie van cafeïne twee uur voor de experimenten (behalve in studie 3) en van de consumptie van alcohol op de dag van de test en de vorige dag.

De studie werd uitgevoerd in overeenstemming met de Verklaring van Helsinki, voldoet aan alle relevante ethische voorschriften voor het werken met menselijke deelnemers en werd goedgekeurd door de lokale of nationale ethische commissies (Registernummers: 4974-14, 18-8020-BO en 057617/2015/OTIG aan respectievelijk de Ruhruniversiteit Bochum, het Universitair Ziekenhuis Essen en ETT TUKEB Hongarije). Alle deelnemers gaven schriftelijke geïnformeerde toestemming vóór het testen.

Imaging en kwantitatieve sensorische testen (QST) werden uitgevoerd op dezelfde dag in studie 1 en gemiddeld 2-3 dagen uit elkaar in studies 2 en 3 (zie aanvullende tabel 1 voor details). MRI metingen gingen altijd vooraf aan de QST sessie.

Metingen-functionele MRI

Hoge-resolutie anatomische en open-oog rust-staat fMRI metingen werden verkregen van alle deelnemers. De scanparameters (inclusief apparatuur) verschilden per centrum en staan vermeld in tabel 4. Tijdens de metingen werden de deelnemers geïnstrueerd om stil en ontspannen te liggen, zonder in slaap te vallen, en elke beweging te vermijden. Er werden schuimkussens, en in studies 1 en 2, pneumatische kussens gebruikt om hoofdbewegingen te beperken. Alle anatomische MRI-metingen werden gescreend op incidentele bevindingen.

Metingen-QST

Hitte (HPT), koude (CPT) en mechanische (MPT) pijndrempels werden verkregen volgens het QST-protocol28. Warmte (WDT), koude (CDT) en in studie 2 en studie 3, mechanische (MDT) waarnemingsdrempels werden verkregen als extra controlemaatregelen. Alle sensorische metingen werden verkregen van de palmaire linker onderarm, proximaal aan de pols crest. In het kader van QST worden de thermische drempels bepaald met behulp van een grenswaardenmethode. Daartoe werden stijgende en dalende temperaturen op de huid aangebracht met een MSA thermische stimulator (Somedic, Hörby, Zweden) in studie 1 en Pathway thermische stimulatoren (Medoc Ltd., Ramat Yishai, Israël) in studies 2 en 3. In alle studies werden ATS thermoden gebruikt op een huidoppervlak van 30 × 30 mm, met een basistemperatuur van 32 °C. Deelnemers werden geïnstrueerd om het begin van de pijn aan te geven door op een knop te drukken. Voor alle thermische drempels 6, in plaats van 3 (zoals in het oorspronkelijke protocol)28, stimulus herhalingen werden uitgevoerd om tussen-subject variantie te verminderen. Bovendien werd de eerste meting als teststimulus bij de analyse buiten beschouwing gelaten. HPT en CPT werden berekend als de rekenkundige gemiddelden van de vijf resterende drempeltemperaturen. MPT’s en MDT’s werden bepaald met behulp van een trapmethode. Vijf oplopende en vijf aflopende series van pinprikstimuli (MRC Systems, Heidelberg, Duitsland) werden afwisselend toegepast op de linker handpalm onderarm, terwijl de deelnemer de instructie kreeg om de stimuli te categoriseren als al dan niet pijnlijk. Mechanische detectie drempel werd beoordeeld analoog met von Frey filament stimuli. MPT en MDT werden berekend als de log-getransformeerde geometrische gemiddelde kracht bepaald in vijf oplopende en aflopende trap-drempel-runs.

Aanvullende maatregelen

Leeftijd, geslacht, zelf-gerapporteerde lengte, gewicht en, voor vrouwelijke deelnemers, de datum van de eerste dag van de laatste menstruatie en het gebruik van voorbehoedsmiddelen, werd voorafgaand aan alle metingen genoteerd. Bovendien werden de zelfgerapporteerde wekelijkse alcoholconsumptie en het opleidingsniveau (lagere school, middelbare school, universiteit) geregistreerd voor onderzoek 1 en 2. Vóór de QST vulden de deelnemers de Pain Sensitivity Questionnaire (PSQ)56, de Pain Catastrophizing Scale (PCS)57, de State-Trait Anxiety Inventory (STAI)58, de korte Duitse versie van de Depression Scale (ADS-K, Center for Epidemiologic Studies)59 en, aanvullend in studies 2 en 3, de Pittsburgh Sleep Quality Index (PSQI)60 en de waargenomen stress vragenlijst (PSQ20)61 in. In studies 2 en 3 werd de bloeddruk zowel vóór de MRI- als de QST-metingen gemeten. Bovendien waren voor monster 1 T50-waarden beschikbaar uit een parallel experiment dat was uitgevoerd op de dag vóór de fMRI-tests. T50 vertegenwoordigt de temperatuur (in °C) die nodig is om een warmte-pijn waardering van 50 te induceren (op een schaal van 0, geen pijn tot 100 ondraaglijke pijn). T50 waarden werden verkregen uit een niet-lineaire (tweede-orde polynomiale) interpolatie van beoordelingen verkregen in antwoord op 15 tonische warmte-pijn stimuli (duur: 16 s) tussen 42,5 °C en 48 °C, gepresenteerd in pseudo-gerandomiseerde grid-search mode.

Berekening van pijngevoeligheid

De doelvariabele voor de voorspelling was een enkele samengestelde maat van individuele pijngevoeligheid die HPT, CPT en MPT samenvatte zoals gedefinieerd in ref. 8.

In studie 1 werden HPT, CPT en MPT Z-getransformeerd (gemiddeld gecentreerd en gestandaardiseerd) en werden HPT, evenals MPT omgekeerd (vermenigvuldigd met -1), zodat hogere Z-waarden een hogere pijngevoeligheid aangaven. Vervolgens werd het rekenkundig gemiddelde van de Z-getransformeerde variabelen berekend voor elke deelnemer en gedefinieerd als pijngevoeligheidsscore. In Studie 2 en 3 werd dezelfde procedure toegepast, behalve dat de Z-transformatie gebaseerd werd op het populatie-gemiddelde en de standaardafwijking van Studie 1, om er zeker van te zijn dat dezelfde schaal werd gebruikt in alle studies. Extreme QST-waarden werden gedefinieerd met gebruikmaking van de normatieve 95%-percentielen zoals gerapporteerd in ref. 28; deelnemers met extreme HPT-, CPT- of MPT-waarden in ten minste twee van de drie modaliteiten werden uitgesloten. Deze screening resulteerde in het uitsluiten van 0, 3 en 2 deelnemers in Samples 1, 2 en 3, respectievelijk (Supplementary Table 2).

fMRI preprocessing

Aangezien fMRI-gebaseerde functionele connectiviteit gevoelig is voor in-scanner bewegingsartefacten62,63, is de juiste preprocessing en signaal opschoning de sleutel tot succesvolle connectiviteit-gebaseerde voorspelling. Resting-state functionele MRI-gegevens werden op identieke wijze voorbewerkt in alle drie studies. De toegepaste, nipype-gebaseerde workflow is afgebeeld op aanvullende Fig. 1. Het maakte gebruik van derden neuroimaging software, code aangepast van de software tools C-PAC53 en niworkflows54, en in-house python routines.

Brain extractie van zowel de anatomische en de structurele beelden, maar ook als weefsel-segmentatie van de anatomische beelden werd uitgevoerd met FSL inzet en fast64. Anatomische beelden werden lineair en niet-lineair co-registered aan de 1 mm-resolutie MNI152 standaard hersenen template hersenen met ANTs65 (zie https://gist.github.com/spisakt/0caa7ec4bc18d3ed736d3a4e49da7415 voor de broncode).

Functionele beelden werden co-registered aan de anatomische beelden met de grens-gebaseerde registratie techniek van FSL flirt. Alle resulterende transformaties werden opgeslagen voor verder gebruik. De pre-processing van functionele beelden gebeurde in de oorspronkelijke beeldruimte, zonder resampling. Op herschikking gebaseerde bewegingscorrectie werd uitgevoerd met FSL mcflirt. De resulterende zes schattingen van de hoofdbeweging (3 rotaties, 3 translaties), hun gekwadrateerde versies, hun afgeleiden en de gekwadrateerde afgeleiden (bekend als de Friston-24-expansie66) werden berekend en opgeslagen voor hindercorrectie. Bovendien werd de beweging van het hoofd samengevat als framewise verplaatsing (FD) tijdreeks, volgens de methode van Power63, om te worden gebruikt bij het censureren en uitsluiten van gegevens. Na bewegingscorrectie werden uitbijters (bv. bewegingspieken) in de tijdreeksgegevens verzwakt met AFNI despike67. De unie van de geërodeerde witte-stof kaarten en ventrikel maskers werden getransformeerd naar de inheemse functionele ruimte en gebruikt voor het extraheren van ruis-signaal voor anatomische CompCor correctie68.

In een hinderlijke regressie stap, 6 CompCor parameters (de 6 eerste hoofdcomponenten van de ruis-regio tijdreeksen), de Friston-24 beweging parameters en de lineaire trend werden verwijderd uit de tijdreeks gegevens met een algemeen lineair model. Op de residuele gegevens, werd temporele bandpass filtering uitgevoerd met AFNI’s 3DBandpass om de 0.008-0.08 Hz frequentieband te behouden. Het voorafgaande gebruik van AFNI’s despike wordt verwacht om aliasing van resterende beweging artefacten in de naburige time-frames tijdens bandpass filtering69 te verzachten. Om de invloed van bewegingsartefacten verder te verzwakken, werden mogelijk met beweging verontreinigde tijdframes, gedefinieerd door een conservatieve FD > 0,15 mm drempel, uit de gegevens verwijderd (bekend als scrubbing van de gegevens)70. Deelnemers werden uitgesloten van verdere analyse als de gemiddelde FD groter was dan 0,15 mm, of als meer dan 30% van de frames werden gewist. Dit resulteerde in de uitsluiting van respectievelijk 4, 8 en 7 deelnemers in Steekproeven 1, 2 en 3 (Supplementary Table 2). Kwaliteitscontrole (registratie-check, tapijt-plots, zie bijvoorbeeld aanvullende Figs. 2-4) werd uitgevoerd gedurende de workflow.

Functionele connectiviteit analyse

De 122-parcel versie van de MIST71 multi-resolutie functionele hersenen atlas en grijze-stof maskers verkregen uit de anatomische beeld werden getransformeerd naar de native functionele ruimte. Deze atlas (geconstrueerd met behulp van de BASC-methode, dat wil zeggen bootstrap analyse van stabiele clusters) werd onlangs aangetoond dat goed presteren in connectiviteit gebaseerde voorspellende modellering72. Native-space atlas regio’s werden gemaskeerd met de grijze-materie maskers die werden verkregen uit de anatomische afbeelding en getransformeerd naar functionele ruimte eerder. Met deze atlas-individualisatie techniek, zal de uiteindelijke regionale signaal afkomstig zijn – met een hoge waarschijnlijkheid – van grijze-stof voxels voor elk onderwerp (die we zorgvuldig gecontroleerd handmatig voor alle onderwerpen), terwijl met de conventionele methode, een variabele verhouding van grijze en witte-stof voxels zijn opgenomen voor elk onderwerp. Daarom is de input van informatie van het weefsel-segmentatie proces zal naar verwachting subject-to-subject variabiliteit te verminderen (zie Supplementary Fig. 5 voor voorbeelden). Voxel-tijdenreeksen werden gemiddeld over deze geïndividualiseerde MIST regio’s en, samen met het gemiddelde grijze-stof signaal, behouden voor grafiek-gebaseerde connectiviteitsanalyse.

Regionale tijdreeksen werden geordend in grootschalige functionele modules (gedefinieerd door de 7-parcel MIST atlas) voor visualisatie doeleinden (Fig. 1). Partiële correlatie werd berekend over alle paren van regio’s (en globale grijze stof), zoals geïmplementeerd in de nilearn73 python module. Gedeeltelijke, in plaats van eenvoudige correlaties werden gebruikt om indirecte connectiviteit uit te sluiten74. Onze grafiek-modellering aanpak zorgde ervoor dat de globale grijze stof signaal wordt behandeld als een confound tijdens de berekening van de partiële correlatie coëfficiënten, maar, op hetzelfde moment, ook beschouwd als een signaal van belang, zoals het kan vertegenwoordigen waakzaamheid gerelateerde processen75. Partiële correlatiecoëfficiënten werden geordend tot 123 bij 123 (122 regio’s + globale grijze-stofsignaal) symmetrische connectiviteitsmatrices. De bovenste driehoek van deze matrices werd gebruikt als de feature-ruimte voor machine-learning-gebaseerde voorspellende modellering.

Predictieve model training en validatie

Whole-brain resting-state functionele connectiviteitsgegevens van studie 1 (N1 = 35, na alle uitsluitingen, zoals in ref. 8, Supplementary Table 2) werden gebruikt als de input feature-ruimte (P = 7503 features per deelnemer) om individuele pijngevoeligheidsscores te voorspellen, wat leidt tot een grote P-small N setting.

We construeerden een machine-learning pijplijn (https://github.com/spisakt/RPN-signature/blob/master/PAINTeR/model.py) in scikit-learn76, bestaande uit robust feature scaling (verwijdert de mediaan en schaalt met data quantiles), pre-selectie van features77, het selecteren van de K beste features met de sterkste relaties met de doelvariabele en een Elastic Net regressiemodel78 (een lineair model met gecombineerde L1 en L2-normen als regulariser). Het gebruik van elastisch net was een beslissing die werd genomen vóór de analyse. Onze belangrijkste motivatie om voor elastisch net te kiezen was dat het de mogelijkheid biedt om de sparsity (L1 vs. L2 regularisatie) als hyperparameter te optimaliseren, zodat we geen a-priori veronderstellingen hoefden te maken over de sparsity van de discriminatieve ground truth (zie ref. 79 voor de rationale). Samengevat waren de vrije hyperparameters van de machineleerpijplijn het aantal voorgeselecteerde kenmerken (K), de verhouding van de L1/L2-regularisatie en het gewicht (alpha) van de regularisatie. De hyperparameters werden geoptimaliseerd met een rasterzoekprocedure en een negatieve gemiddelde kwadratische fout als kostenfunctie. De waarden voor K varieerden van 10 tot 200 met stappen van 5, en voor alpha was de L1/L2-verhouding inbegrepen. De hyperparameteroptimalisering werd uitgevoerd in een leave-one-participant-out kruisvalidatie (interne valideringsfase). De kruisvalidatie omvatte de volledige pijplijn voor machinaal leren om te voorkomen dat er afhankelijkheden zouden ontstaan tussen de training- en de testmonsters. Merk op dat de fMRI voorbewerking onafhankelijk was tussen proefpersonen, dus niet opgenomen in de kruisvalidatie. De optimale hyperparameters bleken K = 25, L1/L2-ratio = 0.999 en alpha = 0.005.

Externe validatie werd uitgevoerd door toepassing van de RPN-handtekening op de fMRI-gegevens van studies 2 en 3 (N2 = 37, N3 = 19, na uitsluitingen, aanvullende tabel 2), eenvoudig door toepassing van de kenmerktransformatie (schaling) verkregen op monster 1 en vervolgens berekening van het dot-product tussen de individuele connectiviteit matrices en de niet-nul kenmerk gewichten verkregen in monster 1. De resulterende voorspellingen werden vergeleken met de voorspellingen van de RPN-handtekening. De resulterende voorspellingen werden vergeleken met de waargenomen QST-gebaseerde pijngevoeligheidsscores door het berekenen van de gemiddelde absolute fout (MAE), de gemiddelde gekwadrateerde fout (MSE) en de verklaarde variantie. Permutatie-gebaseerde p-waarden werden verkregen voor alle drie maten, met behulp van het mlxtend python pakket. Bovendien werd bootstrapping met conditionele dekking80 gebruikt om p-waarden te verkrijgen voor voorspellende connectiviteitsgewichten om de interpretatie te helpen. We construeerden 10000 bootstrap samples (met vervanging), met een grootte gelijk aan de originele sample, bestaande uit gepaarde hersen- en uitkomstdata. Het voorspellende model met de optimale hyperparameters werd aan elk monster aangepast. Ongecorrigeerde P-waarden werden berekend voor elke geselecteerde verbinding op basis van het aandeel van de gewichten onder of boven nul, zoals in bv. ref. 30. Merk op dat de interpretatie van deze p-waarden en betrouwbaarheidsintervallen (supplementaire tabel 4) beperkt blijft, aangezien zij afhankelijk zijn van de kenmerk-selectieprocedure.

Visualisatie van het voorspellende netwerk

De voorspellende interregionale verbindingen gemarkeerd door de niet-nul regressiecoëfficiënten van de RPN-signatuur werden weergegeven als een lint plot met behulp van het R-pakket circlize (Fig. 3). Overeenkomstige geïndividualiseerde hersengebied maskers werden terug getransformeerd naar standaard ruimte om een studie-specifieke regionale waarschijnlijkheid kaart (als gevolg van co-registratie nauwkeurigheid en individuele variabiliteit in morfologie) te creëren. Waarschijnlijkheid kaarten werden vermenigvuldigd met de som van de overeenkomstige regressiecoëfficiënten tot een regionale voorspellende sterkte kaart, die vervolgens werd gevisualiseerd met FSLeyes en MRIcroGL.

(https://www.mccauslandcenter.sc.edu/mricrogl) (Fig. 3) te creëren. De analyse van grootschalige resting-state netwerk-involvement (zoals gedefinieerd door de MIST71 hersenen atlas) werd uitgevoerd door het samenvatten en Z-transformeren van de voxel waarden over de zeven regio’s van belang. Polaire plot werd gemaakt met het R-pakket ggplot2.

Beschikbaarheid van software

De RPN-signatuur scores kunnen worden berekend op basis van structurele en resting-state functionele datasets door de software tool met dezelfde naam. De RPN-signature software tool bestaat uit de beschreven MRI processing pipeline en het functionele connectoom-gebaseerde predictieve model. Het is als broncode beschikbaar op https://github.com/spisakt/RPN-signature. Aangezien de software de Brain Imaging Data Structure (BIDS)82 en de BIDS-App specificatie volgt, biedt het een standaard command line interface en vertrouwt het op Docker-technologie. De docker image is gedeponeerd op Docker Hub: (https://cloud.docker.com/repository/docker/tspisak/rpn-signature) en is niet afhankelijk van enige software buiten het container image. Dit, samen met de volledig transparante continue integratie-gebaseerde ontwikkeling en geautomatiseerde tagging en versiebeheer, verbetert de beschikbaarheid van software en ondersteunt de reproduceerbaarheid van de resultaten van RPN-handtekeningen.

Rapportagesamenvatting

Verder informatie over de onderzoeksopzet is beschikbaar in de Nature Research Reporting Summary die aan dit artikel is gekoppeld.