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

194. ヒト頸部の光と光音響波の伝播解明と画像診断応用
藤井 宏之
北海道大学 大学院工学研究院 機械・宇宙航空工学部門 熱流体物理学研究室
Key words：近赤外光，光音響イメージングと医療応用，生体模擬溶液，初期圧力，光散乱・熱弾性

緒 言
600～1,100 nm の波長帯域である近赤外光は、他の波長帯域の光に比べて生体による吸収が小さく、生体深部に伝播することが可能であり、生体イメージングに利用されている［1］。光イメージングでは、光の吸収係数と酸素化・脱酸素化ヘモグロビンやメラニンなどの化学組成と関係に基づいて、がん部位の推定や脳機能イメージングなどに応用されている。光イメージングの研究課題として、光の散乱過程により、イメージングのコントラストや解像度が著しく低下することが挙げられる。コントラストや解像度の向上に向けて、様々な試みが行われている。その中で、光音響イメージングが、近年目覚ましく発展している。光音響イメージングでは、近赤外のレーザーを照射し、熱膨張と収縮（熱弾性）を経て発生する光音響波を検出する。光音響イメージングは、光イメージングの利点である高いコントラストと音響イメージングの利点である高い解像度の両方を有している［2］。
光音響イメージングの中でも、生体内部を対象とした光音響トモグラフィが開発されている。光音響トモグラフ
ィでは、光音響波の伝播を記述する数理モデルである光音響波動方程式を用いる。トランスデューサーより検出された圧力と光音響波動方程式による計算値の差が小さくなるように、生体内部の初期圧力や光吸収量の空間分布を再構成する。光吸収量は光吸収係数と関係しているため、生体内部の酸素飽和度を取得することができる。
従来のトモグラフィは生体浅部を対象としているが、生体深部の情報を定量的に評価することを目的とした、定量的光音響トモグラフィの開発が進められている［3］。光吸収量は吸収係数と光強度の積であり、同じ吸収係数の値であっても、光強度が生体浅部と深部では異なるため、光吸収量が異なる。従って、光吸収量から吸収係数と光強度を分離することが定量的評価では必要である。分離のためには、光強度を計算する必要があり、光伝播を記述する数理モデルである光伝播モデルが定量的光音響トモグラフィに導入される。
光音響トモグラフィは様々な生体部位のイメージングに利用されているが、ヒト頸部の甲状腺腫瘍診断に近年応用されている［4］。甲状腺腫瘍の非侵襲診断では主に音響イメージングが使用されている。音響イメージング
では構造特性の違いを基準として、腫瘍の良性と悪性を判別するが、構造特性の変化が小さい腫瘍の種類もあり、その場合の検出は困難になる。他方、光音響トモグラフィではヘモグロビン濃度の違いを基準として判別するため、音響イメージングでは難しいケースを補うことが期待される。先行研究では、音響イメージングに比べて、光音響トモグラフィのイメージング結果の向上が報告されている。そのため、定量的光音響トモグラフィを甲状腺腫瘍診断に適用すれば、更なるイメージング結果の向上が期待されるが、ヒト頸部における光と圧力波の伝播モデルの構築が難しく、未だ適用されていない。ヒト頸部前方では、空洞領域である気管があり、光は気管境界面や頸部外部境界面で反射・屈折し、圧力波に影響を及ぼすと考えられる。しかし、圧力波に対する影響がどの程度であるのかについては、未だ明らかではない。圧力波に対する影響を調べることは、光と圧力に関する計算精度と計算効率に優れた伝播モデルを構築するために重要である。光の反射・屈折をモデルに組み込むことは、計算コストが高い。そのため、その影響が小さければ、計算モデルに考慮しないことが望ましい。これら背景を踏まえて、ヒト頸部モデルや立方体媒体における、光の反射・屈折が光伝播や圧力伝播に及ぼす影響を明らかにすることを本研究の目的の一つとした。
生体深部を対象とした定量的光音響トモグラフィにおいて、初期圧力における熱弾性、光散乱、光伝播による1

寄与を明らかにすることが重要である。各物理過程の寄与が明らかになれば、初期圧力分布の再構成結果から、光吸収係数を抽出することが可能となる。生体は非常に複雑な媒体で個体差も大きいため、生体を模擬した溶液であるコロイド溶液が基礎研究で使用されている［5，6］。コロイド溶液は、粒子の体積分率や粒径を調整することで、熱弾性、光散乱、光伝播を生体の平均値に調整することができる。また、再現性が高い。先行研究では、光伝播モデルとして輻射輸送論、圧力伝播モデルとして光音響波動論が適用されているが、これらの理論では光散乱特性や熱弾性特性がパラメータであり、実験値が使用されている。しかし、特性を計測するための多くの実験が必要であり、体系的な理解は容易ではない。光散乱特性や熱弾性特性を分子モデルより計算し、特性間の関係性や特性と光・圧力伝播の関係性を包括的に明らかにすることは基礎研究として重要である。これら背景を踏まえて、コロイド溶液に対して、光散乱が及ぼす初期圧力への影響を定量的に明らかにすることを本研究のもう一つの目的とした。

方 法
1．ヒト頸部モデリング
ヒト頸部に対する光伝播モデルについて、以下 3 項目を適切に考慮した：（A）空間不均一性、（B）光の異方性、（C）頸部外部境界及び気管境界における光の反射・屈折。項目（A）について、研究代表者らは、成人女性
頸部 MR 画像（図 1a）に対して、解析ソフト ImageJ を用いて、6 つの領域に分割した（図 1b）：背景組織（筋肉、脂肪、骨の平均媒体）、気管、脊髄、脊椎、血管、甲状腺。気管は空洞領域であり、光は殆ど散乱や吸収されず、直進する。一方で、気管以外の部分は光を強く散乱させ、血管部分では散乱に加えて吸収も強い。項目（B）
について、気管部分は等方的に散乱されるが、他の領域では強く前方に散乱される。項目（C）について、外部境
界と気管境界では、空気と接しており、生体と空気では屈折率が異なる（生体の屈折率は水に近く 1.4 程度、空気は 1）。従って、光は境界面で反射・屈折する。甲状腺がんのモデル化として、新生血管より、がん領域ではヘモグロビン濃度が高い場合を考え、光吸収係数を他の領域より高く、1.0/cm に設定した。

図 1．ヒト頸部 MR 画像とヒト頸部モデルa）成人女性ヒト頸部MR画像（スケールは図bに記載）。
b）ヒト頸部モデル（MR画像からの領域抽出）。

ヒト頸部に対する圧力伝播モデルについて、熱弾性特性が必要となる。気管以外の各部位における熱弾性特性は部位依存が小さいことから、気管では空気の特性値、気管以外の部位では水の特性値を用いた。また、ヒト頸部モデルに対して、メッシュ生成ソフト Gmsh を用いて三角形メッシュを作製した。
2．光伝播モデル
生体部位の巨視的なスケールでは、光をエネルギー粒子としての輻射（ふくしゃ）と捉えることができる。この光伝播を記述する方程式が輻射輸送方程式となる。輻射輸送方程式では境界面における光の反射・屈折をフレ2ネル則より考慮する。輻射輸送方程式の計算負荷は膨大であり、計算負荷低減のために、輸送方程式を拡散近似した、光拡散方程式が広く使用されている。光拡散方程式では光の反射・屈折を拡散反射より考慮した。本研究では、輻射輸送方程式と光拡散方程式を同じ媒体で計算し、光の反射・屈折が光伝播に及ぼす影響を解析した。
ヒト頸部モデルに対しては、Streamline Upwind/Petrov Galerkin 法（高精度の有限要素法）を適用した。立方体媒体に対しては、3 次精度風上差分法と Delta-Eddington-Galerkin 求積法（研究代表者らが開発した方法）を用いた。
3．光音響波の圧力伝播モデリング
光音響波の圧力伝播は光音響波動方程式によって記述される。University College London の Cox らの研究グループによって開発されている計算パッケージ、k-Wave［7］を用いて、圧力伝播を数値計算した。ここでは、熱と応力の閉じ込め条件が満たされていると仮定した。これにより、初期圧力の発生は、各媒体微小要素より局所的に発生することを仮定した。
4．コロイド溶液における初期圧力モデリング
光散乱と熱弾性をサブナノスケールよりモデル化した［8］。光散乱については電磁波理論の一つ、干渉散乱理論を適用した［9，10］。熱弾性については従来の溶媒相と粒子相の二相モデルでは考慮されていなかった、粒子間相互作用を組み込んだモデルを開発した。光散乱特性の結果を光伝播モデルに組み込んだ。サブナノスケール
では、コロイド溶液を球形粒子が分散している媒体として、ミリメートルスケールでは、溶液を連続体（立方体）としてモデル化した。

結果および考察
1．ヒト頸部における光音響波の圧力伝播
光吸収体が甲状腺領域に存在する場合における、輻射輸送方程式と光拡散方程式に基づいた光音響波動方程式の数値計算結果を図 2 に示す。レーザーを照射した時刻 0 マイクロ秒において、光源と光吸収体から生じる圧力が他の領域に比べて高くなった。
他方、光源や吸収体から離れた領域においての初期圧力は非常に小さくなった。
これは、光が指数関数的に減衰して内部に伝播するためである。初期圧力は積分光強度と吸収係数の積に比例するため、生体深部では積分光強度の減衰に伴って非常に小さくなる。光音響波は、パルス的に、各点より発生し、時間が経つにつれて重なり合い、境界面に設置された検出点に到達する。時刻 0.71 マイクロ秒において、光源からの圧力波と吸収体からの圧力波が重なり合った。時刻 1.41 マイクロ秒から 4.95 マイクロ秒にかけて、気管に到達した圧力波の位相が反転した。これは気管と背景組織との密度差によって固定端反射が生じたことによっており、検出点には、この逆位相の圧力波が遅れて到達することになる。従って、光吸収体による圧力と逆位相の圧力の到達時間に差がある方が分離する上で望ましいと考えられた。光伝播モデルの差について、明確な違いは観察されなかった。輻射輸送方程式に比べて、光拡散方程式の計算負荷は小さいが、光源近傍の光伝播を適切に記述することができない。しかし、吸収体によって生じる圧力波のピークが鋭ければ、そのモデルの差による影響は小さいため、光拡散方程式の適用は、実質上問題ないと考えられた。

3

図 2．ヒト頸部モデルにおける光音響波の圧力伝播、最大値で規格化a） 光伝播モデルとして輻射輸送方程式を使用した場合。
b） 光伝播モデルとして光拡散方程式を使用した場合。

2．3 次元媒体における光伝播の数値解析
2 次元ヒト頸部モデルでは、輻射輸送方程式と光拡散方程式の計算結果に大きな差異はなかったが、3 次元媒体でも同様な結果が得られるのかどうか検証する必要がある。3 次元ヒト頸部モデルへ向けた初期段階として、3 次元立方体媒体に対して、輻射輸送方程式と光拡散方程式の数値計算結果を比較した。両計算モデルの相対誤差は、媒体の散乱係数で異なるが、全体的に 5％以下であった。よって、3 次元媒体においても光拡散方程式の適用が可能であると考えられた。
3．生体模擬溶液における光音響波の圧力伝播
光散乱特性について、干渉散乱理論と独立散乱理論の結果を比較し、体積分率が 5％以上では、粒子群の相互作用が粒子群の平均配置や光散乱に強く影響を及ぼすことがわかった。両理論の結果の相対誤差は 15％に達することがわかった。熱弾性特性について、粒子群の相互作用を考慮しない従来モデル、研究代表者らが開発した相互作用を考慮したモデル、実験値と比較し、3 者が殆ど一致する結果を得た。この結果より、熱弾性特性には粒子群の相互作用は殆ど影響を及ぼさないことが明らかとなった。
溶液内部に存在する吸収体から発生する初期圧力における、熱弾性と光散乱・光伝播の寄与を数値解析した。
図 3 に、シリカ溶液における初期圧力の体積分率依存性を示す。熱弾性による寄与Γと光散乱・光伝播による寄与Φを区別している。初期圧力において、熱弾性特性による影響は限定的であり、初期圧力の 90％以上が、光散乱・光伝播による寄与であることが明らかとなった。この結果より、熱弾性のモデル化に比べて、光散乱のモデル化が重要であることが示唆された。

4

図 3．溶液における初期圧力 Pa の体積分率依存性
体積分率 1％の結果より規格化、初期圧力は熱弾性による寄与Γと光散乱・光伝播による寄与Φの積で与えられる。

4．今後の課題
3 次元ヒト頸部モデルへの適用は研究期間に実行することができなかった。研究代表者が開発した計算手法を3 次元立方体媒体に適用し、従来法よりも計算精度と効率を各段に向上させることに成功したが、3 次元ヒト頸部モデルに適用するためには、更なる計算手法の発展が必要であることを認識した。境界面における光の反射・屈折が光伝播や圧力伝播に及ぼす影響は小さいことを踏まえて、光伝播モデルとして光拡散方程式を適用することが有効であると考えられた。コロイド溶液に対する初期圧力マルチスケールモデルについて、今回は単波長で考察した。今後は、多波長を考え、初期圧力の波長依存性を解析することが課題である。

共同研究者・謝辞
本研究の共同研究者は、浜松医科大学光尖端医学教育研究センター生体医用光学研究室の星詳子教授である。
本研究の遂行にあたり、多大なご支援を頂きました上原記念生命科学財団に心より感謝申し上げます。

文 献1） Ntziachristos V. Going deeper than microscopy: the optical imaging frontier in biology. Nat Methods. 2010
Aug;7(8):603-14. doi: 10.1038/nmeth.1483. Epub 2010 Jul 30. PMID: 20676081.
2） Wang LV, Hu S. Photoacoustic tomography: in vivo imaging from organelles to organs. Science. 2012 Mar23;335(6075):1458-62. doi: 10.1126/science.1216210. PMID: 22442475.
3） Cox B, Laufer JG, Arridge SR, Beard PC. Quantitative spectroscopic photoacoustic imaging: a review. J
Biomed Opt. 2012 Jun;17(6):061202. doi: 10.1117/1.JBO.17.6.061202. PMID: 22734732.
4） Dima A, Ntziachristos V. In-vivo handheld optoacoustic tomography of the human thyroid.
Photoacoustics. 2016 Jun 27;4(2):65-69. doi: 10.1016/j.pacs.2016.05.003. PMID: 27766210.
5） Petrova E, Liopo A, Oraevsky AA, Ermilov SA. Temperature-dependent optoacoustic response andtransient through zero Grüneisen parameter in optically contrasted media. Photoacoustics. 2017 Jun
23;7:36-46. doi: 10.1016/j.pacs.2017.06.002. PMID: 28725558.

56） Pogue BW, Patterson MS. Review of tissue simulating phantoms for optical spectroscopy, imaging anddosimetry. J Biomed Opt. 2006 Jul-Aug;11(4):041102. doi: 10.1117/1.2335429. PMID: 16965130.
7） Treeby BE, Cox BT. k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic
wave fields. J Biomed Opt. 2010 Mar-Apr;15(2):021314. doi: 10.1117/1.3360308. PMID: 20459236.
8） Fujii H, Terabayashi I, Kobayashi K, Watanabe M. Modeling photoacoustic pressure generation incolloidal suspensions at different volume fractions based on a multi-scale approach. Photoacoustics. 2022
May 14;27:100368. doi: 10.1016/j.pacs.2022.100368. PMID: 35646589.
9） Fujii H, Tsang L, Zhu J, Nomura K, Kobayashi K, Watanabe M. Photon transport model for dense
polydisperse colloidal suspensions using the radiative transfer equation combined with the dependent
scattering theory. Opt Express. 2020 Jul 20;28(15):22962-22977. doi: 10.1364/OE.398582. PMID:
32752548.
10）Fujii H, Ueno M, Inoue Y, Aoki T, Kobayashi K, Watanabe M. Model equations of light scattering
properties and a characteristic time of light propagation for polydisperse colloidal suspensions atdifferent volume fractions. Opt Express. 2022 Jan 31;30(3):3538-3552. doi: 10.1364/OE.447334. PMID:
35209609.

6