痛みのない安静時の脳機能的結合が個人の痛み感受性を予測する

General considerations

研究デザインは、神経画像バイオマーカー50（ニューロンマーカー）に関する最近の勧告、要件、標準を慎重に考慮しながら、以下の考えに基づいて設定されました。

Maximise predictive performance

十分な効果量は臨床的有用性の基本要件であるため、ニューロマーカーの最適な感度を確保するために標準化された前処理パイプラインを採用した50。 fMRIの時系列データを抽出する際には、個々の解剖学的構造を取り入れた高精度の画像アライメントを使用した。 さらに、アーチファクト低減に関する最近の勧告とプロトコル51を採用し、コネクトームベースの解析の特別なニーズを満たすためにワークフローを最適化しました。 これは、既存のニューロイメージングソフトウェアに統一的なインターフェースを提供するコミュニティベースのPythonプロジェクトであるnipype52に基づき、一部、C-PAC53とniworkflows54のオープンソースプロジェクトのコードを再利用したもので、自社開発のオープンソースPythonソフトウェアライブラリPipelines Utilising a Modular Inventory（PUMI、https://github.com/spisakt/PUMI）を使用しました。 予測モデリング（機械学習）アプローチは、安静時脳機能ネットワークによって提供される豊富なデータを利用し、潜在的にはfMRIのhyperacuity55を活用するために利用された。

Assessing predictive power under realistic conditions

We used a preregistered, external validation strategy, that strictly separate model training and performance assessment. モデルのトレーニングには、研究1のデータのみを使用した。 検証のために、異なる研究センター、異なる装置、異なる研究スタッフで、2つの独立したサブ研究（研究2および3）を実施した。 また、一般化可能性を確保するために、検証手順にも適度な異質性を持たせ、研究センター間の手順、装置、画像処理シーケンス、参加者と研究者のコミュニケーションの言語などに適度な異質性を持たせ、研究設定の自由な調整を行った。

予測が神経信号によって駆動され、疼痛感受性に特異的であることを確認する

提案された疼痛感受性のマーカーが、実際に疼痛感受性に関連する神経信号によって駆動されることを確認するために、予測されたスコアと、事前に定義された（そして事前に登録した）種々の交絡因子およびバリデータ変数の相関を評価する。

Ensure accessibility of results

We applied a comprehensive pre-registration and made the source code of the method open-source and free available for the community. さらに、プラットフォームに依存しない、使いやすいドッカーコンテナを提供し、研究製品50として我々の予測モデルを使用する機会を提供し、任意の適切な画像データセットからすぐに痛みの感度予測を得ることができる。 参加者の年齢と性別は補足表1に示した。 研究1にはN1=39人の参加者がいた（参考文献8と同じサンプル）。 これはMZとTSWによってRuhr University Bochum（ドイツ）で行われ、機械学習による疼痛感受性予測のためのトレーニングサンプルとして、さらに予測の内部検証のための基礎として使用された。 研究2および3（N2＝48、N3＝29）は、それぞれFSとTSがUniversity Hospital Essen（ドイツ）で、BKとTKがUniversity of Szeged（ハンガリー）で実施し、外部検証用サンプルとして使用した。 参加基準と除外基準は3施設ともほぼ同じであり、表3に示すとおりである。 7201>

金属インプラント、取り外せないピアス、ピースメーカー、頭・首の位置にあるタトゥー、妊娠または既知の閉所恐怖症はMR測定の禁忌と見なされた。 参加者は実験の2時間前（研究3を除く）にカフェインの摂取を控え、実験当日と前日にアルコールの摂取を控えるよう求められた。

本研究はヘルシンキ宣言に基づき実施され、ヒトを対象とした作業に関するすべての関連倫理規定を遵守し、地方または国の倫理委員会（登録番号: No: ルール大学ボーフム、大学病院エッセン、ETT TUKEBハンガリーでそれぞれ4974-14、18-8020-BO、057617/2015/OTIG)。

画像検査と定量的感覚検査（QST）は、研究1では同日に、研究2および3では平均2-3日間隔で行われた（詳細は補足表1参照）。

Measures-functional MRI

高解像度の解剖学的および開眼安静時のfMRI測定が，すべての参加者から行われた。 スキャンパラメータ（装置を含む）は各施設で異なり、表4に示すとおりである。 測定中、参加者は眠らずにじっとリラックスして横になっているように指示され、一切の動きを避けた。 頭部の動きを制限するために、発泡パッドと、研究1と2では空気枕が使用された。

Measures-QST

QST protocol28に従って熱（HPT）、冷（CPT）、機械（MPT）の痛み閾値が取得された。 温感（WDT）、冷感（CDT）、試験2および試験3では機械的（MDT）検出閾値は、追加の対照尺度として取得された。 すべての感覚測定は、左前腕掌部、手首の稜線近傍で行われた。 QSTの枠組みでは、熱的閾値は限界の方法を用いて決定される。 このため、研究1ではMSA熱刺激装置（Somedic, Hörby, Sweden）、研究2および3ではPathway熱刺激装置（Medoc Ltd., Ramat Yishai, Israel）を用いて皮膚に温度の上昇と下降を適用しました。 すべての試験で、ATSサーモードを30×30mmの皮膚表面に使用し、ベースライン温度を32 °Cに設定した。 参加者は、痛みの発生をボタン操作で示すよう指示された。 被験者間の分散を減らすため、すべての熱閾値6について、3ではなく（元のプロトコルと同様に）28、刺激の繰り返しが行われた。 さらに、最初の測定は試験刺激として解析から切り離された。 HPT と CPT は残りの 5 つの閾値温度の算術平均値として算出した。 MPTとMDTは階段法を用いて決定した。 左前腕掌部にピンクリック刺激（MRC Systems, Heidelberg, Germany）を5回増加、5回減少させながら交互に与え、参加者はその刺激を不快か非快感かに分類するよう指示された。 機械的検出閾値はvon Freyフィラメント刺激と同様に評価した。 MPTとMDTは、5回の階段昇降閾値ランで決定された対数変換幾何平均力として計算された。

追加測定

年齢、性別、自己申告による身長、体重、女性参加者の場合は最終月経初日と避妊薬の使用について、すべての測定前に記録された。 さらに、研究1と2では、自己申告による1週間のアルコール摂取量と教育レベル（小学校、中学校、大学）が記録された。 QSTの前に、参加者は、疼痛感受性質問票（PSQ）56、疼痛破局感尺度（PCS）57、状態像不安尺度（STAI）58、ドイツ版うつ病尺度（ADS-K、疫学研究センター）59を記入し、さらに、第2、3研究ではピッツバーグ睡眠の質指数（PSQI）60と知覚ストレス質問票（PSQ20）61も記入した。 研究2および3では、MRIおよびQST測定前に血圧が測定された。 さらに、試料1では、fMRI検査の前日に行われた並行実験からT50値を入手することができた。 T50は、50の熱痛評価（0、無痛から100の耐え難い痛みまでのスケール）を誘発するのに必要な温度（℃）を表している。 T50値は、42.5℃から48℃の間で15回の緊張性熱痛刺激（持続時間：16秒）に反応して得られた評価の非線形（2次多項式）補間から得られ、疑似無作為化グリッドサーチ方式で提示された。

研究1では、HPT、CPT、MPTをZ変換（平均中心化、標準化）し、HPTとMPTを反転（-1倍）させ、Z値が高いほど疼痛感度が高いことを示した。 そして、Z変換した変数の算術平均を各参加者について計算し、疼痛感受性スコアと定義した。 研究2および3では、研究間で同じ尺度が使用されるように、研究1の母平均と標準偏差に基づいてZ変換した以外は、同じ手順が適用された。 極端なQST値は、文献28で報告されている標準的な95％百分位を使用して定義した。3つのモダリティのうち少なくとも2つで極端なHPT、CPTまたはMPT値を示す参加者は除外された。 このスクリーニングの結果、サンプル1、2、3ではそれぞれ0、3、2人の参加者が除外された（補足表2）。

fMRI前処理

fMRIベースの機能的結合はスキャナ内運動アーティファクト62、63に影響を受けやすいため、結合に基づく予測には適切な前処理と信号洗浄が鍵である。 安静時機能的MRIデータは、3つの研究すべてで同じように前処理された。 ニタイプベースのワークフローは、補足図1に示されている。 また、このワークフローは、サードパーティのニューロイメージングソフトウェア、C-PAC53とniworkflows54のソフトウェアツールから採用したコード、および社内のパイソンルーチンを利用した。 解剖学的画像はANTs65（ソースコードはhttps://gist.github.com/spisakt/0caa7ec4bc18d3ed736d3a4e49da7415参照）を用いて1mm解像度のMNI152標準脳テンプレートに線形および非線形に共置換した。

機能画像はFSL flirtの境界ベース登録技術で解剖学的画像に共置換した。 得られた変換はすべて今後の使用のために保存された。 機能画像の前処理は、リサンプリングを行わず、ネイティブな画像空間で行われた。 FSL mcflirt を用いて、リアライメントに基づく動き補正を行った。 その結果、6つの頭部運動推定値（3つの回転、3つの並進）、それらの2乗バージョン、それらの導関数、2乗導関数（Friston-24-expansion66として知られている）が計算され、厄介な補正のために保存された。 さらに、頭部運動はPowerの方法63に従ってフレーム単位の変位（FD）時系列としてまとめられ、データの打ち切りおよび除外に使用された。 動き補正後、AFNI despike67を用いて時系列データ中の外れ値（モーションスパイクなど）を減衰させた。 7201>

ノイズ回帰ステップでは、6つのCompCorパラメータ（ノイズ領域時系列の6つの第一主成分）、Friston-24運動パラメータおよび線形トレンドが一般線形モデルで時系列データから除去された。 残差データに対して、AFNIの3DBandpassで時間帯域フィルタリングを行い、0.008〜0.08Hzの周波数帯域を保持した。 AFNIのデスパイクを事前に使用することで、バンドパスフィルタリング中に隣接する時間フレームへの残留モーションアーチファクトのエイリアシングを減衰させることが期待される69。 モーションアーチファクトの影響をさらに減衰させるために、保守的なFD > 0.15 mmの閾値で定義された潜在的にモーション汚染された時間フレームをデータから削除した（データのスクラビングとして知られている）70。 平均FDが0.15mmを超えた場合、または30%以上のフレームがスクラビングされた場合、参加者はさらなる分析から除外された。 その結果、試料 1、2、3 において、それぞれ 4 名、8 名、7 名の参加者が除外された（補足表 2）。

機能的結合解析

MIST71 multi-resolution functional brain atlasの122区画版と解剖学的画像から得られた灰白質マスクは、ネイティブな機能空間へ変換された。 このアトラス（BASC法、すなわち安定したクラスタのブートストラップ分析によって構築）は、最近、結合性に基づく予測モデリングで良好に機能することが示された72。 ネイティブ空間アトラス領域は、解剖学的画像から得られたグレイマターマスクでマスクされ、あらかじめ機能空間へ変換されたものである。 このアトラス個別化手法では、最終的な領域信号は、各被験者のグレイマターボクセルから高い確率で発生する（これは全ての被験者について手動で慎重に確認した）。一方、従来の方法では、被験者ごとにグレイマターボクセルとホワイトマターボクセルの比率が変動していた。 したがって、組織セグメンテーションの情報を入力することで、被験者間のばらつきを減らすことが期待できる（例として補足図5を参照）。 ボクセル時系列は、これらの個別化されたMIST領域上で平均化され、平均グレイマター信号とともに、グラフベースの結合性解析のために保持された。

領域時系列は、可視化の目的で大規模機能モジュール（7区画MISTアトラスによって定義）に並べられた（Fig.1）。 偏相関は、パイソンモジュールnilearn73で実装されたように、領域（およびグローバル灰白質）のすべてのペアで計算された。 間接的な結合を除外するために、単純な相関ではなく、偏相関を使用した74。 我々のグラフモデリングアプローチは、偏相関係数を計算する際に、グローバル灰白質信号が交絡として扱われることを保証するが、同時に、それは警戒関連のプロセスを表すかもしれないので、関心のある信号として考慮される75。 偏相関係数は123×123（122領域＋グローバルグレイマター信号）の対称的な連結行列に整理された。 これらの行列の上三角形は、機械学習ベースの予測モデリングのための特徴空間として使用された。

予測モデルの学習と検証

研究1の全脳安静時機能結合データ（N1 = 35、すべての除外後、文献8、補足表2の通り）は、個々の疼痛感受性スコアを予測するための入力特徴空間として使用（被験者あたりP = 7503特徴）、大きなP-小さなN設定となった。

scikit-learn76で機械学習パイプライン（https://github.com/spisakt/RPN-signature/blob/master/PAINTeR/model.py）を構築し、ロバスト特徴スケーリング（中央値を削除し、データ分位でスケーリング）、特徴の事前選択77、ターゲット変数との関係が最も強いK個の最適特徴の選択、Elastic Net回帰モデル78（正規化としてL1およびL2ノームを結合した線形モデル）から構成されていた。 Elastic Netの使用は、分析に先立って決定されたものである。 elastic netを選んだ主な動機は、ハイパーパラメータとしてスパース性（L1対L2正則化）を最適化できるため、判別可能なグランドトゥルースのスパース性について事前に仮定する必要がないことです（根拠は文献79をご参照ください）。 要約すると、機械学習パイプラインの自由なハイパーパラメータは、事前に選択された特徴の数（K）、L1/L2正則化の比率、正則化の重み（α）であった。 ハイパーパラメータはグリッドサーチ法により最適化され、コスト関数として負の平均2乗誤差が用いられた。 Kの値は10〜200で5刻み、αはL1/L2比を含む。 ハイパーパラメータの最適化は、一人抜きクロスバリデーション（内部検証段階）で行われた。 クロスバリデーションは、トレーニングサンプルとテストサンプル間の依存性を避けるために、完全な機械学習パイプラインを組み込んだ。 なお、fMRIの前処理は被験者間で独立していたため、クロスバリデーションには含まれなかった。 最適なハイパーパラメータは、K = 25, L1/L2-ratio = 0.999, alpha = 0.005 であった。

外部検証は、研究2および3のfMRIデータ（N2 = 37, N3 = 19、除外後、補足表2）にRPN署名を適用して、サンプル1で得られた特徴変換（スケーリング）、そして個々の結合度行列とサンプル1で得られたゼロ以外の特徴重み間のドットプロダクトを簡単に算出し実施した。 得られた予測値は、平均絶対誤差（MAE）、平均二乗誤差（MSE）、説明される分散を計算することにより、観察されたQSTに基づく疼痛感受性スコアと比較された。 並べ替えに基づくp値は、mlxtend pythonパッケージを使用して、すべての3つの測定値について得られた。 さらに、条件付きカバレッジ80によるブートストラップは、解釈を助けるために、予測連結性重みのp値を提供するために使用された。 我々は、元のサンプルと同じ大きさの、対になった脳と結果のデータからなる10000個のブートストラップサンプル（置換あり）を構築した。 最適なハイパーパラメータを持つ予測モデルを各標本に当てはめた． 未補正のP値は、例えば文献のように、0以下または0以上の重みの割合に基づいて、選択された各接続について計算された。 30. なお、これらのP値および信頼区間（補足表4）は、特徴選択手順に条件付けされているため、その解釈は限定されたままである。

交絡因子分析

潜在的な交絡因子を探索するために、予測された疼痛感受性スコア（またはサンプル1の場合には交差検証された予測値）は、平均および中央値FD、スクラビング体積の割合、MRIおよびQST測定前の収縮期および拡張期血圧（血圧は機械的疼痛に対する感受性と関連すると以前に報告81されているので）と対照された。 MRI検査とQST検査の間の時間的な遅れ（予測の時間的安定性をテストするため）、年齢、性別、BMI、最終月経初日からの日数、アルコール摂取量（単位/週）、教育レベル、状態不安と特性不安（STAI）。 うつ症状のスコア（ADS-K）、自己報告による疼痛感受性（PSQ）および疼痛異化（PCS）、知覚ストレス（PSQ20）、睡眠の質（PSQI）、非痛覚性QST検出閾値（CDT、WDT、MDT、利用可能な場合）。 さらに、研究1では、予測値を痛み処理脳領域のT50値、MRスペクトロスコピーに基づくGABAおよびグルタミン酸／グルタミンレベルと比較した（詳細は文献8を参照）。

Visualization of the predictive network

RPNシグネチャの非ゼロ回帰係数によって強調された予測領域間結合は、Rパッケージのcirclizeを用いてリボンプロットとして表示した(Fig. 3)。 また、個人別に作成した脳領域マスクを標準空間に逆変換し、研究別の領域確率マップを作成した（co-registrationの精度と形態の個人差を反映した）。 確率マップに対応する回帰係数の和を乗じて領域予測力マップを作成し、FSLeyesとMRIcroGL.

(https://www.mccauslandcenter.sc.edu/mricrogl) で可視化した（Fig.3）。 MIST71脳アトラスで定義された大規模な安静時ネットワーク関与の解析は、7つの関心領域にわたるボクセル値を要約し、Z変換することで行われた。 Rパッケージのggplot2を用いて極座標プロットを行った。

Software availability

RPN-signature scoreは、同名のソフトウェアツールによって、構造および静止状態の機能データセットに基づいて計算することができる。 RPN-signatureソフトウェアツールは、説明したMRI処理パイプラインと機能的なコネクトームベースの予測モデルから構成される。 ソースコードはhttps://github.com/spisakt/RPN-signatureで公開されています。 このソフトウェアはBrain Imaging Data Structure (BIDS)82 とBIDS-Appの仕様に準拠しているため、標準的なコマンドラインインターフェイスを提供し、Docker技術に依存している。 Dockerイメージは、Docker Hubに寄託されている。 (https://cloud.docker.com/repository/docker/tspisak/rpn-signature) に寄託され、コンテナイメージ以外のソフトウェアには一切依存しない。 これは、完全に透明な継続的インテグレーションに基づく開発、自動化されたタグ付けとバージョン管理とともに、ソフトウェアの可用性を高め、RPN署名の結果の再現性をサポートします。

Reporting summary

研究デザインに関するさらなる情報は、この記事にリンクされている Nature Research Reporting Summary で入手可能です。