生体信号の移動エントロピー解析に基づく自律神経機能の評価に関する研究

生体信号の移動エントロピー解析に基づく自律神経

機能の評価に関する研究

著者

坂本 真

学位授与機関

Tohoku University

修士学位論文

生体信号の移動エントロピー解析に基づく

自律神経機能の評価に関する研究

東北大学大学院医工学研究科

医工学専攻

坂本

真

平成 24 年 2 月 9 日

目 次

第1章 緒言 1 1.1 研究の背景と目的 . . . . 1 1.2 本研究のアプローチ . . . . 2 第2章 血圧調整システムとその評価方法 5 2.1 生体循環制御系と血圧調整システムのメカニズム . . . . 5 2.1.1 生体循環制御系 . . . . 5 2.1.2 血圧調整システムのメカニズム . . . . 7 2.2 血圧調整システムの評価方法 . . . . 10 2.2.1 従来の評価方法 . . . . 10 2.2.2 従来の評価方法の問題点 . . . . 12 第3章 システム内の信号間における因果関係の推定 13 3.1 従来の2つの信号間における因果関係の推定方法 . . . . 13 3.2 移動エントロピーの概要 . . . . 14 3.3 結合写像格子系を用いた移動エントロピーのシミュレーション . . . . 18 第4章 生体循環系モデルにおけるシミュレーション 22 4.1 シミュレーションに用いた生体循環系モデルの概要 . . . . 22 4.2 シミュレーションの概要 . . . . 26 4.3 シミュレーションの結果 . . . . 31 第5章 実際の生体データを用いた血圧調整システム内の信号間における因果関係の推定 33 5.1 血管系の信号 . . . . 33 5.1.1 脈波の振幅 . . . . 33 5.1.2 脈波伝播時間 . . . . 34 5.1.3 血圧と脈波伝播時間の関係式 . . . . 35

5.1.4 血管コンプライアンスと総末梢血管抵抗 . . . . 36 5.2 実験の概要 . . . . 36 5.3 解析方法 . . . . 37 5.4 解析結果と考察 . . . . 42 5.4.1 血圧調整システム内の3つの信号間における因果関係 . . . . 42 5.4.2 MTEDの再現性の検証. . . . 52 第6章 結言 55 6.1 本研究の総括 . . . . 55 6.2 今後の課題 . . . . 56 参考文献 57 謝辞 61

第

1

_{章 緒言}

1.1

研究の背景と目的

生体循環系は，非常に多くの要素間の相互作用により成り立つ非常に複雑なシステムである． この循環系内には血液を循環させるポンプ役の心臓と，送り出された血液を体内の各所に運搬 するための血管，および流れる血液そのものの力学的特性を基本として，それらを適切に制御 する体液性，神経性など複数のフィードバックループを有する[1]．複雑なフィードバックルー プは一般に非線形・非定常性を含み，生体におけるゆらぎやカオス的な振る舞いの原因となっ ていると考えられている．このため，循環系を厳密に解析することは非常に困難であり，完全 な定量的モデルは構築されていない．循環系を血行力学的特性と，中枢神経系などが担うそれ を制御する部分との大きく２つに分けて考えると，前者は制御対象，後者は制御器とみなすこ とができる．これまでに，計測の容易さなど臨床的な立場から，前者の血行力学的特性をモデ ル化する試みが多くなされてきた．これらのモデルは簡単のため線形近似などを用いており， 安静状態の循環系ダイナミクスの解析や人工心臓制御などにも応用されている[2]． 一方，生体循環系は生命維持に関わる重大な任務を負った制御システムであることから，循 環系に関わる疾患には重大な疾患となるものが多数ある．その例として，日本人の三大死因で もある脳血管障害（脳卒中）や心臓疾患が挙げられる．これらの疾患は高血圧症と関係が深い とされており，図1.1に示したように，収縮期血圧が高くなるほど脳卒中や虚血性心疾患によ る死亡リスクが上昇するという報告がされている[3]．また，Takagiら[4]は18年間の追跡調 査（1819例）を行い，図1.2に示したように，生存率は正常血圧者，境界域高血圧患者，高血 圧患者の順に低下していたことを報告している．なお，高血圧症は世界保健機構（WHO）と国 際高血圧学会（ISH）が共同で発表したWHO/ISHガイドラインの基準によれば，拡張期血圧 が90mmHg以上，または収縮期血圧が140mmHg以上の症状であり，その中でも拡張期血圧 が94mmHg以下かつ収縮期血圧が149mmHg以下となる症状を境界域高血圧としている（図 1.3参照）．これらの報告から分かるように，高血圧症と死亡率の高さとの関連から，高血圧症 の原因の解明が喫緊の課題となっている． 高血圧症の原因としては大動脈弁閉鎖不全症（心拍出量に関するもの）や粥状動脈硬化（循

図 1.1 脳卒中や虚血性心疾患による死亡リスクに対する収縮期血圧と年齢の影響（ [3]より 引用）．収縮期血圧が高くなるほど脳卒中や虚血性心疾患による死亡リスクが上昇している． 環血液量に関するもの）が挙げられるが，このように原因がわかるものは全体の5∼10%に過ぎ ず，90%以上のものについては原因が明確にはなっておらず，本態性高血圧と呼ばれている[6]． 本態性高血圧の原因については様々な検討がされ，その1つとして，圧受容器反射の異常が報 告されている [7, 8]．圧受容器反射とは，頸動脈洞と大動脈弓にある圧受容器により血圧を常 にモニターして中枢神経に伝え，心臓や血管の働きを介して血圧の変動を急速に修正する循環 反射であり，血圧調整システムの大部分を担っている．本態性高血圧患者では，この圧受容器 反射がうまく機能していないか，その反射作用が弱いために血圧調整が正常に行われない可能 性がある．したがって，この血圧調整システムの機序を詳しく解析することによって，高血圧 症の原因の解明につながる知見が得られる可能性がある．

1.2

本研究のアプローチ

血圧調整システムを評価する場合，従来では，血圧調整の大部分を担う圧受容器反射機能を 評価することが多かった．第2章で述べるように，圧受容器反射は，心臓血管系を介して血圧 を制御する反射であり，血圧から血圧への閉ループ系である．しかし，この閉ループ系では， 各信号に生じた変動の間の因果関係が複雑であり，閉ループ特性を推定することが困難である

図 1.2 収縮期血圧／拡張期血圧分類別の累積生存率の推移（[5]より引用）． 図 1.3 WHO/ISHガイドラインによる血圧の分類． とされてきた．したがって，閉ループの影響を考慮せず，圧受容器反射を血圧変動から心拍数 変動への開ループ系とみなし，線形性に基づき ∆心拍数/∆血圧 などの指標を用いたりして 圧受容器反射機能の評価が行われていた．だが，このような従来の評価方法では，非線形・非 定常性を含む血圧調整システムの評価が精度良く行えない可能性がある．また，前述したよう に，従来の血圧調整に関する研究では，ほとんどが心拍数変動との関係を調査したものであっ たが，実際には末梢血管抵抗などの心拍数以外の生理信号も大きく関わっていることが知られ ている．したがって血圧調整システムの機序をより詳しく解明するには，血圧変動と複数の生 理信号間との関係も評価する必要がある． そこで本研究では，2つの信号間における確率的因果関係を解析することができる移動エン トロピー（tranfer entropy: TE）解析に着目した．TEは確率変数によって定義される情報量 を用いて計算され，線形解析だけでなく，非線形やカオス性を含む時系列の解析にも適してい る．加えて，TEを用いた解析で遅延時間を考慮することにより，時系列相互間の情報の移動

方向，すなわち因果関係の方向を推定することができる．したがって，本研究では，血圧変動 と多くの要素から構成される血圧調整システム内の複数の信号間においてTE解析を行い，そ れぞれの信号間の因果関係を調べることで，血圧調整システムの機序の解明を目指した． 以下，本論文の構成を示す． 第2章では，心臓・血管系における血圧調整システムのメカニズムについて述べるとともに， その評価方法に関するこれまでの研究について述べる． 第3章では，従来の2つの信号間における因果関係を推定する方法について述べるとともに， 本研究で血圧調整システムを評価するために着目した移動エントロピーについて，その概要と 計算方法について述べる． 第4章では，生体循環系モデルを用いたシミュレーションを行い，移動エントロピー解析に よる血圧調整システムの評価の有効性について述べる． 第5章では，実際の生体データを用いて，血圧調整システム内の信号間における因果関係の 推定を行った結果について述べ，考察を行う． 第6章は結言であり，本論文の総括を行い，今後に残された課題について述べる．

第

2

_{章 血圧調整システムとその評価方法}

2.1

生体循環制御系と血圧調整システムのメカニズム

2.1.1

生体循環制御系

生体循環系は，人体に不可欠な血液の運搬を任務とした大切なシステムである．この循環器 系は，血液を身体のすみずみまで運搬するためのポンプである心臓と，ポンプによって送り出 された血液の通路である血管から成っており，心臓を中心とした複雑な閉鎖系をなしている． 心臓は，洞房結節で発生する約1Hzのパルスによって周期的に収縮し，全身に血液を循環させ ている．この収縮の1分間あたりの回数を心拍数（heart rate: HR），もしくは心収縮率と呼 び，心臓から送り出される血液の量を心拍出量（cardiac output: CO）と呼ぶ．心臓から拍出 された血液のエネルギーは動脈壁に位置エネルギーとして蓄えられ血圧（blood pressure: BP） を発生させる．血圧は，心臓の収縮期に最高値へ達し，拡張期に最低値へ至るが，これらはそ れぞれ，収縮期血圧（systric blood pressure: SBP）と拡張期血圧（diastolic blood pressure: DBP）と呼ばれ，両者の差（=SBP–DBP）を脈圧（pulse pressure: PP）と呼んでいる．この 血圧の情報は圧受容体によって検出され，心臓血管中枢へ伝えられる．心臓血管中枢では，こ の情報をもとに適当な操作量を決定し，自律神経系の活動を介して心臓の制御を行い，また血 管平滑筋に伝わることにより全末梢血管抵抗（total peripheral vascular resistance: TPR）を 変化させている． 生体にとって最も基本的な機能である循環，呼吸，消化，代謝，分泌，体温維持，排泄，生 殖などは自律機能と呼ばれ，自律神経が平滑筋，心筋，腺を支配することで，調節が行われて いる．体性神経系が随意的な制御を受けるのに対し，自律神経系は随意的な制御を受けないた め，植物神経系や不随意神経系とも呼ばれる．図2.1に示したように自律神経系の遠心路は胸 腰髄から起始する交感神経系と脳幹および仙髄から起始する副交感神経系に大別される[9]． 内臓器官の多くは交感神経と副交感神経の二重支配を受けており，心臓も自動能が変調をう ける形で2つの神経に制御されている．心臓を支配する自律神経は心臓交感神経，心臓副交感 神経と呼ばれ，それぞれ心臓の興奮や収縮に対して促進性，抑制性に働く．つまり，心臓交感 神経の活動によってノルアドレナリンが神経伝達物質として分泌され，心臓側のβアドレナリ

ン作動性受容体に作用して心収縮率を上昇させたり，心筋収縮性や弛緩速度を増加させる．一 方，心臓迷走神経が活動する場合は，アセチルコリンが分泌され，コリン作動性受容体に作用 して心収縮率を減少させる． 血管も心臓と同様に自律神経の支配を受けている．血管を構成する平滑筋は血管運動神経支 配や液性調整を受け，ある程度の緊張状態（基礎緊張，basal tone）を保っている．血管運動 神経には交感神経性血管収縮繊維，交感神経性血管拡張繊維，副交感神経性血管拡張繊維があ る．このうち，前一者は最も広範囲に分布しており，常に調節活動を行っている． この様に，自律神経系による調節は交感神経系と副交感神経系のバランスによって成り立っ ているが，このバランスが崩れると自律神経失調状態となる．20世紀初頭にEppingerとHess は，自律神経失調症を交感神経緊張症と副交感神経緊張症とに分類した[10]．前者では高血圧， 頻脈，瞳孔散大，便秘，一方，後者では低血圧，胃酸過多，下痢，喘息など失調状態になる部 位によって異なった症状が現れる．

(a) Sympathetic Nervous (b) Parasympathetic Nervous

(a) Sympathetic Nervous (b) Parasympathetic Nervous

2.1.2

血圧調整システムのメカニズム

生体では，循環系の物理的平衡により血圧を一定に保つよう調整されたり，圧受容器反射な どの反射性調節や体液量調節により血圧の調整が行われたりする． 圧受容器反射は受容器（動脈圧受容器）により血圧を常にモニターして中枢神経系に伝え， 血圧の変動を秒や分のオーダーで急速に修正する循環反射である[12]．動脈圧受容器は頸動脈 洞と大動脈弓にあり，それぞれ頸動脈洞圧受容器と大動脈圧受容器と呼ばれる．動脈圧受容体 は伸展受容器であり，動脈壁の中膜に近い外膜に分布する受容器が動脈圧により伸展されてイ ンパルスを発生する． 図2.2は，圧受容器反射による血圧調整システムのメカニズムについてまとめたものである． 血圧が上昇して圧受容器が刺激されると，心臓および血管を支配する交感神経の活動が抑制さ れ，相反的に心臓迷走神経（心臓副交感神経）の活動が促進される．同時に副腎髄質からのカ テコールアミン分泌が減少して，交感神経の抑制と同じ効果をもたらす．このため，心臓では 心拍数が下がり，心臓収縮性が低下するために心拍出量が減少する．血管では血管収縮繊維の 活動が抑制されるため，血管が拡張して総末梢血管抵抗が低下する．血圧は総末梢血管抵抗と 心拍出量の積であるため，血圧は下降し，元の値付近に戻ることになる．また，血管収縮繊維 は静脈系も支配し，血圧が上昇すると血管コンプライアンスは増加し平均体循環圧は減少する． これも，心拍出量を減少させ，血圧を下降させる．さらに，圧受容器反射はバゾプレッシン分 泌を減少させ，血管収縮を和らげることで血液量を減少させる．これも血圧の下降に貢献する． これらの働きのほとんどは固定プログラムによる反射によって制御されているが，血圧変動 がこのような反射システムの固定プログラムだけでは是正しきれないような場合には，さらに 高位の制御システムが参加して，反射機能は修正される．このように，心臓・血管による血圧 調整システムは多数のフィードバックによって複雑にコントロールされているが，さらにこの システムは呼吸性の胸郭内圧の外乱等にさらされており，複雑な様相を呈している．図2.3は， この心臓血管・呼吸系による制御メカニズムを記述した，Mulderによる循環系の血圧制御定 性モデルである[13]．

۝טଇঞ ढ़੤ڥը ¸ଇঞ ݗ࠰ଇঞ Ϝࢄຓઑ ۝טڈ 8K GK B EB 6 ۝ ט ڥըࢌࢹ३ډԩ௫ ৌ ฝ ճ ی ໙ ढ़ಈࢽ໷ ढ़ಈ ГӖಈࢽ໷ ڥϜ ؖ൉ఝϜ ܉Օ३ډԩ௫ ܉Օ३ډԩ௫ ൐܉Օ३ډԩ௫ ढ़੤ ੒งऊ୰ܨ ó ó ô ô 図 2.3 Mulderによる循環系の血圧制御定性モデル [13]（M:心筋細胞，PM:洞房結節，RV: 抵抗血管，CV:容量血管，A:動脈圧）．

2.2

血圧調整システムの評価方法

2.2.1

従来の評価方法

従来における血圧調整システムを評価する方法の1つとして，前節で述べた圧受容器反射機 能の評価が行われてきた．圧受容器反射機能は，血圧の変化に対するRR間隔（隣り合うR波 の間の時間間隔）の変化のしやすさとして定義される． 主な圧受容器反射機能の評価方法としては，次のような方法がある． a） Valsalva（息こらえ）法[14] 息をこらえることで胸腔内圧を変化させて血圧を動かし，それに対応する心拍反応の大 きさを観測する方法である． b） 血管収縮・拡張薬物注入法 [15] フェニレフリンやニトログリセンなどの血管収縮性および拡張性薬物を静注し，収縮期 血圧の上昇とそれに引き続くRR間隔の関係より1次回帰直線を求め，その傾きから評 価する方法である． c） neck chamber 法[16] 頸部に装着したネックチャンバーに圧力を加え，頸動脈を圧迫あるいは吸引して圧受容 器を刺激した際のRR間隔と平均血圧を求め，刺激直前の値からの差同士の比から評価 する方法である． しかし，これらの方法で圧受容器反射機能の評価を行う場合，侵襲的な試験を被験者に課す こととなり，被験者の負荷が大きくなる．そこで，非侵襲的に圧受容器反射機能を評価する方 法としてα-indexやρmaxが提案されている． d） α-index 圧受容器反射を血圧変動から心拍数変動への開ループ系とみなし，その伝達関数のゲイ ンを圧受容器反射の感度とするものである[17]．α-indexは低周波成分（0.04∼0.15Hz） や高周波成分（0.15∼0.45Hz）におけるコヒーレンスが高い周波数領域に限って計算さ れ，それぞれαLF，αHF と表される．しかし，呼吸の影響を避けるため，αLFに限って 計算されることが多い [18]．αLF は(2.1)式のように計算される．

αLF = √ SRRI SSBP (2.1) なお，SRRI，SSBP はそれぞれRR間隔，収縮期血圧のパワースペクトルにおける低周 波成分のパワーを表す． e） ρmax 著者らの研究グループにおいて，血圧と心拍数変動の関係に基づく評価指標としてρmax を提案している[19–21]．ρmaxは血圧と心拍数変動のMayer波帯域（0.05∼0.15Hz）に おける最大相互相関係数で表され，圧受容器の感受性が高いほど，値は1に近づく．他 の方法と異なり，2変数間の時間差を考慮しているのが特徴である．次にρmaxの解析方 法について詳しく述べる． 1）測定された心電図波形から RR間隔を算出し，その逆数に60 を乗算して心拍数 （HRV）を求める．連続血圧波形からは拍内平均血圧（M BP）を求める． 2）拍毎のデータであるHRV とM BP に3次のスプライン補間を施し，再サンプリン グする．再サンプリング周波数は2Hz以上とする． 3）両データのMayer波帯域成分をディジタルフィルタを用いて抽出する． 4）結果表示の便宜上，HRV に–1を乗算する． 5）処理後のHRV とM BP をそれぞれx(t), y(t)としたとき，ある時刻t[s]において， その前後60秒分，合計120秒間にわたりハミング窓をかけてから，(2.2)式のよう に計算される2乗平均値で規格化された相互相関関数（相互相関係数）ρxy(τ )を計 算する． ρxy(τ ) = ϕxy(τ ) √ ϕxx(0) ϕyy(0) (2.2) ここで，ϕxy（τ）はx(t)からy(t)への相互相関関数，ϕxx(0)，ϕyy(0)はそれぞれx(t)， y(t)の自己相関関数である． 6）以上の手順を経て算出されたρxy(τ )をもとに，最大相互相関係数ρmaxは， ρmax= max 0s<_{=τ <=10s}ρxy(τ ) (2.3) のように定義される．

以上の処理を，再サンプリング時間毎に施すことにより，ρmaxを経時的に求めることが できる．ここでHRに負号を付けて解析を行うのは，値が大きいほど感受性が高いと表 現するためである．このρmaxは，血圧変動から心拍変動への線形結合の強さを時間領域 で求め，最大値が1になるように規格化した指標であるといえ，0.1Hz付近のコヒーレ ンス関数の値とほぼ同等である．

2.2.2

従来の評価方法の問題点

2.2.1項で述べた圧受容器反射機能の評価方法では，実際には閉ループを成している血圧変 動と心拍数変動の間を開ループ系とみなしているため，正確な評価を行っているとは言えず， 基本的に信号間の因果関係を知ることはできない．さらに，血圧調整システムでは血管調節な どの心拍数以外の生理信号も存在しており，システムの機構を理解するためには血圧と心拍数 との関係だけでは不十分である可能性がある． そこで本研究では，統計的手法に基づき，2つの信号間における因果関係を解析することが できる手法を用いて圧受容器反射機能の評価を行うこととした．詳細については，第3章で述 べる．

第

3

_{章 システム内の信号間における因果関係の}

推定

本章では，血圧調整システム内の信号間における関係性のうち，とりわけ因果性を推定する ために本研究で採用した方法について述べる．また，生体循環系モデルのシミュレーションで 得られる各信号に対し採用した解析方法を適用した結果を示し，その有効性について考察を 行う．

3.1

従来の

2

つの信号間における因果関係の推定方法

従来，2つの信号間における関係性を調べる方法としては，相互相関関数，コヒーレンス関 数などの線形解析法や，相互情報量，相対エントロピーなどを用いる非線形解析法が挙げられ る．これらは，2つの信号間における関係の強さの程度を評価するのに適しているが，2つの 信号に対して対称な形をしているため，因果性という情報の流れる方向まで含めて解析するこ とには適していない．相互相関関数でも，情報の流れが一方向のみであれば位相差からその方 向を推定できるが，2つ以上の成分があり双方向に流れている場合には，その検出は困難であ る[22]． これに対し，因果関係の推定に適した方法としては，遅れ時間を考慮した相互情報量である遅 延相互情報量（time delayed mutual information: TDMI）[23]や移動エントロピー（transfer entropy: TE）[24]などがある．しかし，是永ら[25]はTDMIは周期的な変化をする信号間で は正しく因果関係が推定できないと報告している．本研究の解析対象である生体におけるシス テム内の信号は周期的変動をするため，TDMIによって生体システム内の信号間の因果関係を 推定することは難しいと考えられる．一方，TEは周期的変動をする信号間であっても正しく 因果関係を推定できることが分かっている [25]．また，Schreiber [24]はTEを用いることで 結合写像格子系を用いた情報移動方向の推定を行っており，TEは非線形的な変動をする信号 間にも適用できることが分かっている．そこで本研究では，信号間の因果関係の推定法として TEを採用することにした．

3.2

移動エントロピーの概要

2つの変数XとY について，時刻tにおける実現値をそれぞれx(t)，y(t)とすると，変数X から変数Y へ時間τ 後に移動するエントロピーT E(X, Y, τ )は(3.1)式のように定義される． なお，<・>tは時間平均を表す． T E(X, Y, τ ) = ⟨

log₂ pY|Y X(y(t + τ )|y(t), x(t))

pY|Y(y(t + τ )|y(t))

⟩

t

(3.1)

ここで，p_{Y|Y X}(y(t + τ )|y(t), x(t))とp_Y|Y(y(t + τ )|y(t))は，条件付き確率密度関数であり，

Bayesの定理を用いて結合確率密度関数で表すことができる．(3.2)式にp_{Y|Y X}(y(t+τ )|y(t), x(t))

の場合を示す．さらに(3.2)式において，pY Y X(y(t + τ ), y(t), x(t))とpY X(y(t), x(t))は結合確 率密度関数であり，(3.3)式で示すカーネル推定法を用いて計算する．N はデータのサンプル 数，||・||maxは最大距離，Θは(3.4)式で表されるステップ関数である． pY|Y X(y(t + τ )|y(t), x(t)) = pY Y X(y(t + τ ), y(t), x(t)) pY X(y(t), x(t)) (3.2) pY Y X(y(t + τ ), y(t), x(t)) = 1 N − τ N_∑−τ t′₌₁ Θ     r− y(t′+ τ )− y(t + τ) y(t′)− y(t) x(t′)− x(t) max      (3.3) Θ(z) =      1, for z >_{= 0} 0, for z < 0 (3.4) また，(3.3)式のr は精度半径であり，その決定方法は，SROT法 [27]，SNR法 [28]，OS

法[29]，Plug-In法[30]などがあるが，本研究では，SROT法を採用した．SROT法は，(3.5)

式のように定義される．ただし，Aは適用する時系列データの標準偏差と四分位範囲を1.34で 除した値のうち小さい方の値とし，nはデータのサンプル数を表す． r = 0.45A n15 (3.5) さらに，(3.6)式に示すように，0s <_{= τ <= 5s}におけるT E(X, Y, τ )の最大値をT Emax(X, Y ) と定義した．

T Emax(X, Y ) = max 0s<_{=τ <=5s}T E(X, Y, τ ) (3.6) 2変数X，Y 間で，正味の情報移動が有意に存在する，すなわち(3.7)式が成立するならば， Xは原因的で，Y は結果的と推察することができる．これは，情報理論の観点からは，因果関 係は情報移動の非対称性であるとみなされるためである．このようにして，TEを用いて因果 関係の推定を行った．

T Emax(X, Y ) > T Emax(Y, X) > 0 (3.7) 次に，TEを図を用いて説明する．変数Xのエントロピーを(3.8)式のように定義すると， TEは(3.9)式のように表される．なお，Y′は遅延τ後のY を表す． H(X) =−∑ t p(x(t)) log₂p(x(t)) (3.8) T E(X, Y, τ ) = ⟨

log₂pY Y X(y(t + τ ), y(t), x(t))・pY(y(t))

pY X(y(t), x(t))・pY Y(y(t + τ ), y(t))

⟩ t =−H(X, Y, Y′) + H(X, Y ) + H(Y, Y′)− H(Y ) (3.9) これらの関係を図示したものが図3.1である．この中で，T E(X, Y, τ )ははじめはY には無 く，かつ遅延τ後にXからY に移動した情報を意味する．このように，TEはある時間経過後 に，ある変数から別の変数に移動した情報の絶対量を表していることが分かる． なお，TEは対象とする変数を3変数以上に拡張でき，それを内因性移動エントロピー（ in-trinsic transfer entropy : ITE）と呼ぶ．ITEは(3.10)式のように定義され，因果関係を調べ たい2変数以外の変数による情報の回り込みを除くことができる．また，TEと同様，(3.11)

式に示すように，0s <_{= τ <= 5s}におけるIT E(X, Y, τ ; Z)の最大値をIT Emax(X, Y ; Z)と定義 した．

IT E(X, Y, τ ; Z) =

⟨

log₂ pY|Y ZX(y(t + τ )|y(t), z(t), x(t))

pY|Y(y(t + τ )|y(t), z(t)) ⟩ t (3.10) IT Emax(X, Y ; Z) = max 0s<_{=τ <=5s}IT E(X, Y, τ ; Z) (3.11) 一方，1つの変数が持っている情報量はその変数によって異なるため，異なる個体から得ら れる生理的信号のTEの値そのものを単純に比較することはできない．したがって，(3.12)∼

図3.1 変数X, Y, Y′それぞれのエントロピーH(X), H(Y ), H(Y′)とT E(X, Y, τ )の関係図． 赤い部分がT E(X, Y, τ )を表す． (3.15)式のように対象とする全ての変数の結合エントロピーで除することでTE，ITEの規格 化を行った（図3.2）．この規格化により，関連する信号全ての結合エントロピーに対し，その 中で移動した情報量の割合がどの程度であるかを示すことができる．なお，4章以降において TEやITEを表す場合は，この方法により規格化されたものを用いている．

normT E(X, Y, τ ) = T E(X, Y, τ )

H(X, Y, Y′) × 100[%] (3.12)

normIT E(X, Y, τ ; Z) = IT E(X, Y, τ ; Z)

H(X, Y, Y′, Z) × 100[%] (3.13) normT Emax(X, Y ) = T Emax(X, Y ) H(X, Y, Y′) × 100[%] (3.14) normIT Emax(X, Y ; Z) = IT Emax(X, Y ; Z) H(X, Y, Y′, Z) × 100[%] (3.15)

図 3.2 変数X, Y, Y′の結合エントロピーH(X, Y, Y′)とT E(X, Y, τ )の関係図．青の斜線部

3.3

結合写像格子系を用いた移動エントロピーのシミュレーション

3.2節を踏まえ，実際にTEを用いたシミュレーションを行った．具体的には，Schreiber [24] による結合写像格子系を用いたTEのシミュレーションの追試を行った． まず，一方向の結合写像の一次元格子を(3.16)式のように定義する．過程mが1増加する 度に新たな時系列が生成される．nはステップ数または時間を表し，ϵは結合強度を表す．し たがって，xm n は過程mにおけるnステップ目の格子点を表す．図3.3に示したように，この 結合写像格子系は1ステップ前の値xm_n−1とxm_n を用いて次のステップの値xm_n+1を決定してお り，mが増加する方向にのみ情報が移動するようになっている．写像関数fは(3.17)式で定義 されるテント写像を用いた．図3.4は，テント写像関数をグラフにしたものである． xm_n+1= f (ϵxm_n−1+ (1− ϵ)xm_n) (3.16) f (x) =      2− 2x, for 0.5 <_{= x <= 1} 2x, for 0 <_{= x < 0.5} (3.17) 図3.3 結合写像格子系．mは過程を表し，nはステップ数または時間を，ϵは結合強度を表す． 次に，この結合写像格子系において，各過程及び各ステップの初期値に0から1までの乱数 を与えて，時系列を生成した．ある過程における時系列の例を図3.5に示した．その後，結合 写像格子系によって生成された時系列に対し，次の条件を用いて粗視化を行った．

0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x f(x) 図3.4 テント写像関数．      if 0.5 <_{= x <= 1,} then x = 1 if 0 <_{= x < 0.5,} then x = 0 そして，過程mでの時系列をImとし，τ = 1におけるIm−1からImへのTE（以下，T E(Im−1 →Im)と表記する）を計算した．T E(Im−1→Im)を計算する際には，m = 100における105 ステップ経過後の105ステップ間の時系列を用いた．この計算を10回繰り返し，その平均値 を図3.6に示した．図3.6の青線はT E(Im−1→Im)の平均値とその標準誤差，黒線はT E(Im →Im−1)の平均値とその標準誤差，赤線はSchreiberによる理論曲線，水色の実線は追試結果 における最小二乗法によって求めた近似曲線である．理論曲線と近似曲線を，それぞれ(3.18)， (3.19)式に示した． Schreiberによる理論曲線：T E(Im−1→Im) = 0.8554ϵ2 (3.18) 追試結果による近似曲線：T E(Im−1→Im) = 0.9929ϵ2− 0.0073ϵ + 0.4895×10−4 (3.19) 以上より，Schreiberによるシミュレーション結果と同等の結果が得られたことが分かる．

1 1.0001 1.0002 1.0003 1.0004 1.0005 1.0006 1.0007 1.0008 1.0009 1.001 x 105 0 0.2 0.4 0.6 0.8 1 n x 図 3.5 テント写像関数を用いた結合写像格子系によって生成された時系列の例（ϵ = 0.01）． また，ϵが増加すると，1つ前の過程の値x_nm−1によるxm_n+1への寄与が大きくなるため，Im−1 からIm_{への情報の移動量が増加し，T E(I}m−1_→Im_{)も増加する．一方，I}m_からIm−1_にお いて，情報の移動は生じていないので，T E(Im→Im−1)は図3.6のように，ほぼ0に近い値と なっていた．この結果から，TEは2つの信号間における情報の移動量と移動方向を推定でき ることが示された．

0 0.01 0.02 0.03 0.04 0.05 0 0.5 1 1.5 2 2.5 x 10 -3 ε TE[bits] 図 3.6 テント写像関数を用いた結合写像格子系によって生成された時系列から計算された TEのシミュレーションの追試結果．青線はT E(Im−1→Im)の平均値とその標準誤差，黒線は T E(Im→Im−1)の平均値とその標準誤差，赤線はSchreiberによる理論曲線，水色の実線は追 試結果における近似曲線である．

第

4

_{章 生体循環系モデルにおけるシミュレー}

ション

本章では，Seidelら[31]による生体循環系モデルを用いて血圧や心拍数などの生体データを 取得し，ITEを用いることで，血圧調整システム内の信号間における因果関係の推定が行える かを検証した．

4.1

シミュレーションに用いた生体循環系モデルの概要

第2章で述べたように，生体循環系は中枢性調節機構を総動員しており，その上，個々の機 構が協調して血圧をある一定の状態になるように制御している．心臓では図4.1に示した通り， 交感神経と副交感神経により制御が行われており，血管は，主に交感神経のみにより制御が行 われている．Seidelらによる生体循環系モデルでは，これらを模擬したモデルとなっている． 図4.1 圧受容器反射の関係図 [31]．心臓は交感神経と副交感神経の支配を受け，血管は主に 交感神経のみの支配を受けている．

その生体循環系モデルの概要は図4.2に示した通りであり，一つ一つの要素についてそれぞれ

(4.1)∼(4.14)式のように，微分方程式や飽和要素，遅延要素を含んだ式で定義されている [32]．

以下，(4.1)∼(4.14)式について詳細を述べる． 2.1.2項で述べたように，圧受容器は血圧の変化を神経活動に変換している．したがって，血 圧の時間変化（1階微分）を用いて圧受容器の感受性vbを(4.1)式のように定義している．な お，pは血圧であり，k1，k2，p(0)は個体によって異なる定数である． vb= k1(p− p(0)) + k2 dp dt (4.1) 髄質は圧受容器の活動と呼吸による神経活動を受け，交感神経および副交感神経活動を生成 することから，交感神経活動vs，副交感神経活動vpを(4.2)，(4.3)式のように定義している． なお，frは呼吸周波数であり，v(0)s ，ksb，ksr，∆ϕrs，v (0) p ，kbp，krp，∆ϕrpは個体によって異なる 定数である． vs= max(0, v(0)s − kbsvb+ krs| sin(πfrt + ∆ϕrs)|) (4.2) vp = max(0, vp(0)+ kpbvb+ kpr| sin(πfrt + ∆ϕrp)|) (4.3) 交感神経の伝達物質（ノルアドレナリン）のダイナミクスは比較的低速である．この特性を 踏まえ，ノルアドレナリンの濃度変化を(4.4)，(4.5)式のように定義している．ccN a は心臓の ノルアドレナリン濃度を，cvN aは血管のノルアドレナリン濃度を表す．なお，θcN a，θvN a は遅 延時間であり，τcN a，k s c_{cN a}，τvN a，k s c_{vN a} は個体によって異なる定数である． dccN a dt =− ccN a τcN a + k_cs cN avs(t− θcN a) (4.4) dcvN a dt =− cvN a τ_{vN a} + k s c_{vN a}vs(t− θvN a) (4.5) 心臓のノルアドレナリン濃度と迷走神経活動は洞結節の位相速度によって決定されることか ら，洞結節の位相ϕを(4.6)式のように定義している．なお，T(0)は個体によって異なる定数 である． dϕ dt = 1 T(0)fsfp (4.6)

ここで，fs，fpはそれぞれ(4.7)，(4.8)式で与えられる．なお，θpは遅延時間であり，k_ϕcN a， ˆ ccN a，ncN a，k p ϕ，vˆp，npは個体によって異なる定数である． fs= 1 + kϕcN a ( c_{cN a} + (ˆc_{cN a}− c_{cN a}) (ccN a) n_{cN a} (ˆccN a) n_{cN a + (c} cN a) n_{cN a} ) (4.7) fp = 1− kpϕ ( vp(t− θp) + (ˆvp− vp(t− θp)) (vp(t− θp))np ˆ vnp p + (vp(t− θp))np ) F (ϕ) (4.8) また，洞結節の感度は位相効果曲線と呼ばれる関数Fによって(4.9)式のように定義される． この関数は，Reinerら[33]による洞結節の複雑なイオンモデルを用いて得られた関数である． F (ϕ) = ϕ1.3(ϕ− 0.45) (1− ϕ) 3 (1− 0.8)3_{+ (1}− ϕ)3 (4.9) 心拍出量は心臓のノルアドレナリン濃度と1拍前の心周期間隔に影響を受ける．このことを 踏まえて，心収縮性を(4.10)，(4.11)式のように定義している．なお，Ti−1は1拍前の心周期 間隔であり，S(0)，k_Sc，kt_S，Sˆ，nSは個体によって異なる定数である． S_i′ = S(0)+ k_ScccN a+ k t STi−1 (4.10) Si= Si′+ ( ˆS− Si′) (S_i′)nS (S_i′)nS + ( ˆS)nS (4.11) 大動脈は，心収縮による血圧の指数関数的な減衰を表したWindkesselモデルによって表現 される [34]．ノルアドレナリンの血中濃度変化が末梢抵抗に影響を及ぼすことを表すために， Windkesselモデルの時定数τvを(4.12)式のように定義している．なお，τv(0)，τ¯v，cˆvN a，nvN a は個体によって異なる定数である． τv = τv(0)− ¯τv ( cvN a+ (ˆcvN a− cvN a) (cvN a) n_{vN a} (ˆcvN a) n_{vN a + (c} vN a) n_{vN a} ) (4.12) 血圧は，収縮期血圧と拡張期血圧に分けて定義している．本研究では，収縮期の時間をτsys とした．心収縮性に応じた脈圧によって血圧が増加するように，収縮期血圧p1は滑らかな曲 線を用いて(4.13)式のように定義している．なお，di−1は収縮開始時の血圧値であり，tiはそ の時刻である．

p1 = di−1+ Si t− ti τsys exp { 1−t− ti τsys } (4.13) Windkesselモデルに基づき，拡張期血圧p2は指数関数的に減少するよう(4.14)式のように 定義している． dp2 dt =− p2 τv(t) (4.14) 微分方程式を解く際には，4次のRunge-Kutta法を用いて計算を行った [35]．Runge-Kutta 法とは，常微分方程式の近似解を求める手法の1つであり，ある点において，既知の値と傾き から次の値を計算する数値的解法である．次に4次のRunge-Kutta法の計算方法について述 べる． 1階の常微分方程式 dy dt = f (t, y) (4.15) において，y(t0) = y0，tの増分をhとすると， yn+1 = yn+ 1 6(k1+ 2k2+ 2k3+ k4) (4.16) k1 = hf (tn, yn) (4.17) k2 = hf (tn+ h 2, yn+ k1 2) (4.18) k3 = hf (tn+ h 2, yn+ k2 2) (4.19) k4 = hf (tn+ h, yn+ k3) (4.20) のようにして，次の値yn+1を求めることができる．

4.2

シミュレーションの概要

4.1節で説明した生体循環系モデルの定数を変化させ，若年者と高齢者の生体循環系を模擬 したモデルを作成した。そして血圧，心拍数，血管運動と関連する信号の3つの信号間のITE を求めることで，血圧調整システム内における各信号間の因果関係を調べ，若年者モデルと高 齢者モデルとで違いが見られるかを検証した． 若年者と高齢者の生体循環系を模擬したモデルを作成するにあたり，圧受容器反射機能は年 齢に伴い低下する[36, 37]ことを考慮した．

まず，心臓系においては，交感神経活動と副交感神経活動を表すvs，vpにおいて圧受容器反 射機能のゲインを表すkb_s，kb_pを若年者モデルの方が大きくなるように設定した．そして，血 管系においては，Windkesselモデルの時定数τvを高齢者モデルの方が大きくなるように設定 した．若年者モデルと高齢者モデルにおける各定数の設定値を表4.1，表4.2に示した． 表4.1 若年者を模擬した生体循環系モデルの各定数の設定値．赤字は若年者を模擬するため に設定値を変更した定数である． p(0) 50.0mmHg θp 0.5s k1 0.02mmHg−1 S(0) 25.0mmHg k2 0.00125smmHg−1 kSc 40.0mmHg vs(0) 0.8 k_St 10.0mmHgs−1 kb_s 0.75 nS 2.5 kr_s 0.1 Sˆ 70.0mmHg fr 0.2s−1 T(0) 1.1s ∆ϕr s 0.0 kϕcN a 1.6 vp(0) 0.0 ˆccN a 2.0 kb p 0.35 ncN a 2.0 kr_p 0.1 k_ϕp 5.8 ∆ϕr p 0.0 ˆvp 2.5 τcN a 2.0s np 2.0 ks_c cN a 1.2 τ (0) v 1.8s θ_{cN a} 1.65s ¯τv 1.25s τ_{vN a} 2.0s ˆc_{vN a} 10.0 ks_c vN a 1.2 nvN a 1.5 θ_{vN a} 1.65s τsys 0.125s ITEの対象とした信号は，血圧は拍内平均血圧M BP を，心臓系では心拍数HRを，血管 系では収縮期血圧と拡張期血圧の差である脈圧P P とした．血管系の信号としてP P を用いた 理由は，以下で述べるように総末梢血管抵抗T P Rと関連しているためである[38]． (4.21)式のように，総末梢血管抵抗T P Rは拍内平均血圧M BP を心拍出量COで除するこ とで算出できる． T P R = M BP CO (4.21) 血管コンプライアンスCは血圧の変化∆P とその際の血管の容積変化∆V を用いて，C = ∆V /∆P と表すことができ，また∆P は脈圧P P，∆V は一回拍出量SV で近似できること

表4.2 高齢者を模擬した生体循環系モデルの各定数の設定値．赤字は高齢者を模擬するため に設定値を変更した定数である． p(0) 50.0mmHg θp 0.5s k1 0.02mmHg−1 S(0) 25.0mmHg k2 0.00125smmHg−1 kSc 40.0mmHg vs(0) 0.8 k_St 10.0mmHgs−1 kb_s 0.6 nS 2.5 kr s 0.1 Sˆ 70.0mmHg fr 0.2s−1 T(0) 1.1s ∆ϕr_s 0.0 k_ϕcN a 1.6 vp(0) 0.0 ˆccN a 2.0 kb_p 0.2 ncN a 2.0 kr_p 0.1 k_ϕp 5.8 ∆ϕr_p 0.0 ˆvp 2.5 τ_{cN a} 2.0s np 2.0 ks_c cN a 1.2 τ (0) v 2.4s θ_{cN a} 1.65s ¯τv 1.0s τ_{vN a} 2.0s ˆc_{vN a} 10.0 ks_c vN a 1.2 nvN a 1.5 θ_{vN a} 1.65s τsys 0.125s から， C = ∆V ∆P≒ SV P P (4.22) と表される． (4.23)式のように，心拍出量COは一回拍出量SV と心拍数HRの積で表される． CO = SV × HR (4.23) 以上，(4.21)∼(4.23)式から， M BP = C× P P × HR × T P R (4.24) ∴ P P = M BP T P R× C × HR (4.25) となる． シミュレーションにより得られたデータのうち，定常状態時におけるM BP，HR，P P の 500秒間のデータを解析に用いた．シミュレーションにより得られた血圧波形，M BP，HR，

P P の一例を図4.3に示した．なお，ITE解析前の処理として，拍ごとのデータである各信号 に，再サンプリング周波数を2Hzとして3次のスプライン補間を施し，補間後の各信号を0.04 ∼0.15Hzを通過域とする帯域通過ディジタルフィルタ（3次バターワース型）に通した．そし てフィルタ処理後の各信号を平均0，分散1に正規化し，M BP，HR，P P の3つの信号間にお けるIT Emaxを求めた．M BP，HR，P P の3つの信号間の関係性を見る場合，M BP →HR， HR→M BP，HR→P P，P P→HR，P P→M BP，M BP →P P の計6通りの方向がある． その6通りのIT Emaxの中で最も大きい値をとる方向を最大移動エントロピー方向（maximam

transfer entropy direction: MTED）とすると，MTEDが3つの信号間において最も確率的な 因果関係が強い方向であるとみなすことができる．

図4.3 若年者モデルでのシミュレーションにより得られた血圧波形，平均血圧（M BP），心 拍数（HR），脈圧（P P）の例．

図4.4 高齢者モデルでのシミュレーションにより得られた血圧波形，平均血圧（M BP），心 拍数（HR），脈圧（P P）の例．

4.3

シミュレーションの結果

若年者モデルと高齢者モデルそれぞれ10回シミュレーションを行い，MTEDの割合を図4.5 に，図4.5を基に3つの信号間における因果関係を表したものを図4.6に示した．図4.5から 分かるように，若年者モデルでは，100%の割合でMTEDはM BP →HRであった．これは， 圧受容器反射機能のゲインを表す定数を大きくしたことにより，すなわち圧受容器の感受性を 高めたことにより，HRがよりM BP を強く反映するようになったためであると考えられる． 高齢者モデルにおいては，圧受容器反射機能のゲインを小さくし，圧受容器の感受性を低下さ せたため，M BP →HRのIT Emaxが小さくなり，相対的にM BP がHRの影響をより受け やすくなったため，MTEDがHR→M BP となったと考えられる．このように，年齢による 圧受容器反射機能や血管特性の変化によって，血圧調整システム内における信号間の因果関係 が変わってくることを，ITEを用いることで評価できる可能性が示唆された． 図 4.5 若年者モデルと高齢者モデルそれぞれにおけるMTEDをとる割合．

図4.6 （a）若年者モデルと（b）高齢者モデルのM BP，HR，P P間における因果関係．太 線は6通りの情報の流れの中で最も強い因果関係があった方向を表す．

第

5

_{章 実際の生体データを用いた血圧調整シス}

テム内の信号間における因果関係の推定

第4章では，生体循環系モデルを用いて血圧調整システム内の信号間における因果関係を推 定するシミュレーションを行ったが，実際の生体データの場合について検証するために，若年 者と高齢者に対して行われた呼吸統制実験におけるデータを用いて解析を行った．

5.1

血管系の信号

心臓の心室の収縮によって大動脈基部に血液が急激に送り込まれると，この部分の内圧は 上昇し，血管壁は引き伸ばされて膨れる．こうして大動脈基部とその隣接部の間に圧勾配が 生じる．大動脈基部を膨らませていた血液は，この圧勾配によって加速され，隣接部に流れ 込んでこの部分の血管壁を膨らます．このような変化が弾性管である血管に沿って末梢へと 伝わっていく．すなわち，内圧上昇と壁の伸展が波動として末梢に伝播していく．これを脈波 と呼ぶ．脈波は，光を利用して体表面付近の血流量変化を非侵襲的に測定する光電容積脈波 （photoplethysmography: PPG）のように，容易に測定を行うことができる．脈波の波形や伝播 の速度は，心臓の拍出特性のみならず血管の形状や血管壁，血液の物理的特性に依存する[39]． 第4章でのシミュレーションでは血管系の信号として脈圧PPを用いていたが，PP以外に もいくつか考えられる．これらについて，以下で述べる．

5.1.1

脈波の振幅

脈波の振幅（pulsatile amplitude: PA）は，脈波の拍動成分であり，心臓の収縮期における 組織血液量の増大に依存し，動脈の血管コンプライアンス，動脈抵抗によって決定される．図 5.1に典型的な脈波波形とPAの定義を示した．

図 5.1 脈波波形とPAの定義．

5.1.2

脈波伝播時間

脈波伝播時間（pulse transmission time: PTT）は大動脈基始部で脈波が発生してから指先 に到達するまでの時間として定義され，血圧と負の相関にあるとされる．これは，血圧の変動 によって血管のコンプライアンスが変化し，脈波の伝達特性が変化するためである[40]．PTT は，図5.2に示したように心電図（electrocardiogram: ECG）のR波ピークから末梢での脈波 であるPPGの立ち上がり点までの時間を求めることで算出することができる．PPGの立ち上 がり点は明確な定義がなされておらず，その決定法は複数存在する[41]．以下にその代表的な ものを述べる． • ボトム法：PPGの値が最小となる点 • トップ法：PPGの値が最大となる点 • 1次微分法：速度脈波（PPGの1階微分）の最大点 • 2次微分法：加速度脈波（PPGの2階微分）の最大点 このうち，本研究では1次微分法を用いてPTTの算出を行った．

図5.2 脈波伝播時間PTTの定義．

5.1.3

血圧と脈波伝播時間の関係式

脈波が血管内を伝播する速度v[m/s]は血管の寸法と血管壁の柔らかさに依存する．解析的 な関係式は(5.1)式で示されるMoens-Kortewegの式によって表される[42]． v = √ ghE ρd (5.1) ここでg[m/s2]は重力定数，Eは血管のヤング率，ρ[g/m3]は血液の比重，h[m]は血管壁厚， d[m]は血管の直径である． 血圧が上昇すると，hは減少し，dは増加するので，vは減少するように思えるが，Hughes らによればEは次式で表されるように血圧P [mmHg]に指数関数的に増加し，この影響が形状 変化の影響を上回るので，実際にはvは増加する． E = E0eγP (5.2) ここで，E0はP = 0におけるヤング率，γ[mmHg−1]は血管に依存する定数であり，その値 の範囲は0.016∼0.018[mmHg−1]である． 脈波伝播時間P T T [s]は血管経路長さをL[m]とすれば(5.3)式のように表すことができる． P T T = L v (5.3)

(5.1)∼(5.3)式より， ghE0e γP ρd = L2 P T T2 (5.4) ∴ d h · ρL2 gE0 = eγPP T T2 (5.5) となる．このように，血圧と脈波伝播時間における関係は，(5.5)式のように表すことがで きる．(5.5)式の左辺は，血管の直径dと血管壁厚hを含んでおり，血管運動と関連する信号 と考えられる．したがって，血管系の信号M Kを(5.6)式のように定義する． M K = eγPP T T2 (5.6)

5.1.4

血管コンプライアンスと総末梢血管抵抗

血管コンプライアンスCと総末梢血管抵抗T P Rの積CRは(4.25)式から以下のように表さ れ，M BP，HR，P P から近似的に算出することができる． CR = C× T P R≒ M BP P P × HR (5.7) このCRは，Windkesselモデルにおける時定数を表しており，血管運動をよく反映している と考えられる．

5.2

実験の概要

被験者は椅子に座った状態で十分な安静を取った後，安静座位の状態で7分間，生体信号を 連続的に測定した．生体信号は，心電図，連続血圧および光電容積脈波の3種類であり，表5.1 に示した測定装置等を用いて測定を行った．被験者は7分間の測定中に，自由呼吸と2種類の 呼吸統制を行うが，そのプロトコルの詳細を図5.3に示した．まず測定開始から3分間は自由 呼吸を行い，その後の2分間は10秒周期（呼気6秒，吸気4秒）の呼吸統制を行う．さらに最 後の2分間は4秒周期（呼気2秒，吸気2秒）の呼吸統制を行い，7分間の測定を終える．呼 吸統制は視覚と音でタイミングを合わせることが可能なPC上のソフトウェアを用いた．測定 のサンプリング周波数は全て1kHz とした．

被験者は若年者59名（男性48名，女性11名；25.7±6.3歳），高齢者96名（男性93名， 女性3名；69.8±4.0歳）の計155名であった．1人の被験者につき日および時間帯を変えて 最大で4回の測定を行い，総データ数は381例であった． 表 5.1 測定条件 測定信号 測定装置 アンプ 心電図 エーカークリップ_{（フクダ電子） （}ECG100C_BIOPAC） 連続血圧 Portapres （TNO-TPD Biomedical Instrumentation）

-光電容積脈波 Nellcor DS-100A_{（Envitec） （}PPG100C_BIOPAC）

図 5.3 呼吸統制実験のプロトコル．

5.3

解析方法

M BP，HRの他に，以下に示した8種類の血管系の信号の中から1つを用いて，3つの信号 間におけるIT Emaxを求めることにより，解析を行った．なお，解析前の処理として，拍ごと のデータである各信号に，再サンプリング周波数を2Hzとして3次のスプライン補間を施し， 補間後の各信号を0.04∼0.15Hzを通過域とする帯域通過ディジタルフィルタ（3次バターワー ス型）に通した．そしてフィルタ処理後の各信号を平均0，分散1に正規化し，IT Emaxを計 算した． 1）P P：収縮期血圧と拡張期血圧の差である脈圧． 2）P A：5.1.1項で述べた脈波の振幅．

3）P T T：5.1.2項で述べた脈波伝播時間． 4）M K1：(5.6)式において(5.8)式のようにP = M BP として算出した信号． M K1 = eγM BPP T T2   (5.8) 5）M K2：M K1の自然対数をとった信号．(5.9)式のように定義される． M K2 = ln M K1 = γM BP + 2 ln P T T   (5.9) 6）M K3：M K1の逆数をとった信号．(5.10)式のように定義される． M K3 = 1 M K1 = e −γMBP_{P T T}−2_(5.10) 7）CR：(5.7)式で定義された信号． 8）1/CR：CRの逆数をとった信号．(5.11)式のように定義される． 1/CR = P P × HR M BP    (5.11) ある信号の自然対数や逆数も対象とした理由は，そのような処理を行ってもその信号が有す る情報量は変化しないことを検証するためである． これらの各血管系の信号についてM BP とHRと共に，ある被験者の自由呼吸時，10秒呼 吸時，4秒呼吸時における信号を図5.4∼図5.6に示した．なお，各信号は平均0，分散1に正 規化して表示した．

40 50 60 70 80 90 100 -5 0 5 PP 40 50 60 70 80 90 100 -5 0 5 PA 40 50 60 70 80 90 100 -5 0 5 PTT 40 50 60 70 80 90 100 -5 0 5 MK1 40 50 60 70 80 90 100 -5 0 5 MK2 40 50 60 70 80 90 100 -5 0 5 MK3 40 50 60 70 80 90 100 -5 0 5 CR 40 50 60 70 80 90 100 -5 0 5 1/CR 時間[s] MBP HR 血管系の信号 図5.4 自由呼吸時における各信号の例．青線がM BP，緑線がHR，赤線が血管系の信号で ある．

220 230 240 250 260 270 280 -5 0 5 PP 220 230 240 250 260 270 280 -5 0 5 PA 220 230 240 250 260 270 280 -5 0 5 PTT 220 230 240 250 260 270 280 -5 0 5 MK1 220 230 240 250 260 270 280 -5 0 5 MK2 220 230 240 250 260 270 280 -5 0 5 MK3 220 230 240 250 260 270 280 -5 0 5 CR 220 230 240 250 260 270 280 -5 0 5 1/CR 時間[s] MBP HR 血管系の信号 図5.5 10秒呼吸時における各信号の例．青線がM BP，緑線がHR，赤線が血管系の信号で

320 330 340 350 360 370 380 -20 2 PP 320 330 340 350 360 370 380 -20 2 PA 320 330 340 350 360 370 380 -20 2 PTT 320 330 340 350 360 370 380 -20 2 MK1 320 330 340 350 360 370 380 -2 0 2 MK2 320 330 340 350 360 370 380 -20 2 MK3 320 330 340 350 360 370 380 -20 2 CR 320 330 340 350 360 370 380 -2 0 2 1/CR 時間[s] MBP HR 血管系の信号 図5.6 4秒呼吸時における各信号の例．青線がM BP，緑線がHR，赤線が血管系の信号で

5.4

解析結果と考察

5.4.1

血圧調整システム内の 3 つの信号間における因果関係

まず，若年者と高齢者の血圧調整システム内の3つの信号間における6方向のIT Emaxの中 で最も大きい値の方向MTEDをとる割合を，8種類の血管系の信号別に図5.7∼図5.14に示し た．図5.10∼図5.14を見ると，M K1，M K2，M K3において，またCR，1/CRにおいてほ ぼ似たような結果が得られ，予想した通り，ある信号の自然対数や逆数をとってもその信号が 有する情報量はほとんど変化しないことが分かった． 次に，ある1つの方向における若年者と高齢者のMTEDの割合の差の絶対値を求め，計6 方向の総和をとったものを8種類の血管系の信号ごとに図5.15に示した．これにより，8種類 の血管系の信号の中でどれを用いた場合が，最も若年者と高齢者との間における違いを見られ るかが調べられる．図5.15を見ると差の絶対値の総和はP P を用いた場合が最大であったの で，M BP，HR，P P の3つの信号間のMTEDをとる割合を調べることで，若年者と高齢者 における血圧調整システムの違いを調べられる可能性が高いと考えられる．したがって，これ 以降はP P を用いた場合のITEについて考察を進めていく． 図5.7 M BP，HR，P P の3つの信号間の若年者と高齢者それぞれにおけるMTEDをとる 割合．

図5.8 M BP，HR，P Aの3つの信号間の若年者と高齢者それぞれにおけるMTEDをとる 割合．

図5.9 M BP，HR，P T T の3つの信号間の若年者と高齢者それぞれにおけるMTEDをと

図 5.10 M BP，HR，M K1の3つの信号間の若年者と高齢者それぞれにおけるMTEDを とる割合．

図 5.11 M BP，HR，M K2の3つの信号間の若年者と高齢者それぞれにおけるMTEDを とる割合．

図 5.12 M BP，HR，M K3の3つの信号間の若年者と高齢者それぞれにおけるMTEDを とる割合．

図 5.13 M BP，HR，CRの3つの信号間の若年者と高齢者それぞれにおけるMTEDをと る割合．

図5.14 M BP，HR，1/CRの3つの信号間の若年者と高齢者それぞれにおけるMTEDを とる割合．

図5.15 8種類の血管系の信号ごとの若年者と高齢者におけるMTEDの割合の差の絶対値の 総和．

まず，シミュレーションの結果と比較する．図4.5，図5.7を見ると，シミュレーションと実験 データのどちらの場合においてもMTEDをとる割合は，若年者はM BP →HRの方向が最も 大きく，高齢者はHR→M BP の方向が最大であった．このパターンに当てはまった被験者は， 若年者については，圧受容器反射機能が高い被験者であり，高齢者については，圧受容器反射機 能が低い被験者である可能性が高いと考えられる．また，シミュレーションではM BP →P P の方向が最も大きくなることは無かったが，実験データでは若年者においてM BP →HRと 同程度の割合であった．M BP →HRの因果関係というのは，圧受容器反射により心拍数を調 整していることを意味し，M BP →P P の因果関係というのは，圧受容器反射により血管を調 整していることを意味すると考えられる．したがって，この2つの方向は主に圧受容器反射に よる働きを表していると考えられ，実際に，圧受容器反射機能が低下している高齢者において はM BP →HRとM BP →P P のいずれにおいても強い因果関係がほとんど見られなかった． しかし，高齢者においても圧受容器反射は機能しているので，M BP →HRとM BP →P P の方向の情報移動は存在していると考えられる．高齢者においてHR→M BP やHR→P P の方向のMTEDとなる割合が大きかった理由は，循環系の物理的要因による因果関係を表す HR→M BP やHR→P P のIT Emaxが相対的に大きくなったためであると考えられる． 次に，若年者，高齢者別に従来の圧受容器反射感度の指標であるαLF で被験者を分類（0 <₌ αLF < 5，5 <_{= α}LF < 10，10 <_{= α}LF）した場合におけるMTEDをとる割合をそれぞれ図 5.16，図5.18に，3つの信号間における因果関係を表したものをそれぞれ図5.17，図5.19に示 した．なお，若年者と高齢者におけるαLF の分布は図5.20のようになっており，また図5.21 に示した通り，両者のαLFの平均値には有意差がみられた（p < 0.001）．図5.16，図5.18を見 ると，αLFが同程度であっても若年者と高齢者とでMTEDをとる傾向が異なっていた．前述 したようにαLF という指標は圧受容器反射機能のゲインを表す指標であり，M BP →HRや M BP →P P の方向の因果関係と同じ傾向になると考えられる．若年者に関してはこの予想通 り，αLF が大きい被験者ほどMTEDがM BP →HRやM BP→P P である割合が大きくなっ ていた．しかし高齢者においては，αLFが大きい被験者であってもMTEDがM BP →HRや M BP →P P になることはほとんどなかった．この理由は，前述したようにM BP →HRと M BP →P P の方向の情報移動は確かに存在しているが，それらの方向よりもHR→M BP やHR→P P のIT Emaxが相対的に大きくなったためであると考えられる．また，高齢者で αLF が大きかった被験者は圧受容器反射機能のゲインが大きいのではなく，不整脈等により心 拍数変動が大きくなった結果，(2.1)式におけるSRRIが大きくなってしまったことが影響して いる可能性も考えられる．

さらに図5.16において注目すべき点として，若年者ではαLF が大きい群がM BP →HRと M BP →P P に分かれているという点がある．前述したように，M BP →HRの因果関係とい うのは，圧受容器反射による心拍数調整を表し，M BP →P P の因果関係というのは，圧受容 器反射による血管調整を表すと考えられる．したがって，αLF では分類することができなかっ た圧受容器反射による心拍数調整と血管調整の優位性の区別を，ITEを用いることにより評価 できる可能性が示唆された．

図 5.16 若年者において被験者をαLF に基づいて3群（0 ≤ αLF < 5，5 ≤ αLF < 10，

10≤ αLF）に分類した場合におけるMTEDをとる割合．

図 5.17 若年者において被験者をαLF に基づいて3群（0 ≤ αLF < 5，5 ≤ αLF < 10，

10≤ αLF）に分類した場合のM BP，HR，P P 間における因果関係．線が太くなるほど強い 因果関係があったことを表す．

図 5.18 高齢者において被験者をαLF に基づいて3群（0 ≤ αLF < 5，5 ≤ αLF < 10，

10≤ αLF）に分類した場合におけるMTEDをとる割合．

図 5.19 高齢者において被験者をαLF に基づいて3群（0 ≤ αLF < 5，5 ≤ αLF < 10，

10≤ αLF）に分類した場合のM BP，HR，P P 間における因果関係．線が太くなるほど強い 因果関係があったことを表す．

図 5.20 若年者と高齢者のαLF の分布．緑の直線は分類（0 ≤ αLF < 5，5 ≤ αLF < 10，

10≤ αLF）の境界線を表す．