Numerical Simulation of Photoacoustic Effect and Its Possibility of Applications to Diagnostic Imaging and Treatment Support

Abstract

著者らの研究グループでは光音響イメージングに関する基礎研究から臨床研究までを幅広く実施しており，その一部として，光音響波を発生した光吸収体の濃度の定量的イメージングや金ナノ粒子を造影剤として用いる際の光音響波の評価のためにコンピュータを用いたシミュレーション（数値計算）を行っている．この数値計算技術を発展させることで，皮膚科領域での光音響作用による治療の際の治療計画に寄与することができる可能性がある．本稿では光音響波の数値計算とそれを応用した研究を解説する．

Translated Abstract

The authors have been studying the photoacoustic imaging from its basic to clinical applications. As a part of the wide-ranged studies, the numerical simulations of the photoacoustic pressure wave have been carried out for the image reconstruction of the concentration of the photon absorber and for the evaluation of gold nanoparticle as a contrast agent. The numerical simulation can be applied to planning safe and effective treatments using photoacoustic effect in dermatology. This article reviews the numerical simulations and the applications of them.

1.  はじめに

レーザー光を用いた診断・治療技術は光と組織の相互作用を巧みに利用している．本稿ではレーザー光と生体の相互作用の中から，細胞・組織中の色素がパルスレーザー光を吸収することによって応力波が生じる，光音響的（光機械的）作用に注目する．皮膚科領域では，熱緩和時間が短い色素の沈着等を選択的に治療できるナノ秒等の短いパルス時間幅（パルス幅）を持つパルスレーザー光が用いられているが，この時に光熱的作用のみならず光音響的作用が治療効果に深く関与していると考えられている．近年ではピコ秒パルスレーザー光照射で生じる光音響的作用によってより効果的で選択的な刺青治療などが行われている¹⁾．レーザー治療をさらに高度化するために，Andersonは多光子励起蛍光顕微鏡や光音響イメージング等のイメージング技術の発展に言及し，近年はそれらの技術を臨床に応用する好機である²⁾と述べている．

光音響イメージング³⁾では，ナノ秒パルスレーザーを用いて，治療の際よりも非常に小さいフルエンスの光を照射することで，非侵襲的光音響的作用による熱弾性過程を経て生じる光音響波を超音波トランスデューサで検出し，画像再構成を行って光音響波を発した色素の分布を断層画像化する．ヘモグロビンが周辺組織よりも比較的強く吸収する波長530 nm周辺や750 nm程度の近赤外領域のレーザー光を用いたヒトの手掌や小動物などの微小血管や脳活動に伴うヘモグロビン濃度変化のイメージングに関する研究⁴⁾や，臨床において前立腺がん⁵⁾，子宮頚部病変⁶⁾や乳がん⁷⁾の病状が進んだ際に変化する微小血管のイメージングに関する研究が行われている．また，波長532 nmなどのレーザー光によってメラノーマなどのがん細胞のイメージングに関する研究^8,9)があり，その高い空間分解能により注目を集めている．さらに高いコントラストを得るために光を強く吸収する性質を持つ金ナノ粒子等を造影剤として用いる研究も行われている．

近年ではコンピュータの性能向上に伴って，著者らが研究しているイメージング対象の定量化を目指した光音響的作用のシミュレーション（光音響波の発生と伝播の数値計算）が比較的手軽に行えるようになってきた．将来的にシミュレーションをさらに発展させることで，治療支援に応用できる可能性がある．大城ら¹⁾によれば特にピコ秒パルスレーザー光を用いた治療においては至適なレーザー光の照射条件等は現状では十分に明らかにされておらず，今後さらに検討していく必要がある．また，患者毎に条件を最適化するテーラーメイド医療への要望が高まる可能性がある．このような課題に対して，光や応力波の作用のシミュレーションを応用すれば，照射条件等の治療計画の立案の支援のために様々な検討が行える．

そこで本稿では，2節において光音響波の発生と伝播に関する理論と数値計算法について説明する．3節においては著者らの金ナノ粒子が発生する光音響波と散乱が光音響波に与える影響に関する研究を，4節では光音響波の数値計算を用いた色素濃度分布の定量的画像再構成法を紹介する．近年では皮膚科領域で光音響的作用の治療メカニズムが議論されているので，5節では皮膚科領域の治療で見られる光音響的作用について述べ，光音響的作用の数値計算の治療支援応用の可能性について議論し，6節で結言を述べる．

2.  生体内の光と光音響波伝播

細胞・組織中の色素がパルス光を吸収すると熱弾性過程やアブレーション（細胞・組織の飛散），光学破壊（プラズマの発生）によって応力波が発生する¹⁰⁾．これ以降では，熱弾性過程によって生じ，組織を大きく損傷しない弱い応力波を光音響波として区別して呼ぶことにし，本節では光の伝播とその後の光音響波の発生と伝播について，その順番で理論を概説し，数値計算例を示す．（光学破壊で生じる衝撃波は必要に応じて応力波と区別する．）

2.1  生体内の光伝播

生体内において細胞膜と細胞質，あるいは細胞質と細胞小器官の間で屈折率が異なり，また，細胞質内の様々な物質の濃度の分布によって屈折率が異なるため，光は様々な方向に散乱される．生体組織が光を散乱する度合いは散乱係数μ_s [mm⁻¹]（[cm⁻¹]も単位としてよく用いられる，1 cm⁻¹ = 0.1 mm⁻¹），あるいは，換算散乱係数μ_s' = (1 − g)μ_sで表す．g（異方散乱パラメータ）は光が前方に散乱されやすい場合には1，等方的に散乱される場合には0，後方に散乱されやすい場合には−1に近い値になる．生体組織中ではg = 0.9程度である．換算散乱係数や後述の吸収係数については解説論文^11,12)に詳しく紹介がある．文献¹¹⁾では皮膚，脳，乳房，骨，脂肪などについて，波長λ [nm]とμ_s' [mm⁻¹]の関係を近似式μ_s' = a(λ/500)^−bで与えている．ここで，aとbは組織毎に異なる定数である．例えば皮膚について，文献¹⁰⁾中の各種論文データの平均的な値はa = 4.60，b = 1.421である．また，Salomatina らの測定¹³⁾の結果から，真皮ではa = 4.53，b = 1.292，上皮ではa = 6.87，b = 1.161となっている．これによれば，真皮，上皮のμ_s'はそれぞれ，波長530 nmで 4.20 mm⁻¹，6.42 mm⁻¹になる．

一方，血液中のヘモグロビンや水，メラニン，ビリルビン，脂質などの様々な内因性の（生体内部に元々存在する）光吸収体や検査等の造影剤として投与された蛍光物質や金ナノ粒子などの外因性の光吸収体が光を吸収する．光が吸収される度合いは吸収係数μ_a [mm⁻¹]がで表す．光吸収体のμ_aはその物質の持つモル吸光係数ε [mm⁻¹/(mol/l)]とその濃度C [mol/l]の積で表され，組織のμ_aは，様々な光吸収体の組織中に占める体積分率で重みづけした光吸収体のμ_aの総和となる．Svaasand らの研究¹⁴⁾やMeglinskiらの研究¹⁵⁾では血液濃度，血液酸素飽和度，メラニン等の影響を考慮してヒトin vivoにおける波長と真皮，上皮のμ_aの関係を定式化している．波長530 nmにおける真皮，上皮のμ_aは，Svaasand らはそれぞれ0.38 mm⁻¹，1.31 mm⁻¹，Meglinskiらは0.283 mm⁻¹，0.17 mm⁻¹としており，血液濃度の見積もり等によって違いがある¹²⁾．この値はCaucasianについてのものであり，Asian，Africanではまた異なる^14,16)．生体内に多く含まれる水とヘモグロビンは近赤外光（波長700～1,000 nm）を吸収しにくい．この波長帯での生体組織のμ_aはおよそ0.02～0.1 mm⁻¹と比較的小さく，近赤外光は生体深部にまで到達する．μ_s'はおよそ0.5～3.0 mm⁻¹の値をとるとされている¹⁷⁾．

以上の散乱と吸収を受けて光エネルギーは光照射位置から離れるほど減衰する．この光エネルギーの伝播はふく射輸送現象であり，位置rにおいて，ある単位方向ベクトルsの方向に進行する光の放射束密度をψ(r, s) [Wmm⁻²sr⁻¹]として，以下のふく射輸送方程式¹⁸⁾で定式化される，

  

$$s\cdot \nabla \psi =-({\mu }_s+{\mu }_a)\psi +{\mu }_s{\displaystyle {\int }_{4\pi }P(s,s\text{'})\psi d\omega \text{'}+S,}$$（1）

ここで，左辺は放射束密度の変化率，右辺第一項は散乱と吸収による光の減少を表す．第二項は散乱によってs'の方向を向いていた光が確率密度関数（散乱位相関数）P(s, s')に従ってsに方向を変えることによる増加を表している．P(s, s')には前述のgをパラメータとするHenney-Greenstein 関数が良く用いられる．第三項Sは，ここでは光照射強度を表している．（1）式は有限差分法等の方法で数値計算法ができる¹⁹⁾が，現状では膨大なメモリと時間が必要になるため，モンテカルロ法を用いて近似的に解くことが一般的である．

モンテカルロ法を用いた光伝播シミュレーションでは光子パケットと呼ばれるある単位の光エネルギーを持つ粒子を想定して，多数の光子パケットの移動を追跡する²⁰⁾．光子パケットが散乱する方向は散乱位相関数に基づいて，また，光子パケットが次の散乱を受けるまでの自由行程はLambert-Beerの法則をもとに確率的に決定する．光子パケットは通過する位置の吸収係数に応じてエネルギーを減少し，その減少分を記録することでその位置で吸収される光エネルギーを決定する．光子パケットのエネルギーが基準値以下になったら追跡を終了する．モンテカルロ法は並列計算²¹⁾やGPU²²⁾を用いて高速で計算できるようになっている．Prahl，Jaques，Wangら²³⁾や，Fang²⁴⁾のモンテカルロ法による光伝播計算のためのソースが公開されている．

十分に散乱された後の光の伝播は，μ_aがμ_s'よりも十分小さく，等方的に散乱されると仮定して，以下の光拡散方程式で近似される¹⁸⁾，

  

$\left\{-\nabla \cdot D(r)\nabla +{\mu }_a\right\}\Phi (r)=q_0(r),$（2）

ここで，フルエンス率Φ(r) [Wmm⁻²]はψを全立体角にわたって積分した値であり，q₀(r)は光照射強度（等方的であると仮定する）である．（2）式では光は等方的に散乱されるものとする．D(r) = 1/(3 μ_s') [mm]は光拡散係数である．（2）式は任意形状の媒体を四面体要素で分割し，有限要素法等によってベクトルと行列を用いた式に書き直して数値計算できる．

2.2  光音響波の発生と伝播

光音響波はその発生する範囲（立方体領域とする）において，光照射で発生する熱の拡散と光音響波の伝播が無視できる程度に短時間にパルスレーザー光が作用することで効率的に発生する．この条件は熱閉じ込め条件，応力閉じ込め条件と呼ばれ，それぞれ以下の（3）式および（4）式で表される²⁵⁾，

  

$$\text{τ}<{\text{τ}}_{\text{th}}=d_{\text{c}}{}^2/{\alpha }_{\text{T}},$$（3）

  

$$\text{τ}<{\text{τ}}_{\text{s}}=d_{\text{c}}/v,$$（4）

ここで，τ[s]はパルス幅であり，τ_th[s]は熱閉じ込め条件の閾値（熱緩和時間），τ_s[s]は応力閉じ込め条件の閾値（応力緩和時間），α_T[mm²/s]は熱拡散係数（組織ではおおよそ0.14 mm²/s），vは音速（軟組織は水と近い値でおよそ1.5 × 10⁶ mm/s）であり，d_cは光音響波を発生させる立方体領域の代表長さである．例えば100 μm（=0.1 mm）の範囲で光音響波を発生させるのであれば，τ_th = 71 ms，τ_s = 67 nsとなり，光音響イメージングや皮膚科における治療にも用いられるパルス幅10 ns程度のナノ秒パルスレーザーで効率的に光音響波を発生することになる．

光音響波の振幅は内因性，外因性の色素等の吸収係数とそこに到達する光のフルエンス率の積に依存する．位置rに光が到達してΦ(r)が決まると，そこで吸収される光エネルギー密度は，H(r, t) = μ_aΦ(r)δ(t) [Wmm⁻³]で表される．光エネルギーの吸収と熱膨張の過程は光音響波の伝播と比較して極めて短時間で経過するため，デルタ関数を用いている．非粘性均質媒体中の光音響波の伝播は以下の光音響波動方程式³⁾で記述される，

  

$${\nabla }^{2}p(r,t)-\frac{1}{v^2}\frac{{\partial }^2}{\partial t^2}p(r,t)=-\frac{\beta }{C_p}\frac{\partial }{\partial t}H(r,t),$$（5）

ここで，vは生体内の音速，βは熱膨張係数，C_pは比熱である．生体軟組織中の音速は水とほぼ等しく，およそ1.5 mm/μsである．フルエンス率と光音響波の振幅は線形の関係にある．無限媒体中の光音響波は，（5）式を解いて，以下のように求まる，

  

$p(r,t)=\frac{\beta }{4\pi C_p}{\displaystyle {\int }_V\frac{\partial H(r,t-t\text{'})}{\partial t}{\Vert r\Vert }^{-1}}dr,$（6）

ここで，t' = ||r||/vである．以上の方程式から光伝播と光音響波伝播をシミュレーションすることができる．（5），（6）式も光伝播の時と同様に有限差分法や有限要素法を用いて計算できる．

Coxらはk-space法による光音響波伝播計算を提案している²⁶⁾．一様な無限媒体中の光音響波伝播計算において，k-space法はある時刻tの光音響波を以下の3つのステップで計算する：（1）時刻t = 0 μsで光音響波（初期光音響場）p₀(r)のフーリエ変換p₀(k)を計算する，（2）p₀(k)にcos(vkt)を掛け算する，（3）（2）を逆フーリエ変換する．ここで，kは波数である．光音響波は時間と位置によって変化する関数であるが，（1）～（3）の計算手順では，時刻tの光音響波の位置情報を波数情報に変換すると，それが数式の上でp₀(k)とcos(vkt)の単純な掛け算で表すことができ，数値計算が簡便になることを利用している．cos(vkt)は(exp(ivkt) + exp(−ivkt))/2とも書けるので，物理的には光音響波を反対の位相をもつ2つの波に分けて計算していると解釈することもできる²⁷⁾．Coxらはこのk-space法を音速が不均一な媒体における光音響伝播の計算にも拡張している²⁸⁾．k-space法のプログラムk-Waveはウェブサイト²⁹⁾において公開されている．

熱弾性過程で発生する光音響波では媒質が弾性的に変位する（元の位置に戻る）が，アブレーションや光学破壊で発生する応力波では弾性的ではないため，応力波の波形が光音響波とは異なる³⁰⁾（（5）式右辺が異なる）．応力波の発生や到達による媒体の変形・変性が十分に小さければ上記のように計算できるが，変形・変性が無視できない場合は媒体の形状の変化や変性による音響的特性パラメータの変化を反映して計算する必要があると考えられる．波動方程式が使用できる範囲ではフルエンスと応力波の関係は線形である．光音響波・応力波の発生・到達によって生じる細胞・組織の損傷については損傷が起きる閾値等を考慮して別途判断する必要がある．衝撃波が生じる場合¹⁰⁾には流体力学における保存則から導出される非線形偏微分方程式を計算する必要がある．

2.3  イメージングのための光音響波の数値計算の例

光音響波に関する数値計算を行っている研究の例として，YasaらのKubelka-Munk近似による光伝播計算を用いた研究³¹⁾や，Coxらによるk-space法による筋や脂肪，骨など音速が異なる組織で構成された腕のような領域での伝播を計算した研究²⁸⁾，Jacquesによる皮膚層と血管をモデリングした媒体に関するモンテカルロ法を用いた研究³²⁾等がある．

著者らは光拡散方程式（2）と光音響波動方程式（5）を有限要素法で計算した³³⁾．ここでは2次元の光散乱・吸収媒体内に色素の濃度が高い高光吸収領域（一辺1 mmの領域）を3ヶ所に配置し，筋肉組織内に血管があるような媒体を模擬した（Fig.1）．パルスレーザー光の波長は生体イメージングでよく用いられる近赤外領域を仮定した．パルスレーザー光は媒体の上部から照射される．パルスレーザー光の伝播時間は無視できるため，Fig.1の時刻t = 0 μsで初期光音響場が決定している．光音響波の振幅はフルエンス率と吸収係数の積に比例し，光は高光吸収領域以外によっても吸収されるため，フルエンス率が最も大きくなる照射位置近傍において大きい光音響波が発生した．また，高光吸収領域においては，フルエンス率は小さいが，吸収係数が大きいために強い光音響波が発生した．それぞれの高吸収領域は同じ吸収係数であるが，光照射位置に最も近く，フルエンス率が比較的大きい光照射位置直下の高光吸収領域において光音響波が最も大きく，光照射位置から離れてフルエンス率が小さくなるほど光音響波は小さくなる．同じ吸収係数を有していても光照射位置からの距離によって発生する光音響波の振幅に違いが生じる．t = 0 μs以降は発生した光音響波が音速1.5 mm/μsで伝播していく様子が観察された．

Fig.1 

Conditions and the propagation of the photoacoustic waves generated by near-infrared light (the images with t = 0, 1, and 2 μs are reprinted from the literature³²⁾ with permission).

3.  金ナノ粒子が発生する光音響波と散乱の影響

光音響イメージングでは様々な金属等のナノ粒子の中でも特に，強い光吸収特性を有する金ナノ粒子を造影剤として用いた光音響イメージング研究が数多く行われている³⁴⁾．金ナノ粒子をポリエチレングリコールによって表面修飾（PEG化）することや，がん細胞に過剰に発現するタンパク質に対する抗体を結合することなどによって，がんの造影剤として用いる研究が行われている^35,36)．

著者らは金ナノ粒子が発生する光音響波のシミュレーションを試みた³⁷⁾．この研究はより効果的に光音響波を発生する金ナノ粒子の形状や大きさを探索するために有用である．まず初めに，Mie散乱理論や，任意の形状の金ナノ粒子内の電界を計算する離散双極子近似法を用いて，球形状や，Teranishiら³⁸⁾の合成した立方体，八面体形状の金ナノ粒子の吸収断面積（吸収係数）と散乱断面積（散乱係数）を求め，それらをもとにモンテカルロ法を用いて波長532 nmのレーザー光の金ナノ粒子の水溶液中の伝播計算し，そこから発生して超音波トランスデューサで検出される光音響波を（6）式を用いて計算した．光音響波は吸収された光エネルギーに依存し，吸収断面積が大きい一辺が50 nmの立方体形状金ナノ粒子は，一粒子あたりに換算して，直径が50 nmの球状の2倍，直径10 nmの球状金ナノ粒子の350倍の振幅を持つ光音響波を発生することがわかった．この結果は実測に一致した．

また，この研究³⁷⁾では周囲媒体の散乱係数の影響も調べた．Fig.2は直径10 nmの球状金ナノ粒子水溶液の散乱係数が変化した場合の光音響波を示している．金ナノ粒子の濃度は一定である．散乱係数が大きくなるにつれて光音響波の振幅が大きくなっている．この結果は実測ともよく一致した．このような振幅の変化は金ナノ粒子水溶液中の吸収された光エネルギーの分布が散乱係数によって変化することで生じる．Fig.3にはモンテカルロ法で計算した吸収光エネルギー分布を示している．散乱係数が大きくなっていくにつれて光の直進成分が減り，光子が媒体中に広がって光路長が伸びるために吸収される光エネルギーが増えて光音響波の振幅が増加することになる．

Fig.2 

Simulated photoacoustic signals from the spherical gold nanoparticle in Intralipid solution with various μ_s (reprinted from the literature³⁶⁾ with permission).

Fig.3 

The distribution of the light energy absorbed by the aqueous solution of the spherical gold nanoparticle in the cylindrical well with a radius of 4 mm for (a) μ_s = 2.4, (b) 14.3, (c) 47.6, and (d) 142.9 mm⁻¹ (reprinted from the literature³⁶⁾ with permission).

生体は光の散乱体であり，組織の種類や状態によって散乱係数は変化する．上記の研究のように，生体内の光音響波の発生に関して直接の原因である吸収だけではなく，散乱について考慮することで，光音響的作用を用いた診断や治療の精度を向上させることができると考えられる．

4.  光音響の数値計算を応用した光音響イメージング

光音響イメージングでは超音波トランスデューサによって検出した光音響波から色素の分布を画像として再構成する．光音響イメージングでヒトに用いられるパルスレーザー光の照射エネルギーは5節で紹介する治療に用いられるもの（0.1～10 J/cm²程度）とは大きく異なることに注意されたい．光音響イメージングにおけるレーザー光照射はIECやJISなどの各種規格にもとづいて，最大許容露光量（Maximum Permissible Exposure: MPE）を下回る範囲で用いられている³⁾．例えば著者らの臨床研究用装置は，IEC60825-1に従って，照射されるフルエンス（波長750 nm）は14.3 mJ/cm²（MPEは25.2 mJ/cm²）である³⁹⁾．一般的にMPE以下の照射条件であれば皮膚組織に障害は生じないと考えられ，著者らの臨床研究でも障害は起きていない．

光音響イメージング研究の多くでは，画像は（5），（6）式のみをもとにしたbackprojection法^40,41)などのように測定された光音響波の到達時間と信号強度から直接的に画像再構成している．2.3で述べたように，光音響波は光照射位置から離れるほど小さくなり，画像中における輝度は小さくなるため，生体浅部と深部の画像の輝度による状態の比較が困難である．そこで，光音響波の数値計算を用いて，検出された光音響波から生体内の色素濃度（吸収係数）分布を定量的に画像再構成する方法が研究されている⁴²⁾．この方法は実測された光音響波と，数値計算によって求めた光音響波の誤差を最小化する吸収係数などの光学特性値を繰り返し計算によって求める方法である．光拡散方程式を用いた方法^43-46)や，輻射輸送方程式を用いた研究⁴⁷⁾，モンテカルロ法を用いた研究^48,49)が行われている．この方法によれば，生体深部における光の減衰を考慮して色素濃度を計算するため，理想的には生体の浅部と深部を比較できる画像が得られる．

著者らがモンテカルロ法を用いて行ったシミュレーション研究では，照射用光ファイバーを設置した曲率半径25 mmの焦点型超音波トランスデューサ（プローブ）で光音響波を検出しながら計測対象表面を0.2 mm間隔で2次元スキャンした．計測対象内の異なる深さ（1 mm，4.5 mm，8 mm）に同じ吸収係数をもつ血管を模した棒状の高光吸収領域を配置した時（Fig.4）にy軸に沿ってスキャンした時に検出される光音響波データをシミュレーションした結果をFig.5左上に示す．この図では横軸にプローブの移動距離，縦軸に光音響波到達時間を示し，光音響波の振幅をカラースケールで示している．実際にイントラリピッドとインク水溶液を用いて作製した生体の光学特性と音速を模擬したファントムにNd:YAGレーザーで励起した光パラメトリック発振器（OPO）を用いて波長755 nmのナノ秒パルスレーザー光を照射した際に得られたデータ（Fig.5右上）とよく一致している．深さ1 mmにある高光吸収領域と比較して，深さ4.5 mmの高光吸収領域が発生する光音響波の振幅は約0.4倍，深さ8 mmの時には約0.2倍となっており，深さに応じて光音響波の振幅が異なっている．Fig.5左上のシミュレーションで得たデータから再構成した画像（Fig.5下）では深さ1 mmの吸収係数が最も大きくなり，4.5 mm，8 mmのところはその75%程度の吸収係数が再構成され，検出された光音響波データでは異なるように見える高光吸収領域が実は深さによらずおおよそ同じ程度の吸収係数をもつことが分かる．ここで示した結果では深い位置では空間分解能が低下し，吸収係数は浅いところよりもやや低く推定されている．著者らの研究³³⁾では，深い位置にある高吸収領域の範囲が大きく推定されることに伴って吸収係数が低く推定される様子が示されていることから，高吸収領域の体積を正確に推定することが深い位置における吸収係数の改善につながる可能性があり，光音響波測定位置を増やすなどの工夫をして領域を正確に特定できるようにしたり，画像のぼやけを軽減する処理⁴⁵⁾を画像再構成に施したりすることによって空間分解能を改善し，定量性を向上させることができる可能性がある．著者らの方法は光伝播と光音響波伝播の方程式を組み合わせて，光音響波から直接に吸収係数分布を求める方法であるが，超音波画像の再構成法を用いて測定された光音響波から先にその発生領域を求めておき，その情報を元に吸収係数（色素濃度）を推定する方法⁴²⁾を採用すれば，超音波画像再構成法によっては高い空間分解能を得て，定量性を改善できる可能性がある．画像再構成の精度改善や計算時間の削減などが定量的な色素濃度分布の画像再構成の今後の課題である．

Fig.4 

Conditions of the numerical simulation and phantom experiment of the quantitative photoacoustic image reconstruction.

Fig.5 

Measured photoacoustic data in the numerical simulation (top left) and phantom experiment (top right), and the image of μ_a reconstructed from the measured photoacoustic data (simulation) (bottom).

以上のように，未だ上述の課題は残されているが，定量的な光音響画像再構成では不均一な生体内の吸収係数分布を画像として再構成する．がん組織においては血管の新生等により正常な組織とは血液濃度や酸素飽和度が異なることがある^50,51)．また，皮膚の色素沈着等では色素濃度が高くなる．いずれの場合も変化した光吸収体の吸収スペクトルに対応した吸収係数に変化を生じる．このように吸収係数や散乱係数といった光学特性値は組織の状態によって変化するため，光学特性値の定量的な画像再構成は診断に有用な情報を提供できる可能性がある．また，5節で議論する数値計算の治療支援応用においては，事前に数値計算に用いる光学特性値の分布を知っている必要があり，この時にヒトin vivo光音響測定データからの定量的な画像再構成が必要になると考えられる．

5.  皮膚科領域における光音響的作用と数値計算の治療支援応用の可能性

著者らは未着手であるが，数値計算を応用した治療機序の解明や治療支援の可能性がある．本節では皮膚科領域における光音響的作用と刺青治療の数値計算例を紹介し，治療支援応用について展望する．

5.1  皮膚科領域における光音響的作用

皮膚科領域におけるパルスレーザーを用いた治療はAnderson⁵²⁾のselective photothermolysisが基礎になっており，皮膚組織における熱緩和時間よりも十分に短いパルス幅のレーザー光照射によってターゲットとなる色素が含まれている皮膚組織を選択的に治療し，熱的な損傷等の副作用を低減することができる．Watanabe⁵³⁾によればこの時に光が皮膚に引き起こす作用は（1）光化学的作用，（2）光熱的作用，（3）光音響的作用であるとされている．大城ら¹⁾によれば，メラニンの熱緩和時間は50 ns，表皮は3～10 msであり，ナノ秒やピコ秒のような短いパルス幅のレーザーを用いることで2節の（3），（4）式で表された熱および応力閉じ込め条件が満たされ，光音響的作用が生じやすくなっている．

実際にナノ秒パルスレーザー光が光音響的作用による応力波の発生を示す例として，Gómezら⁵⁴⁾の角質層アブレーションによる薬剤の経皮的吸収の促進に関する研究ではQ-switched Nd:YAGレーザー（パルス幅約10 ns，波長1,064，532，355 nm）照射によるウサギの耳の角質層のレーザーアブレーションによって生じる応力波をマイクロフォン（Elecret）によって測定している．アブレーションに必要な最小の照射フルエンスは約6.4 J/cm²（波長1,064 nm），7 J/cm²（532 nm），0.9 J/cm²（355 nm）であり，これを超えた約7.6 J/cm²（波長1,064 nm），26 J/cm²（532 nm），1.2 J/cm²（355 nm）まではフルエンスと応力波の振幅は線形的な関係にあった．この場合には2節の（5）式で計算できると考えられる．以下に紹介する研究では波長1,064 nmでの最小フルエンス約6.4 J/cm²を下回る照射条件であり，アブレーションは起きにくく，フルエンスと応力波の振幅の関係は線形的である可能性が高いが，応力波は照射フルエンスだけでなく色素の濃度（吸収係数）や色素が存在する位置（深さ）にも依存するため実際の現象・作用は複雑で複合的であると考えられる．

Friedmanら⁵⁵⁾はQ-switched Nd:YAGレーザー（パルス幅4～6 ns，波長1,064 nm，平均フルエンス3.4 J/cm²，繰り返し周波数10 Hz，スポットサイズ6 mm）を用いたニキビ痕の“nonablative”な治療を行って6ヶ月間にわたって皮膚性状を測定した．副作用は一時的な紅斑や軽度の点状出血だけであり，色素沈着異常や瘢痕は見られなかった．CO₂レーザーやEr:YAGレーザーによる“ablative”な治療においては術後の紅斑や感染，瘢痕等の副作用も少なくない⁵⁶⁾ことと比較すると，光音響的作用の利用による副作用の軽減が示されていると考えられる．Yaghmaiらのphotoagingの治療に関する研究⁵⁷⁾では，Q-switched Nd:YAGレーザー（パルス幅5～7 ns，波長1,064 nm，平均フルエンス3.2 J/cm²，繰り返し周波数10 Hz，スポットサイズ6 mm）の照射の際に，通常の治療で照射するエネルギーを2回に分けて短時間の間に照射することで，刺すような，あるいは焼けるような感覚を減らすことができ，治療効果は同等か改善するような結果が得られている．この効果は光音響的作用を利用していることによると結論付けている．

近年ではナノ秒よりも短いピコ秒のパルス幅をもつレーザー光による治療も行われるようになってきている．大城らによる総説論文¹⁾では刺青治療や真皮メラノサイトーシス⁵⁸⁾，加齢性色素性病変の治療におけるピコ秒パルスレーザーの臨床応用が紹介されており，ピコ秒パルスレーザーによる治療は，光熱的作用ではなく，光音響的作用が支配的になるために副作用が少なく，これまでとは一線を画した新しい技術として期待されている．一方で，文献上でも照射等に関する治療条件に関しての言及がほとんどなく，波長の影響や最適な照射時間等については今後の課題であるとされている．

5.2  刺青治療とその数値計算例

刺青治療ではルビー（波長694 nm，パルス幅28～40 ns）やAlexandrite（波長755 nm，パルス幅50～100 ns），Nd:YAG（波長1,064 nm，パルス幅5～10 ns）などのQ-switchedナノ秒パルスレーザーが用いられていた⁵⁹⁾．刺青には，カーボン（黒）や酸化クロム（緑），硫化水銀（赤）などの金属粒子やアゾ色素のような有機色素など様々な吸収スペクトルの特徴を持つ物質が用いられるため，レーザーの波長やパルス幅の選定が重要である⁶⁰⁾．色素の粒子の大きさは10～100 nm 程度であり，熱緩和時間は0.1～10 ns（10～100 psとの推定もある）と非常に短いため⁶⁰⁾，最近ではピコ秒パルスレーザーの使用がより効果的であるとされてきている⁵⁹⁾．

Rossら⁶¹⁾のピコ秒とナノ秒Q-switched Nd:YAGレーザー（パルス幅35 psもしくは10 ns，フルエンス0.65 J/cm²，繰り返し周波数10 Hz，スポットサイズ1.4 mm）による刺青治療を比較した研究では，ピコ秒パルスレーザーの方が治療効果が高いことを示した．また，カーボン粒子水溶液を用いた実験ではレーザー照射によっての色が薄くなり，粒子が膨張したことや吸熱性steam-carbon reaction（粒子にひびが入り，中に蒸気が入り込んで，粒子を空洞にして無色化する）が起きたことを確認した．また，ピコ秒パルスレーザー光照射では光学破壊（10⁹～10¹¹ W/cm²のフルエンスで発生するのが典型的，単位に注意，W = J/s）によって衝撃波を生じる可能性があり，線維芽細胞死や色素粒子の変化のほかに真皮の炎症等を引き起こす可能性があると述べている．

Hoらの数値計算による研究⁶²⁾ではパルスレーザー光によって色素粒子が破壊される機序を明らかにした．シミュレーションはレーザー光と組織の相互作用を計算するプログラムであるLATIS⁶³⁾を用いている．Mie散乱の理論から刺青の色素（グラファイト粒子）の吸収断面積を求め，臨床で用いられるようなレーザー光照射条件（例えば波長1,064 nm，平均フルエンス0.27 J/cm²，パルス幅10～35 psなど）を用いて，表皮や真皮，脂肪層などの光学特性値を設定，皮膚表面から2 mmの深さに存在する粒子（直径10～100 nm）が吸収する光エネルギーをモンテカルロ法で計算し，粒子の温度上昇，粒子内に発生する応力波と，それによって生じる粒子にかかる抗張力などを計算した結果，粒子の温度が融点に達しない照射条件でも十分大きい抗張力が発生し，粒子が破壊されることが明らかになった．このことはすなわち，刺青色素の粒子は熱的作用ではなく，光音響的作用によって破壊されることを意味している．また，粒子の温度が水の沸点を超えることで周辺組織にキャビテーションが発生し，損傷する可能性が示されている．さらに，粒子の温度上昇の結果，steam-carbon reactionが起きている可能性があるとしている．様々な条件を検討した結果，組織の損傷をできるだけ小さくして治療効果を得るためにレーザー光のパルス幅はおおよそ10～100 psが良いと結論付けている．

5.3  数値計算の治療支援応用への展望

以上のHoらの研究⁶²⁾や光線力学療法の治療計画等にモンテカルロ法を用いた光伝播のシミュレーションを用いている研究⁶⁴⁾のように，本稿で紹介した生体内の光伝播や光音響波の伝播に関する数値計算法を発展させて様々な波長，皮膚組織の光学特性，照射フルエンス等の条件で計算することで，治療計画の策定等の治療支援が可能になると考えられる．このときに計算によって得られる光照射が引き起こす作用に関する知見が，実際の臨床の状況を正確に反映・予測することを担保するために，臨床研究における知見やファントムを用いた実験と比較して，計算精度を慎重に確認する必要がある．また，光照射と応力の発生・伝播に伴って生じる細胞・組織の変化についても数値計算を行う際には光学特性値に加えて熱伝導率や熱膨張係数等の熱力学的特性値，生体組織の引張強度，破壊強度等の力学的特性など，多数のパラメータを必要とする．より高精度に計算するために，光学特性値測定についても引き続き研究が行われることが期待される．皮膚や筋肉，骨などの組織の力学特性の測定⁶⁵⁾や細胞への応力波の影響¹⁰⁾はこれまでにも様々な研究で行われているが，光ピンセット⁶⁶⁾などの技術を用いた細胞レベルでの測定によってもさらに高い精度で測定できる可能性がある．光学特性値分布をヒトin vivoで測定・イメージングできる定量的光音響画像再構成や拡散光トモグラフィー^17-19)のような技術を応用したり，光コヒーレンスエラストグラフィー⁶⁷⁾のような技術でin vivoでの力学特性の測定が高い精度で行えるようになれば，将来的には患者ごとに測定したパラメータを用いた数値計算によって治療効果を予測できるようになり，より高度なテーラーメイド医療が実現できる可能性がある．

6.  おわりに

計算機の性能が引き続き向上しているなかで，今後のさらなる計算時間削減や，高い計算精度の実現によって，数値計算が臨床現場でも役立つ技術になる可能性がある．光音響的作用をはじめとしたレーザー光の生体への作用を利用した診断・治療に関わる数値計算技術が，臨床におけるニーズを理解しながら研究開発が進められることによって，臨床研究や基礎医学・生物学の研究，医工学研究などの関連分野と足並みをそろえて発展して，医療の質の向上に貢献し，医師や患者の役に立つ技術となることを期待する．

利益相反の開示

本論文に関して，開示すべき利益相反はない．

謝辞

本論文で紹介した研究の一部は，日本医療研究開発機構　産学共創基礎基盤研究プログラム（ヒト生体イメージングを目指した革新的バイオフォトニクス技術の構築），科研費（16H02155, 18K04188）および，京都大学化学研究所共同利用・共同研究（#2015-42, 2016-44, 2017-44）の補助を受けて行われた．

引用文献

• 1)   大城  貴史， 大城  俊夫， 佐々木  克己．ピコ秒レーザーによる刺青および良性色素性病変に対する治療．日レ医誌．2017; 37: 427-434.
• 2)   Anderson  RR. Lasers for Dermatology and Skin Biology. J Invest Dermatol. 2013; 133: E21-E23.
• 3)   Li  C,  Wang  LV. Photoacoustic tomography and sensing in biomedicine. Phys Med Biol. 2009; 54: R59-R97.
• 4)   Wang  LV,  Hu  S. Photoacoustic Tomography: In Vivo Imaging from Organelles to Organs. Science. 2012; 335: 1458-1462.
• 5)   Horiguchi  A,  Shinchi  M,  Nakamura  A,  Wada  T,  Ito  K,  Asano  T,  Shinmoto  H,  Tsuda  H,  Ishihara  M. Pilot study of prostate cancer angiogenesis imaging using a photoacoustic imaging system. Urology. 2017; 108: 212-219.
• 6)   Okawa  S,  Sei  K,  Hirasawa  T,  Irisawa  K,  Hirota  K,  Wada  T,  Kushibiki  T,  Furuya  K,  Ishihara  M. In vivo photoacoustic imaging of uterine cervical lesion and its image processing based on light propagation in biological medium. Proc SPIE. 2017; 10064: 100642S.
• 7)   Toi  M,  Asao  Y,  Matsumoto  Y,  Sekiguchi  H,  Yoshikawa  A,  Takada  M,  Kataoka  M,  Endo  T,  Kawaguchi-Sakita  N,  Kawashima  M,  Fakhrejahani  E,  Kanao  S,  Yamaga  I,  Nakayama  Y,  Tokiwa  M,  Torii  M,  Yagi  T,  Sakurai  T,  Togashi  K,  Shiina  T. Visualization of tumor-related blood vessels in human breast by photoacoustic imaging system with a hemispherical detector array. Sci Rep. 2017; 7: 41970.
• 8)   Wang  Y,  Xu  D,  Yang  S,  Xing  D. Toward in vivo biopsy of melanoma based on photoacoustic and ultrasound dual imaging with an integrated detector. Biomed Opt Express. 2016; 7: 279–286.
• 9)   Jathoul  AP,  Laufer  J,  Ogunlade  O,  Treeby  B,  Cox  B,  Zhang  E,  Johnson  P,  Pizzey  AR,  Philip  B,  Marafioti  T,  Lythgoe  MF,  Pedley  RB,  Pule  MA,  Beard  P. Deep in vivo photoacoustic imaging of mammalian tissues using a tyrosinase-based genetic reporter. Nat Photon. 2015; 9: 239-246.
• 10)   Doukas  AG,  Flotte  TJ. Physical characteristics and biological effects of laser-induced stress waves. Ultrasound Med Biol. 1996; 22: 151-164.
• 11)   Jacques  SL. Optical properties of biological tissues: a review. Phys Med Biol. 2013; 58: R37-R61.
• 12)   Lister  T,  Wright  PA,  Chappell  PH. Optical properties of human skin. J Biomed Opt. 2012; 17: 090901.
• 13)   Salomatina  E,  Jiang  B,  Novak  J,  Yaroslavsky  AN. Optical properties of normal and cancerous human skin in the visible and near-infrared spectral range. J Biomed Opt. 2006; 11: 064026.
• 14)   Svaasand  LO,  Norvang  LT,  Fiskerstrand  EJ,  Stopps  EKS,  Berns  MW,  Nelson  JS. Tissue parameters determining the visual appearance of normal and port-wine stain. Lasers Med Sci. 1995; 10: 55-65.
• 15)   Meglinski  IV,  Matcher  SJ. Quantitative assessment of skin layers absorption and skin reflectance spectra simulations in the visible and near-infrared spectral region. Physiol Meas. 2002; 23: 741-753.
• 16)   Tseng  SH,  Grant  A,  Durkin  AJ. In vivo determination of skin near-infrared optical properties using diffuse optical spectroscopy. J Biomed Opt. 2008; 13: 014016.
• 17)   Yamada  Y,  Okawa  S. Diffuse Optical tomography: Present status and its Future. Opt Rev. 2014; 21: 185-205.
• 18)   Arridge  SR. Optical tomography in medical imaging. Inverse Problems. 1999; 15: R41-R93.
• 19)   Klose  AD,  Ntziachristos  V,  Hielscher  AH. The inverse problem based on the radiative transfer equation in optical molecular imaging. J Compt Phys. 2005; 202: 323-345.
• 20)   Wang  LH,  Jacques  SL,  Zheng  LQ. MCML—Monte Carlo modeling of photon transport in multi-layered tissues. Comput Meth Prog Bio. 1995; 47: 131-146.
• 21)   Alerstam  E,  Svensson  T,  Andersson-Engels  S. Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration. J Biomed Opt. 2008; 13: 060504.
• 22)   Alerstam  E,  Lo  WCY,  Han  TD,  Rose  J,  Andersson-Engels  S,  Lilge  L. Next-generation acceleration and code optimization for light transport in turbid media using GPUs. Biomed Opt Express. 2010; 1: 658-675.
• 23)  Prahl S. Monte Carlo Light Scattering Programs.（参照2018年9月30日），https://omlc.org/software/mc/.
• 24)  Fang Q. MCX open-source light transport simulator.（参照2018年9月30日），http://mcx.sourceforge.net/mmc.
• 25)   Zhou  T,  Yao  J,  Wang  LV. Tutorial on photoacoustic tomography. J Biomed Opt. 2016; 21: 061007.
• 26)   Cox  BT,  Beard  PC. Fast calculation of pulse photoacoustic fields in fluids using k-space methods. J Acoust Soc Am. 2005; 117: 3616-3627.
• 27)  Cox B, Beard PC. Modeling photoacoustic propagation in tissue using k-Space techniques. Wang LV. (ed.) Photoacoustic Imaging and Spectroscopy. Florida: CRC Press; 2009. p. 25-34.
• 28)   Cox  BT,  Kara  S,  Arridge  SR,  Beard  PC. k-space propagation models for acoustically heterogeneous media: Application to biomedical photoacoustic. J Acoust Soc Am. 2007; 121: 3453-3464.
• 29)  Photoacoustic Imaging Group, Department of Medical Physics and Bioengineering, University College London. Photoacoustic propagation models.（参照2018年9月30日），http://www.medphys.ucl.ac.uk/research/mle/photoacoustic.htm.
• 30)  石原 美弥．バイオメディカルフォトニクス－生体医用光学の基礎と応用―電気学会次世代バイオメディカル・レーザ応用技術調査専門委員会編．1.3節，光の生体組織への作用．東京：オーム社；2009. p. 22-30.
• 31)   Yasa  ZA,  Jackson  WB,  Amer  NM. Photothermal spectroscopy of scattering media. Appl Opt. 1982; 21: 21-31.
• 32)   Jacques  SL. Coupling 3D Monte Carlo light transport in optically heterogeneous tissues to photoacoustic signal generation. Photacoustics. 2014; 2: 137-142.
• 33)   Okawa  S,  Hirasawa  T,  Kushibiki  T,  Ishihara  M. Numerical evaluation of linearized image reconstruction based on finite element method for biomedical photoacoustic imaging. Opt Rev. 2013; 20: 442-451.
• 34)   Wang  S,  Lin  J,  Wang  T,  Chen  X,  Huang  P. Recent advances in photoacoustic imaging for deep-tissue biomedical applications. Theranostics. 2016; 6: 2394-2413.
• 35)   Mallidi  S,  Kim  S,  Karpiouk  A,  Joshi  PP,  Sokolov  K,  Emilianov  S. Visualization of molecular composition and functionality of cancer cells using nanoparticle-augmented ultrasound-guided photoacoustics. Photoacoustics. 2015; 3: 26-34.
• 36)   Zeng  Q,  Iwakuma  N,  Sharma  P,  Moudgil  BM,  Wu  C,  MacNeill  J,  Jiang  H,  Grobmyer  SR. Gold nanoparitcles as a contrast agent for in vivo tumor imaging with photoacoustic tomography. Nanotechnology. 2009; 20: 395102.
• 37)   Okawa  S,  Hirasawa  T,  Sato  R,  Kushibiki  T,  Ishihara  M,  Teranishi  T. Numerical and experimental investigations of dependence of photoacoustic signals from gold nanoparticles on the optical properties. Opt Rev. 2018; 25: 365-374.
• 38)   Eguchi  M,  Mitsui  D,  Wu  HL,  Sato  R,  Teranishi  T. Simple reductant concentration-dependent shape control of polyhedral gold nanoparticles and their plasmonic properties. Langmuir. 2012; 28: 9021-9026.
• 39)   Irisawa  K,  Hirota  K,  Hashimoto  A,  Murakoshi  D,  Ishii  H,  Tada  T,  Wada  T,  Hayakawa  T. Photoacoustic imaging system for peripheral small-vessel imaging based on ultrasound technology. Proc SPIE. 2016; 9708: 970807.
• 40)   Sperl  JI,  Zell  K,  Menzenbach  P,  Haisch  C,  Ketzer  S,  Marquart  M,  Koening  H,  Vogel  MW. Photoacoustic image reconstruction: a quantitative analysis. Proc SPIE. 2007; 6631: 663103.
• 41)   Xu  M,  Wang  LV. Universal back-projection algorithm for photoacoustic computed tomography. Phys Rev E. 2005; 71: 016706.
• 42)   Cox  B,  Laufer  JG,  Arridge  SR,  Beard  PC. Quantitative spectroscopic photoacoustic imaging: a review. J Biomed Opt. 2012; 17: 061202.
• 43)   Laufer  J,  Cox  B,  Zhang  E,  Beard  P. Quantitative determination of chromophore concentrations from 2D photoacoustic images using a nonlinear model-based inversion scheme. Appl Opt. 2010; 49: 1219-1233.
• 44)   Yuan  Z,  Wang  Q,  Jiang  H. Reconstruction of optical absorption coefficient map of heterogeneous media by photoacoustic tomography coupled with diffusion equation based regularized Newton method. Opt Express. 2007; 15: 18076-18081.
• 45)   Okawa  S,  Hirasawa  T,  Kushibiki  T,  Ishihara  M. Reconstruction of the optical absorption coefficient from photoacoustic signals measured by scanning coaxial probe with regularization methods. Proc SPIE. 2014; 8943: 89433W.
• 46)   Okawa  S,  Hirasawa  T,  Kushibiki  T,  Ishihara  M. Effect of the approximations of light propagation on quantitative photoacoustic tomography using two-dimensional photon diffusion equation and linearization. Opt Rev. 2017; 24: 705-726.
• 47)   Saratoon  T,  Tarvainen  T,  Cox  BT,  Arridge  SR. A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation. Inverse Problems. 2013; 29: 075006.
• 48)   Fonseca  M,  Malone  E,  Lucka  F,  Ellwood  R,  An  L,  Arridge  S,  Beard  P,  Cox  B. Three-dimensional photoacoustic imaging and inversion for accurate quantification of chromophore distribution. Proc SPIE. 2017; 10064: 1006415.
• 49)   Okawa  S,  Hirasawa  T,  Tsujita  K,  Kushibiki  T,  Ishihara  M. 3D quantitative photoacoustic image reconstruction using Monte Carlo method and linearization. Proc SPIE. 2018; 10494: 104944Y.
• 50)   Vaupel  P,  Kallinowski  F,  Okunieff  P. Blood flow, oxgen and nutrient supply, and metabolic microenvironment of human tumors: a review. Cancer Res. 1989; 49: 6449-6465.
• 51)   Sorg  BS,  Moeller  BJ,  Donovan  O,  Cao  Y,  Dewhirst  MW. Hyperspectral imaging of hemoglobin saturation in tumor microvasculature and tumor hypoxia development. J Biomed Opt. 2005; 10: 044004.
• 52)   Anderson  RR,  Parrish  JA. Selective Photothermolysis: precise microsurgery by selective absorption of pulse radiation. Science. 1983; 220: 524-527.
• 53)   Watanabe  S. Basic of laser application to dermatology. Arch Dermatol Res. 2008; 300: S21-S30.
• 54)   Gómez  C,  Costela  Á,  García-Moreno  I,  Llanes  F,  Teijón  JM,  Blanco  D. Laser treatment on skin enhancing and controlling transdermal delivery of 5-Fluorouracil. Lasers Surg Med. 2008; 40: 6-12.
• 55)   Friedman  PM,  Jih  MH,  Skover  GR,  Payonk  GS,  Kimyai-Asadi  A,  Geronemus  RG. Treatment of atrophic facial acne scars with the 1064-nm Q-Switched Nd:YAG laser. Arch Dermatol. 2004; 140: 1337-1341.
• 56)   Asilian  A,  Salimi  E,  Faghihi  G,  Dehghani  F,  Tajmirriahi  N,  Hosseini  SM. Comparison of Q-Switched 1064-nm Nd: YAG laser and fractional CO2 laser efficacies on improvement of atrophic facial acne scar. J Res Med Sci. 2011; 16: 1189–1195.
• 57)   Yaghmai  D,  Garden  JM,  Bakus  AD,  Gold  MH,  Saal  BM,  Goldberg  DJ,  Massa  MC. Photodamage therapy using an electro-optic Q-Switched Nd:YAG laser. Lasers Surg Med. 2010; 42: 699-705.
• 58)   Ohshiro  T,  Ohshiro  T,  Sasaki  K,  Kishi  K. Picosecond pulse duration laser treatment for dermal melanocytosis in Asians: A retrospective review. Laser Ther. 2016; 25: 99-104.
• 59)   Naga  LI,  Alster  TS. Laser Tattoo Removal: An Update. Am J Clin Dermatol. 2017; 18: 59-65.
• 60)   Barua  S. Laser-Tissue Interaction in tattoo removal by Q-Switched lasers. J Cutan Aesthet Surg. 2015; 8: 5-8.
• 61)   Ross  EV,  Naseef  G,  Lin  C,  Kelly  M,  Michaud  N,  Flotte  TJ,  Raythen  J,  Anderson  RR. Comparison of responses of tattoos to picosecond and nanosecond Q-switched Neodymium:YAG lasers. Arch Dermatol. 1998; 134: 167-171.
• 62)  Ho DD-M,  London  R,  Zimmerman  GB,  Young  DA. Laser-tattoo removal—A study of the mechanism and the optimal treatment strategy via computer simulations. Lasers Surg Med. 2002; 30: 389-397.
• 63)   London  RA,  Glinsky  ME,  Zimmerman  GB,  Bailey  DS,  Eder  DC,  Jacques  SL. Laser-tissue interaction modeling with LATIS. Appl Opt. 1997; 26: 9068-9074.
• 64)   Davidson  S,  Weersink  R,  Haider  M,  Gertner  M,  Bogaards  A,  Giewercer  D,  Scherz  A,  Sherar  M,  Elhilali  M,  Chinet  J. Treatment planning and dose analysis for interstitial photodynamic therapy of prostate cancer. Phys Med Biol. 2009; 54: 2293–2313.
• 65)   岡  久雄， 入江  隆．皮膚軟部組織モデル化のための生体定数．バイオメカニズム学会誌．1993; 17: 216-222.
• 66)   Ou-Yang  HD,  Wei  MT. Complex fluids: Probing mechanical properties of biological systems with optical tweezers. Annu Rev Phys Chem. 2010; 61: 421-440.
• 67)   Liang  X. Biomechanical properties of in vivo human skin from dynamic optical coherence elastography. IEEE Trans Biomed Eng. 2010; 57: 953-059.