上原記念生命科学財団研究報告集, 35 (2021)

上原記念生命科学財団研究報告集, 35（2021）

82. 遅延光学系と多重時間解析を活用した脳機能定量評価法
小谷 潔
東京大学 先端科学技術研究センター
Key words：時間遅れ力学系，磁気共鳴画像，同期現象，脳血流動態反応

緒 言
時間遅れ力学系は生体、工学、社会システムなどに幅広くみられ、多様なダイナミクスを示すことが知られている。
一方で、時間遅れ力学系は無限次元特性を持つため、解析的に振る舞いを読み解くことは困難である。我々は非線形力学系の縮約理論を時間遅れ力学系に発展させ、モード分解法を提案している。本研究では、基礎的な解析理論として、時間遅れを含む数理モデルの性質を、モード分解を用いた数理解析によって調べ、解析理論をさらに進展させる。特に、近赤外光を照射してその応答から脳機能を読み解くには入出力間の拡散、遅れを解析する必要がある。そのため、入射するレーザ光源については、遅延光学系を活用することで多様なダイナミクスの振動をデザインすることが望まれる。
さらに、近赤外光を用いたヘモグロビン濃度定量評価に向けて、被験者の解剖学的な脳構造を模したシミュレーション
および被験者に近赤外光を照射する実験から頭部の解剖学的な構造と吸収散乱特性評価の関係を解析する。

方法および結果
1．時間遅れ力学系のモード分解による数理解析
1.1 Lang-Kobayashi 方程式の解析
レーザや神経系を含め、時間遅れと非線形効果が相まって周期軌道（リミットサイクル軌道）を生み出す例は多く存在する。しかしながら、そのような系の外力に対する応答を解析するための低次元力学系を導出（縮約）し解析する研究は進んでいない。本研究でははじめに、Lang-Kobayashi 方程式を用いて時間遅れフィードバックによる多重時間強度を持つレーザの数理解析を行った。ここでは簡単のため、文献［1］によって提案されている利得飽和効果などを考慮しない Lang-Kobayashi 方程式を用いる。モデルの式は𝑑𝑑Ε
= (1 + 𝑖𝑖𝑖𝑖)𝑁𝑁Ε + 𝜅𝜅Ε(𝑠𝑠 − 𝜏𝜏) exp(−𝑖𝑖𝜔𝜔0 τ)
𝑑𝑑𝑑𝑑
𝑑𝑑𝑑𝑑 𝑃𝑃 𝑁𝑁 (1 + 2𝑁𝑁)|Ε|2= − −𝑇𝑇
𝑑𝑑𝑑𝑑 𝑇𝑇 𝑇𝑇で与えられ、各パラメータは全て文献［1］に従う。𝜅𝜅は時間遅れフィードバックの強度に関するパラメータである。
はじめに時間遅れフィードバックによる不安定化の様子を図 1 に示す。フィードバックのない𝜅𝜅＝0 ではレーザ強度は一定であり、𝜅𝜅が増えると不安定化され振動状態が発生することが確認された。𝜅𝜅＝1.2×10−3 で発生した振動の周期は約 1.1 ns であり、GHz 帯域での振動が与えられる。さらにこの振動状態を制御するため、𝜅𝜅＝1.2×10−3 のケースにおいて N に対する位相感受関数［2、3］を算出すると図 2 の通りであった。

ここではさらに、弱い注入電流を周期的に振動させることで振動周期のデザインができることを確認する。位相感受関数を用いた理論計算からは𝑃𝑃＝𝑃𝑃0 sin(2Ω𝑡𝑡)において、𝑃𝑃0 ＝10−5 の弱い周期外力に対して 2：1 で同期する領域を求めた。𝑃𝑃0 が十分に小さい場合は位相結合関数を用いた解析から5.67485 × 10−9 ＜Ω＜5.78875 × 10−9 の範囲で同期することが予想される。得られた解析結果を確認するために同期する領域をシミュレーションによって求めたところ、図 3でΩmin ＝5.68595 × 10−9 、Ωmax ＝5.7318 × 10−9 が得られ、理論解析と概ね一致することを確認した。

1

(a)

(b)

(c)

(d)

図 1．時間遅れフィードバックによるレーザダイナミクスの変化（a）𝜅𝜅＝0、（b）𝜅𝜅＝1.2×10－３、（c）𝜅𝜅＝1.5×10－３、（d） 𝜅𝜅＝1.2×10－３の場合について、力学変数 N における位相感受関数。

(a)

(b)

図 2．周期外力に対する同期特性の解析（a）𝜅𝜅 =1.2×10－３における周期外力に対する位相結合関数。
（b）周期外力に対する同期領域。

2

1.2 位相振幅方程式による振幅効果の解析
前節では、外力が十分に弱くシステムがリミットサイクル軌道上で運動していることを前提として解析を行った。一方で、フィードバックレーザの同期制御においてはリミットサイクル軌道をとらない場合や、制御により強い外力を必要とする場合が起こりうる。そのようなより複雑な振動状態での同期制御を行うには、振幅効果を考慮した解析を検討する必要があるが、時間遅れ系に関する振幅効果の解析理論は十分には整備されていない。ここでは理論的な新たに枠組みを拡張するため、モード分解可能な時間遅れ力学系モデルを用いた解析を進め、モード間の非線形な相互作用が同期現象に与える影響を解析した。ここで解析する数理モデルは𝑑𝑑𝑑𝑑(𝑡𝑡)
𝜋𝜋𝜋𝜋 22𝜋𝜋= −𝑥𝑥 (𝑡𝑡 − ) + 𝜖𝜖𝜖𝜖(𝑡𝑡) [1 − 𝑥𝑥(𝑡𝑡)2 − 𝑥𝑥 (𝑡𝑡 − ) ] + 𝐺𝐺0 sin( 𝑟𝑟𝑟𝑟)
𝑑𝑑𝑑𝑑
22𝑇𝑇で表される。本解析ではεをε＜＜1 を満たすパラメータとして与えるが、その場合にはリミットサイクル軌道への緩和が非常に遅い周期振動を有する。そこで、このモデルを用いてモード間の非線形相互作用の影響を評価する。外力の振幅𝐺𝐺0 ＝0 の場合についてはこの微分方程式の解軌道は Functional representation によって𝑥𝑥0 (𝜎𝜎)＝(𝑡𝑡)sin(𝑡𝑡＋𝜎𝜎) (−𝜏𝜏 ≤ 𝜎𝜎 ≤ 0)のリミットサイクル軌道を描き、また数理解析によってその位相感受関数、振幅感受関数は解析的に導出できる。これらの感受関数を用いた位相振幅縮約を行うことで、最終的に位相振幅間の非線形相互作用を取り入れた位相振幅方程式

𝑑𝑑𝑑𝑑
1[𝑎𝑎 𝑅𝑅2 + 𝑎𝑎3 𝑅𝑅3 + 𝑔𝑔0 (𝜓𝜓)]= Δ𝜔𝜔 +𝑑𝑑𝑑𝑑
1 + 𝑅𝑅𝑎𝑎0 2

𝑑𝑑𝑑𝑑
𝑅𝑅(𝑏𝑏0 − 𝜆𝜆1 )[𝑎𝑎2 𝑅𝑅2 + 𝑎𝑎3 𝑅𝑅3 + 𝑔𝑔0 (𝜓𝜓)] + 𝑔𝑔1 (𝜓𝜓)= 𝜆𝜆1 𝑅𝑅 + 𝑏𝑏2 𝑅𝑅2 + 𝑏𝑏3 𝑅𝑅3 −1 + 𝑅𝑅𝑎𝑎0𝑑𝑑𝑑𝑑

この位相振幅方程式の解析から、が得られ、上式の各パラメータが元のパラメータεと外力項𝐺𝐺0 、𝑟𝑟によって表される。

振幅𝐺𝐺0 の周期振動によって同期される領域の広さを解析し、従来の位相方程式による解析と比較を行った。図の赤線が安定平衡点であり、 0.964＜𝑟𝑟＜1.033の範囲で同期することがわかった。一方で位相縮約理論からは同期領域は0.951＜𝑟𝑟＜1.048と解析されるため、時間遅れに誘起された振動状態においても、振幅効果がモード間の非線形相互作用を介して同期状態を不安定化しうることが示された。

前節の位相方程式および本節位相振幅方程式による解析から、時間遅れフィードバックを介したレーザの振動特性を理解するための基礎的な知見を得ることができた。

図 3．非線形位相振幅方程式による同期領域の解析結果
ε＝0.04、𝐺𝐺0 ＝0.05の場合。

3

2．MRI 画像を用いた頭部モデルの作成とシミュレーション
ヒトの頭部に近赤外光を照射しその応答から内部状態を推定するには、個人ごとの解剖学的な構造および複数の組織による異なる光学特性の組み合わせを考慮しなければならない。ここでは磁気共鳴画像（MRI：Magnetic ResonanceImaging）による被験者の解剖画像から頭部モデルを構築し、有限要素法による光散乱特性の数値計算から評価する。
MRI 装置は Siemens の MAGNETOM Prisma を使用した。ここでは通常の T1 強調画像に加え、True-FISP（firstimaging with steady-state precession）シーケンスおよび 3D-FLASH（three-dimension fast low angle shot）シーケンスを用いて解剖学的な頭部構造の検出に活用した。これらの MRI 画像から被験者の頭皮、頭蓋骨、脳脊髄液、灰白質・白質の 4 層を抽出し、生体組織における近赤外光の光伝播解析を行った。生体組織は強い散乱性を示すため、拡散近似を適用し、
1 𝜕𝜕
Φ(𝑟𝑟、𝑡𝑡) − 𝐷𝐷∇2 Φ(𝑟𝑟、𝑡𝑡) + 𝜇𝜇𝑎𝑎 Φ(𝑟𝑟、𝑡𝑡) = 𝑆𝑆(𝑟𝑟、𝑡𝑡)
𝑐𝑐 𝜕𝜕𝜕𝜕

で示される光拡散方程式を有限要素法による数値計算によって評価する。前頭部において、大脳縦裂を跨ぐプローブ配置と跨がない配置について、それぞれ送光プローブに対して受光プローブの位置を 20 mm、40 mm に設置した場合の多点受光プローブによるシミュレーションを行った。頭部プローブの具体的な位置は図 4A、B のとおりとした。

図 4．頭部プローブの配置左が配置 A、右が配置 B。

はじめに波長 830 nm の近赤外光を 1 GHz の振動状態で入力した際に、多点受光プローブから得られる振動の位相差は表 1 の通りであった。表 1 の 1 ch、2 ch ともに入射光との位相差はプローブ配置 A の方が大きく、またその差は2ch の方が大きかった。これは解剖学的な構造の差を反映していると考えられる。一方でその差は 1 ch では微小であるため、実験によって適切な位相差の評価を行うには、送光プローブと受光プローブの距離を離す必要があることが示された。次に、半値幅が 150 ps のパルス光を入力し、多点受光プローブの出力から光拡散方程式による解の波形をCurve Fitting することで平均吸光係数を算出した（図 5）。
表 1．プローブ配置 A および B における位相差（rad）

3．実験による頭部構造に起因する光伝搬特性の違いの評価
前節のシミュレーション結果を検証するために、頭部シミュレーションと同じ位置（図 4A、B の位置）にプローブを設置し、波長 836 nm のパルスレーザを照射した実験を行った。実験においては被験者頭部にプローブを設置した上で安静座位における 300 秒の計測を行い、計測結果を 1 秒間積算したものを 1 ブロックとし、1 ブロックごとに文献4

［4］の手法で Curve Fitting を行い平均吸光係数の実測値を算出した。得られたシミュレーションと実験での平均吸光係数を比較した結果を図 5 に示す。平均吸光係数はプローブ間隔 20 mm（1 ch）の場合よりも 40 mm（2 ch）の方が高く、また配置 A に比べて配置 B の方が高いという結果が得られた。

考 察
本研究においては遅延フィードバックを用いた近赤外レーザの制御に関する理論研究と、被験者頭部モデルの構築とシミュレーション、実験による脳内ヘモグロビン濃度推定に関する評価を行った。以下に考察を記す。
Lang-Kobayashi 方程式の数理解析によって、時間遅れフィードバックを活用することで GHz 帯域での振幅変調を持つレーザ光を設計できること、さらにその周波数は外力によって制御可能であることを示した。また、被験者の頭部形状を模した近赤外レーザシミュレーションからは、GHz 帯域の振幅変調を持つ近赤外光は頭部構造による吸収散乱特性の変化を評価可能であることが示された。さらにパルスレーザを用いた多点受光プローブでの頭部シミュレーションおよび実験から、プローブ間隔 20 mm と 40 mm の比較においては 40 mm の方が平均吸光係数が高くなった。深層の灰白質を経由した帰還光の割合が多いことが平均吸光係数が高くなった原因と考えられ、タスク等に伴う脳血流動態を評価するには 40 mm が適切であることを示唆している。さらに、プローブ配置を変えた結果からは、配置 A は配置 B よりも平均吸光係数が低かった。これは大脳縦裂付近に存在する脳脊髄液の吸光係数が低いためと考えられる。
また実験と頭部シミュレーションともにプローブ間距離、配置による平均吸光係数の変化は同様の傾向を示しており、被験者の頭部構造モデルは適切に実験での吸収散乱特性を評価できることが示唆された。今後、被験者の頭部モデルに対して様々な特性の光を用いたシミュレーションを行うことで、神経活動に起因する灰白質の吸収散乱特性変化を被験者ごとに適切に評価する光源・実験・解析手法のデザインが可能になると考えられる。

図 5．平均吸光係数のシミュレーションと実測値の比較
エラーバーは実測値の標準偏差。

文 献
1） Hohl, Angela, and Athanasios Gavrielides. "Bifurcation cascade in a semiconductor laser subject to opticalfeedback." Physical review letters 82.6 (1999): 1148. （PMID なし）

52） Kotani K, Yamaguchi I, Ogawa Y, Jimbo Y, Nakao H, Ermentrout GB. Adjoint method provides phase response
functionsfordelay-inducedoscillations.

PhysRevLett.

2012Jul27;109(4):044101.

doi:

10.1103/PhysRevLett.109.044101. Epub 2012 Jul 24. PMID: 23006090
3） Yamaguchi I, Ogawa Y, Jimbo Y, Nakao H, Kotani K. Reduction theories elucidate the origins of complex
biological rhythms generated by interacting delay-induced oscillations. PLoS One. 2011;6(11):e26497. doi:
10.1371/journal.pone.0026497. Epub 2011 Nov 7. PMID: 22087228
4） 庭野恭彰、赤尾旭彦、四本裕子、神保泰彦、小谷潔：時間分解分光法によるヘモグロビン濃度の推定と BOLD 信
号との定量的比較、電気学会論文誌 C, 138 巻, 7 号, p. 790-798, 2018, DOI: 10.1541/ieejeiss.138.790 （PMIDなし）

6