JP2018108284A - 情報処理装置、情報処理方法およびコンピュータプログラム - Google Patents

Abstract

【課題】生体に照射する複数の放射線の照射量を高速に算出する。【解決手段】本発明の実施形態としての情報処理装置は、生体に照射する複数の放射線の照射量を表す複数の変数を含む目的関数を、変数ごとの制約条件を満たすように最適化する探索処理部を備える。目的関数は、生体における複数の位置のそれぞれに付与される線量と、複数の位置のそれぞれの目標線量との差異を定量化した値を定めており、複数の位置にそれぞれ付与される線量は、位置ごとに定義された複数の係数と、位置ごとに定義された、複数の変数を含む生物学的効果比とを、複数の変数に乗じて定義され、探索処理部は、現在の暫定解において目的関数のニュートン方向へ複数の変数の値を変更していくことにより、制約条件に違反する変数を含むアクティブ集合を生成し、アクティブ集合に基づいて、解の探索方向を決定し、決定した探索方向に向けて次の改善解を算出する。【選択図】図８

Description

本発明の実施形態は、情報処理装置、情報処理方法およびコンピュータプログラムに関する。

高精度外部照射による放射線治療として、放射線の一種である重粒子線を用いた３次元スポットスキャニング法がある。３次元スポットスキャニング法は、患者の身体的な負担が少ないがん治療として注目され始めている。３次元スポットスキャニング法では、照射部位の内部および付近で三次元的に多数のスポットを設定し、スポットごとに細いビームを照射することで、複雑な形状の線量分布を得ることができる。３次元スポットスキャニング法では、実際の治療を行う前に、ユーザー（医師や医学物理士）による医学的所見などをもとに、治療計画を行う。

治療計画では、ユーザーが生体の照射部位（腫瘍部など）に対する線量処方（目標線量分布）を決定し、線量処方に基づいて、コンピューターを用いて、各スポットの照射量などを算出する。照射量の算出は、例えば、目的関数を最適化することにより行う。目的関数は、例えば照射部位の線量分布と、線量処方との二乗誤差を計算する関数である。目的関数の値を最小化することで、各スポットの照射量を算出する。

算出された解（照射量）から照射部位で吸収される線量の分布を計算し、ユーザーが線量処方に照らして、この分布を承認できない場合は（例えばある箇所の線量が許容値を超えるなど）、線量処方を調整して、再度同様のプロセス（治療計画）を行う。医学的所見をもとに目標線量分布を決定しているのに対し、治療計画では放射線の生物物理的な性質なども考慮しており、このため、ユーザーにより承認できる線量分布が得られるまで、治療計画を複数回行うのが一般的である。よって、治療計画には高速性が重要となる。

目的関数の最適化を高速に行うために、目的関数の関数形の情報を利用して、解を逐次的に更新していくＬ−ＢＦＧＳ−Ｂ法が知られている。しかし、この方法では、目的関数の関数形によっては、更新回数が多くなってしまい、高速化が十分に達成されない場合もある。

特許第5454989号

Ehrgott, M., et al., "Mathematical optimization in intensity modulated radio therapy," 4OR (2008) Byrd, R., H., et al., "A limited memory algorithm for bound constrained optimization," Northwestern university technical report NAM-08 (1994) 加瀬優紀，「組織等価比例計数管を用いた重粒子線ビームの線質測定と生物効果の推定（第6章）」，東京工業大学学位論文，2006 Inaniwa, T., et al., "A robust algorithm of intensity modulated proton therapy for critical tissue sparing and target coverage," Phys. med. biol. 56 (2011)

本発明の実施形態は、生体に照射する複数の放射線の照射量を高速に算出する情報処理装置、情報処理方法およびコンピュータプログラムを提供する。

本発明の実施形態としての情報処理装置は、生体に照射する複数の放射線の照射量を表す複数の変数を含む目的関数を、前記複数の変数のそれぞれの制約条件を満たすように最適化することにより、前記複数の変数の値を決定する探索処理部を備える。前記目的関数は、前記生体における複数の位置のそれぞれに付与される線量と、前記複数の位置のそれぞれの目標線量との差異を定量化した値を定めており、前記複数の位置にそれぞれ付与される線量は、前記位置ごとに定義された複数の係数と、前記位置ごとに定義された、前記複数の変数を含む生物学的効果比とを、前記複数の変数に乗じて定義され、前記探索処理部は、現在の暫定解において前記目的関数のニュートン方向へ前記複数の変数の値を変更していくことにより、前記制約条件に違反する変数を含むアクティブ集合を生成し、前記アクティブ集合に基づいて、解の探索方向を決定し、決定した探索方向に向けて次の改善解を算出する。

本発明の実施形態に従った情報処理装置のブロック図。 生物学的効果比の説明図。 ３次元スポットスキャニング法の説明図。 初期解を逐次的に更新して暫定解を得ることを繰り返す様子を示す図。 暫定解におけるニュートン方向と最急降下方向のイメージを示す図。 アクティブ集合生成ステップの動作のイメージ図。 本実施形態に係る改変Ｌ−ＢＦＧＳ−Ｂ法のイメージを模式的に示す。 探索処理部の動作のフローチャートである。 強度変調放射線治療の説明図。 飽和補正線量平均比エネルギーとＲＢＥの関係を示す図。 本実施形態に係る情報処理装置の第１の変形例を示す図。 本実施形態に係る情報処理装置の第２の変形例を示す図。 本実施形態に係る情報処理装置のハードウェア構成例を示す図。

以下、図面を参照しながら、本発明の実施形態について説明する。

本実施形態は、３次元スポットスキャニング法など、放射線（たとえば重粒子線、陽子線）を用いた高精度外部照射による放射線治療において、ユーザー（医師や医学物理士）による治療計画を高速に行うことを支援する。３次元スポットスキャニング法では、照射部位（腫瘍部などの治療対象部位）の内部および付近で三次元的に多数のスポットを設定し、スポットごとに細い放射線（ビーム）を照射することで、複雑な形状の線量分布を得ることができる。

３次元スポットスキャニング法などの治療計画では、ユーザーが、ＣＴ画像や医学的所見をもとに、生体の照射部位に対する線量処方を決定する。そして、決定した線量処方に基づいて、コンピューターを用いて、目的関数を最適化することにより、各スポットの照射量を算出する。

線量処方は、照射部位の情報と、照射部位に対して付与する線量分布の目標値（目標線量分布）とを含む。具体的に、目標線量分布は、照射部位を直方体のメッシュで分割した複数の３次元セル（以下、セルまたは位置と呼ぶ）ごとの線量の目標値を含む。目的関数は、照射部位の線量分布と、線量処方との差異を定量化した値を定めており、一例として、照射部位の線量分布と、線量処方との二乗誤差を計算する関数である。目的関数を最適化することで、各スポットの照射量を解として得る。

算出した各スポットに照射する照射量から、照射部位に吸収される線量の分布（照射部位の線量分布）を計算し、計算した分布が、線量処方に照らして、ユーザーが承認できない場合は、線量処方を再調整して、同様のプロセス（治療計画）を行う。

ここで、放射線の線種としては、重粒子線や陽子線などの粒子線、およびエックス線やガンマ線などの光線がある。重粒子線には、炭素線や酸素線などの種類がある。以下の説明では、放射線として、主に重粒子線、その中でも特に炭素線を想定するが、陽子線など他の種類の粒子線でもよい。また、粒子線でなく、光線（特にエックス線以外の光線）でも可能である。

本実施形態では、各スポットへの放射線の照射量照射量から、照射部位の線量分布を算出するにあたり、生物学的効果比（ＲＢＥ：ｒｅｌａｔｉｖｅ ｂｉｏｌｏｇｉｃａｌ ｅｆｆｅｃｔｉｖｅｎｅｓｓ）を考慮する。以下、生物学的効果比について説明する。

放射線の種類により、同じ線量を照射しても（同じ線量が吸収されても）細胞の生存率が異なることが知られている。このことを考慮するため、たとえば日本国内における重粒子線の治療計画では、物理線量に生物学的効果比（ＲＢＥ）をかけてエックス線換算した値である臨床線量に基づいて治療計画を行っている。

ＲＢＥは、等しい細胞生存率を生じる２種類の放射線の線量の比として定義され、治療計画で用いる場合は、一方を標準的なエックス線（基準エックス線）の線量とする。これを式で表すと以下のようになる。
ＲＢＥ＝Ｄ_ｓｔ／Ｄ_ｒａｄ
Ｄ_ｓｔ：基準エックス線によって目的とする生存率を起こすのに必要な吸収線量
Ｄ_ｒａｄ：対象とする放射線によって同じ生存率を起こすのに必要な吸収線量

ＲＢＥは、生体の体内の深さに応じて変化する。このことを、図２を用いて説明する。エックス線などの光線は、体表面付近で線量が最も大きく、体内での深さが深くなるほど、線量は少なくなる。一方、重粒子線などの粒子線を用いる場合、加速および方向を制御して、その終点（粒子が止まるところ）を制御する。このように制御された粒子線では、体内に入ってから線量は緩やかに上昇し、その終点付近で線量は急上昇し、終点でピークになる。終点より深くなると、線量は急速に低下する。

ＲＢＥの最大値は、重粒子線の線量がピークになる体内深さ（図の位置Ｐ１）に対応し、そのときのＲＢＥの値は約３である。すなわち、このとき臨床線量は、位置Ｐ１における重粒子線の線量の約３倍である。したがって、重粒子線の場合、位置Ｐ１では、エックス線の約１／３の線量で、エックス線と同じ効果（生存率）が得られる。一方、体表面から、位置Ｐ１で急上昇する前までは、ＲＢＥは約２である。

治療計画では、各スポットへの放射線の照射量から、照射部位の線量分布を算出する。算出する線量分布は、生物効果モデルに基づくエックス線換算の線量分布である。この線量分布を線量処方に一致または近づけるよう、各スポットの照射量を逆計算する。これは、線量分布と線量処方（目標線量分布）との二乗誤差を目的関数として最小化する非線形連続最適化問題に定式化することができる（非特許文献１）。この場合の目的関数の変数は、各スポットの照射量である。算出された線量分布が線量処方と完全に一致することはまれであり、ユーザーは自身が承認可能な線量分布となるまで、線量処方を見直しながら治療計画を繰り返し行う。

線量分布（エックス線換算の線量分布）と、線量処方（目標線量分布）との二乗誤差を表す目的関数の式の例を以下に示す。Ｒ×Ｄ×ｘが線量分布、ｐが線量処方に対応する。

ここで、変数ベクトルｘは各スポットへの照射量（すなわち複数の放射線のそれぞれの照射量）を表すｎ個の非負実数値からなるベクトルである。スポットの個数をｎ、各スポットへの照射量を変数ｘ_１，ｘ_２，…，ｘ_ｎとすると、
である。

スポットの数ｎは照射部位の大きさ次第で１０，０００を超える場合もある。３次元スポットスキャニングの概略構成を図３に示す。イオン源から粒子が放射され、放射された粒子が加速器によって加速され、スキャニング電磁石で方向が、スポットに向けて制御される。これにより、放射線がスポットに照射される。粒子の加速および方向を制御することで、照射先のスポットを切り換えることができる。

定数ｐは、照射部位を直方体のメッシュで区切ったｍ個のセル（位置）ごとに定められた目標線量の値からなるベクトルである。セルの個数をｍ、各セルの目標線量をｐ_１，ｐ_２，…，ｐ_ｍとすると、
である。

定数Ｄは、ｎ個のスポットの照射量から、照射部位中のｍか所のセルごとに吸収される線量を算出するために用いられるｍ×ｎ行列である。行列Ｄは、線量付与行列と呼ばれる。行列Ｄを、下記のように表す。

線量付与行列のｊ番目の縦列を、
と記述すると、これは、ｊ番目のスポットの照射量ｘ_ｊに対する線量分布（ｍか所のセルに付与される線量）に対応する。

Ｒは、照射部位中のｍか所の各セルのＲＢＥを表すｍ×ｍ対角行列である。Ｒは、以下のように表すことができる。ｒ_１（ｘ）〜ｒ_ｍ（ｘ）は、セル１〜ｍのＲＢＥを表す。

ＲＢＥを推定するために、生物効果モデルと呼ばれるモデル式が用いられる。たとえば日本国内の重粒子線治療では改良マイクロドシメトリックキネティックモデル（ＭＫＭ）が使われている（特許文献１）。ＭＫＭは、放射線によってドメインと呼ばれる微小体積内に付与されるエネルギー（比エネルギー）から細胞生存率を予測するモデルである。

生物効果モデルを用いて、変数ベクトルｘ＝（ｘ_１，ｘ_２，…，ｘ_ｎ）^Tと、定数（モデルパラメーター）とを含む、ＲＢＥの関数ｒ_１（ｘ）〜ｒ_ｍ（ｘ）を定義する。モデルパラメーターは、使用するモデルおよび放射線の照射条件に基づき、セル毎に事前に決定しておく。

ｒ_ｉ（ｘ）＝Ｄ_ｓｔ／Ｄ_ｉである。ｒ_ｉ（ｘ）の具体例を以下に示す。

ｘ_ｊ（ｊ＝１，２，...，ｎ）以外のすべてのパラメーターは予め決定されている。
Ｓは、ＲＢＥを計算する際に２種類の放射線について合わせる基準となる細胞生存率を表し、日本国内の重粒子線治療では０．１が用いられている。
Ｄ_ｓｔは、標準的なエックス線（基準エックス線）で生存率Ｓを達成するための線量であり、治療計画に先立ち、細胞照射実験などで定まる。
α_０は、放射線の種類（線種）に依存するパラメーターであり、治療計画に先立ち、細胞照射実験などで定まる。
βは、ＭＫＭに固有の定数である。
ｄ_ｉ,ｊは、線量付与行列の（ｉ，ｊ）成分であり、スポットｊに照射された放射線（ビーム）が単位あたりの照射量で位置ｉに付与する線量を表し、予めペンシルビーム法やモンテカルロ法などの既存の計算法で算出される。
は、飽和補正線量平均比エネルギーと呼ばれ、ｄ_ｉ,ｊと併せて算出される。

式（２）の導出方法については、後述する。

ここで、照射部位の線量分布の各成分（照射部位中のｍか所のセルの線量）を
とすると、
式（１）は、以下のように表すこともできる。

目的関数は、式（１）に限定されるものではない。たとえば非照射部位に対する線量制限を考慮した関数や、照射量ｘに対する複数の関数の加重和を用いることもできる。

ここで、目的関数に対する制約条件として、変数ベクトルｘの各変数（各スポットの照射量）が０以上であることを定義する。具体的に、ｊ番目のスポットの照射量をｘ_ｊと表したとき、制約条件は、式（３）のとおり表される。

式（３）の制約条件を満たしつつ、式（１）を最小化するように、変数ベクトルｘを求めることで、ｎ個のスポットに対応するｎ個の照射量を算出する。本実施形態は、このような大規模な非線形連続最適化問題の解を高速に得る手法を提供する。

図１は、本発明の実施形態に従った情報処理装置のブロック図である。図１の情報処理装置は、線量処方データ記憶部１と、パラメーター記憶部２と、目的関数生成部３と、探索処理部４と、照射量データ記憶部５とを備える。探索処理部４は、探索方向算出部１１と、改善解算出部１２とを備える。

線量処方データ記憶部１は、ユーザーにより生成された線量処方のデータを記憶する。線量処方データは、照射部位の情報と、照射部位に対して付与する目標線量分布とを含む。目標線量分布は、複数（ｍ個）のセルにそれぞれ付与する線量の目標値を含む。線量処方データは、式（１）における定数（ベクトル）ｐに相当する。

パラメーター記憶部２は、式（１）における行列（線量付与行列）Ｄ、および行列（ＲＢＥ行列）Ｒ等の情報を記憶している。ここでは行列の形で記憶しているが、これに限定されない。行列ＤおよびＲを構成するのに必要な情報を記憶する限り、任意の形で情報を記憶すればよい。

前述したように線量付与行列Ｄは、ｎ個のスポットの照射量から、照射部位中のｍか所のセルごとに付与される線量を算出するために用いられるｍ×ｎ行列である。線量付与行列Ｄの各成分は、予め定めた定数である。行列Ｒは、照射部位中のｍか所のセルのＲＢＥを表すｍ×ｍ対角行列である。各セルのＲＢＥ（すなわち各対角成分）は、ｘの関数である。

使用する線種および、放射線の照射設定情報に基づき、事前に行列Ｄおよび行列Ｒ等を決定して、パラメーター記憶部２に記憶させておく。照射設定情報は、たとえば重粒子線による３次元スポットスキャニング法の場合、照射方向や、加速器パラメーター、スポットサイズ、線量分布の計算メッシュサイズなどを含む。パラメーターの算出には、ペンシルビーム法またはモンテカルロ法などの公知の算出方法を用いることができるため、ここで述べる事項は最小限に止める。

目的関数生成部３は、線量処方データ記憶部１から線量処方データｐを読み出し、パラメーター記憶部２から行列Ｄと行列Ｒを読み出して、目的関数を生成する。ここでは、式（１）の目的関数を生成する。また、目的関数生成部３は、目的関数の各変数ｘ_ｊ（ｊ＝１，２，…，ｎ）に対する制約条件を生成する。ここでは、式（３）の制約条件を生成する。制約条件のデータをパラメーター記憶部２に格納しておき、目的関数生成部３は、ここから制約条件を取得してもよい。あるいは、目的関数から変数を特定し、変数の各成分が０以上であることを示す制約条件を自動的に生成してもよい。

探索処理部は、制約条件を満たすように、目的関数を最小化することで、解（変数ｘ）を算出する。解は、図４に示すように、予め与えた初期解を逐次的に更新して暫定解を得ることを繰り返すことで、最終的な解（出力解）を得る。初期解（すなわちｘ_１〜ｘ_ｎの初期値）は任意でよい。

本実施形態は、解の算出方法として、目的関数の偏微分係数を利用した解の算出方法であるＬ−ＢＦＧＳ−Ｂ法（非特許文献２）を改変した方法を用いる。以下では、本実施形態に係る方法を、“改変Ｌ−ＢＦＧＳ−Ｂ法”と呼ぶ。改変Ｌ−ＢＦＧＳ−Ｂ法を説明する前に、Ｌ−ＢＦＧＳ−Ｂ法を説明する。

Ｌ−ＢＦＧＳ−Ｂ法は、準ニュートン法をベースにしている。準ニュートン法は、下記のステップＡ１およびステップＡ２を繰り返すことにより、解を逐次的に更新していく算出方法である。

ステップＡ１：目的関数の二次近似関数を擬似的に生成し、生成した二次近似関数の最小解の方向（ニュートン方向）を、探索方向として算出する。最小解とは、関数の値が最小になるときの変数の値である。

ステップＡ２：探索方向（ニュートン方向）に基づいて改善解（現在の解よりも目的関数値が小さくなる解）を算出する。

目的関数の二次近似関数とは、現在の解（初期解または暫定解）
における目的関数値
，一次偏微分係数（ベクトル）∇ｆおよび二次偏微分係数（行列）∇^２ｆにより
と定義される二次関数のことである。なお、目的関数値、一次および二次の偏微分係数が与えられたとき、当該一次および二次の偏微分係数を有する二次関数は一意に決まる。簡単な例として、変数がｘ_１とｘ_２の２つの場合、目的関数値はｆであり、一次微分係数は、∂ｆ／ｄｘ_１， ∂ｆ／∂ｘ_２であり、二次微分係数は、∂^２ｆ／∂ｘ_１ ^２，∂^２ｆ／ｄｘ_２ ^２，∂^２ｆ／∂ｘ_１ｄｘ_２，∂^２ｆ／∂ｘ_２ｄｘ_１の４つであり、これらの目的関数値および一次および二次の偏微分係数を有する二次関数は一意に決まる。ただしＬ−ＢＦＧＳ−Ｂ法では∇^２ｆを計算せず、疑似的に生成された行列Ｂで置き換えている。この行列Ｂは正定値であることを仮定する。

探索方向における改善解の算出には、ラインサーチと呼ばれる公知の算出方法を用いることができる。たとえば、非特許文献２に記載されたバックトラッキング法を用いることができる。これにより、目的関数の関数形が、二次関数に近いとき、少ない更新回数で良好な解を得ることができる。暫定解におけるニュートン方向を模式的に図５に示す。図に記載されている最急降下方向は後述する。

Ｌ−ＢＦＧＳ−Ｂ法は、準ニュートン法を、式（３）のような、変数の下限制約、上限制約またはこれらの両方がある問題への適用も可能としたものである。

Ｌ−ＢＦＧＳ−Ｂ法では、下記のステップＢ１、Ｂ２を繰り返し行う。

ステップＢ１：現在の暫定解において目的関数の傾きである最急降下方向に進んでいった場合（当該方向に沿って各変数の値を変更していった場合）に、制約条件に違反する（式（３）の制約条件の場合、値がゼロ未満となる）変数からなるアクティブ集合を生成する（アクティブ集合生成ステップ）。前述した図５に、暫定解における最急降下方向（目的関数の傾き）を模式的に示す。

ステップＢ２：アクティブ集合に含まれる変数は制約条件を違反しない値に固定して（すなわちアクティブ集合に含まれる変数を定数と見なして）、現在の暫定解における二次近似関数の最小解を算出する。最小解の方向を探索方向として、改善解を算出する。

これにより、制約条件を違反せずに、目的関数値が低下する可能性の高い探索方向を算出する。

ここで、アクティブ集合生成ステップ（ステップＢ１）の詳細を説明する。図６にアクティブ集合生成ステップの動作を説明するための図を示す。アクティブ集合生成ステップは、以下のステップＣ１〜Ｃ３を含む。

ステップＣ１：まず、現在の解（暫定解）ｘ^（０）から開始して、最急降下方向に進む。このとき、最初にいずれかの制約条件が等号で成り立つ解がｘ^（１）であり、かつ、ｊ_１番目の制約条件が等号で成り立ったものとする。２つの解ｘ^（０）、ｘ^（１）を両端とする線分をｓ_１と定める。

ステップＣ２：続いて解ｘ^（１）から開始して、ｊ_１番目の変数の値は（ゼロに）固定したまま、最急降下方向に進む。最初にいずれかの制約条件が等号で成り立つ解がｘ^（２）であり、かつｊ_２番目の制約条件が等号で成り立ったものとする。２つの解ｘ^（１）、ｘ^（２）を両端とする線分をｓ_２と定め、同様にして線分ｓ_３からｓ_ｎを定義する。

ステップＣ３：その上でｓ_１からｓ_ｎを順番に調べ、線分上に二次近似関数の極小解（線分上をいずれの方向に進んでも目的関数の値が大きくなる解）が存在する線分を特定する。このとき、最初に見つかる線分をｓ_ｋとする。そして、ｊ_１からｊ_ｋ−１番目の制約条件に含まれる変数の集合をアクティブ集合として出力する。

なお、Ｌ−ＢＦＧＳ−Ｂ法のさらに詳細なステップは、非特許文献２に記載されているため、ここで述べる事項は最小限に留める。

本実施形態に係る改変Ｌ−ＢＦＧＳ−Ｂ法は、上述したＬ−ＢＦＧＳ−Ｂ法において、ステップＣ１で現在の暫定解から進む方向として、最急降下方向ではなく、ニュートン方向を用いることを特徴とする。

改変Ｌ−ＢＦＧＳ−Ｂ法を説明するための図を図７に示す。図６の最急降下方向が、ニュートン方向に代わっている。ニュートン方向は、前述したように、現在の暫定解において、目的関数の二次近似関数の最小解の方向である。最急降下方向の代わりに、ニュートン方向を用いることで、高速な解の探索が可能になることを、本発明者は知得した。これは、放射線として重粒子線を用いた場合には、目的関数の関数形が二次関数に近いと考えられること、またＬ−ＢＦＧＳ−Ｂ法または改変Ｌ−ＢＦＧＳ−Ｂ法により決定される探索方向はニュートン方向に近くなるためであると考えられる。

本発明者らは、放射線として、重粒子線（炭素線）を用いて、改変Ｌ−ＢＦＧＳ−Ｂ法と、Ｌ−ＢＦＧＳ−Ｂ法を比較するシミュレーションを行った。ただし、照射対象は均質な水の円柱であることを仮定した。また、線量分布を計算するために必要なパラメーター（線量付与行列、ＲＢＥのモデル式等）は、ビームの方向について特許文献１を参考に作成し、ビームの広がりを正規分布で表現した模擬データを用いた。治療計画問題の問題規模は、スポット数ｎが約３３，０００、セル数ｍが約１４，０００であり、ｎの方が大きいのは、３方向から照射する設定で計算したことによる。非特許文献４を参考にして定めた目標精度まで計算するためにかかった時間は、Ｌ−ＢＦＧＳ−Ｂ法では２８０秒、改変Ｌ−ＢＦＧＳ−Ｂ法では８０秒であり、大きく改善された。また、反復回数（解の更新回数）は、Ｌ−ＢＦＧＳ−Ｂ法が１２５回、改変Ｌ−ＢＦＧＳ−Ｂ法は３３回であり、反復回数が約１／４になり、改変Ｌ−ＢＦＧＳ−Ｂ法によって反復回数が大きく減少した。これらにより、高速な解の探索が可能になったことが確認された。

探索処理部４は、探索方向算出部１１と、改善解算出部１２を用いて、制約条件を満たすように、目的関数の解を探索する。探索方向算出部１１は、改変Ｌ−ＢＦＧＳ−Ｂ法に従って、探索方向を算出する。改善解算出部１２は、算出された探索方向において、改善解を算出する。

図８は、探索処理部４の動作のフローチャートである。まず、探索処理部４は、目的関数生成部３から、目的関数のデータと、制約条件のデータを取得する。

探索方向算出部１１は、目的関数の変数ｘ＝（ｘ_１，ｘ_２，…，ｘ_ｎ）^Ｔについて現在の暫定解（最初は任意に決めた初期解）での二次近似関数を疑似的に生成し、その最小解の方向であるニュートン方向を算出する（Ｓ１０１）。目的関数の二次近似関数は、上述した式（４）によって定義される関数である。

探索方向算出部１１は、ステップＳ１０１で決定した方向へ各変数の値を更新していった場合に制約条件に違反する（式（３）の制約条件の場合、値がゼロ未満となる）変数からなるアクティブ集合を生成する（Ｓ１０２）。アクティブ集合を生成する処理の詳細は、前述した通り、ステップＣ１〜Ｃ３に示される。

探索方向算出部１１は、アクティブ集合に基づき、解の探索方向を決定し、決定した探索方向に探索を行うことで、次の改善解を見つける。より詳細には、ステップＳ１０３、Ｓ１０４を行う。

ステップＳ１０３では、アクティブ集合に含まれる変数は制約条件を違反しない値（ここではゼロ）に固定して、現在の暫定解における二次近似関数の最小解を算出し、最小解の方向を探索方向として決定する（Ｓ１０３）。

ステップＳ１０４では、改善解算出部１２は、ステップＳ１０３で決定した探索方向において、改善解（現在の解よりも目的関数値が小さくなる解）を算出する。改善解の算出には、前述したように、ラインサーチと呼ばれる公知の算出方法を用いればよい。

探索処理部４は、終了判定を行い（Ｓ１０５）、処理を継続する場合は、ステップＳ１０４で算出した解を、暫定解として、ステップＳ１０１に戻る。処理を終了する場合は、ステップＳ１０４で算出した解を最終解として出力する。照射量データ記憶部５は、探索処理部４から出力された最終解を受け取り、受け取った最終解を、各スポットの照射量データとして内部に記憶する。

終了判定の例として、解の値、または、解の更新回数、または、改善解の算出までに要した累積時間の少なくとも１つを用いて、処理を終了するかを判定する。例えば、現在の改善解と、直前の改善解からの目的関数値の減少幅（変動幅）を算出し、減少幅が閾値以下であれば、処理の終了を決定し、閾値より大きければ、処理を継続してもよい。放射線治療で用いられるＤ９５などの線量分布の指標を求め、指標が基準値に達していれば処理の終了を決定し、達していなければ処理を継続してもよい。また、解の更新回数が所定値に達した場合は、処理の終了を決定し、所定値に達しない場合は、処理を継続してもよい。または、累積時間が所定値に達した場合は、処理の終了を決定し、達していない場合は、処理を継続してもよい。ここで述べた例のいくつかを組み合わせてもよい。ここで述べた以外の判定方法を用いてもよい。

上述した動作により、制約条件を違反せずに、目的関数値が低下する可能性の高い探索方向を算出することができる。

なお、探索方向算出部１１は、初期解を複数生成してもかまわない。この場合、探索処理部４は、複数の初期解に基づいて算出されたそれぞれの改善解のうち少なくとも目的関数値が最も小さい１つを出力する。

上述した実施形態では、放射線として粒子線を想定したが、また、粒子線でなく、ガンマ線などの光線（特にエックス線以外の光線）でも可能である。エックス線以外の光線であれば、ＲＢＥが１と異なる値になるため、本実施形態の効果を得ることが期待できる。光線の場合、粒子線とは異なり、マルチリーフコメリータを用いて、強度変調放射線治療を行う。この概要を、図９に示す。図９（Ａ）には、マルチリーフコメリータの例が示される。マルチリーフコメリータは、複数のリーフ（金属板）を備え、各リーフを機械的に駆動することで、任意の形状の照射領域を形成する。このマルチリーフコメリータを、図９（Ｂ）に示すように、照射対象に対して回転させながら、回転の各位置で、照射領域を調整しつつ、照射を行う。前述した粒子線の場合、各スポットの照射量を計算したが、光線の場合は、マルチリーフコメリータの全照射領域を複数のピクセルに分割し、回転における各位置におけるピクセルごとの照射量を算出する。

以下では、照射量ｘ_ｊ（ｊ＝１，２，…，ｎ）から位置（セル）ｉ におけるＲＢＥを求める計算式（２）の導出について詳細に説明する。

前述したように、ＭＫＭは生物学的効果比（ＲＢＥ）の計算のために用いられる生物効果モデルの一種である。他の生物効果モデルとして、たとえば欧州で用いられるローカルエフェクトモデルを挙げることができるが（非特許文献３）、今日において日本国内の重粒子線治療ではＭＫＭが用いられていると考えてよい。ＭＫＭによりＲＢＥを計算する方法の例が、特許文献１に開示されている。以下では、特許文献１に基づき、使用する記号（および語句）を本実施形態に合わせて変更しつつ、ＭＫＭを用いたＲＢＥの計算法を示す。

特許文献１の式（１２）に示されるように、着目位置ｉでの細胞生存率Ｓｉは、以下の式（５）で表される。

式（５）を以下のように式変形して、Ｄｉの式（６）が得られる（なお、特許文献１の式（１３）は式変形上の誤りであると考えられる）。

線量Ｄ_ｉはゼロ未満になり得ないことから、式（７）が得られる。

ここで、
である。（Ｚ^＊ _{１Ｄｍｉｘ}）_ｉは、混合照射場中の着目位置ｉにおける飽和補正線量平均比エネルギーの線量加重平均値である。

特許文献１の式（２２）に記載されているように、
である。

したがって、ＲＢＥの計算式は、式（８）のようになる。

以上のとおり、ＲＢＥを求めるための計算式では、細胞照射実験により決定された値が用いられていることから、計算結果には臨床上（実際の治療）の生物学的効果とのかい離がある。そのため臨床上では、経験的な係数をかけて補正した値がＲＢＥとして用いられる。仮にこの係数を１．４５とした場合、最終的にＲＢＥの計算式は、式（９）になる。この式は、先に示した式（２）と同じである。

式（９）の数値的な確認（妥当性の確認）を行うと、以下の通りである。

単一のビーム（スポット）だけを照射する場合、ＲＢＥの計算式は、以下のように、照射量ｘに依存しない値をとる。

特許文献１の図３によれば、飽和補正線量平均比エネルギー（Ｚ^＊ _１Ｄ）_ｉ,１は概ね２０以下の値をとり、単一のビームでは、特に線量のピーク付近において、飽和補正線量平均比エネルギーも高くなることがわかる。単一のビームにおけるＲＢＥを表す式（１０）を、（Ｚ^＊ _１Ｄ）_ｉ,１の関数式と見なし、飽和補正線量平均比エネルギーとＲＢＥの関係をプロットすると、図１０のようになる。この際、他のパラメーターは次のとおり設定した。すなわち、Ｄ_ｓｔの値は特許文献１の段落０１５２より４．０４と定めた。α_０の値は、特許文献１の段落０１５４より０．０７０８と定めた。βの値は、特許文献１の段落０１５３より０．０６１５と定めた。図１０に示されるように、（Ｚ^＊ _１Ｄ）_ｉ,１が大きくなるほど、ＲＢＥも大きくなり、ピーク付近
において、ＲＢＥは概ね３となる。これらのことは単一ビームのＲＢＥに関する知見と当てはまっていることから、導出された式（９）の妥当性が見込まれる。

目的関数は式（１）に限定されない。生体における複数の位置のそれぞれに付与される線量と、複数の位置のそれぞれの目標線量との差異を定量化した値を定めたものであれば、他のものでもかまわない。例えば、他の目的関数として、以下に示す、特許文献１に記載されている式（１５）を用いてもよい。
関数[ｔ]_＋はｔ＞０のときｔをとり，それ以外では０となる関数である。
Ｄ＿ｉ（ｗ）が位置ｉの臨床線量（ＲＢＥをかけた線量）を表す。
Ｄ＿Ｐ＾ｍａｘ、Ｄ＿Ｐ＾ｍｉｎ、Ｄ＿Ｏ＾ｍａｘは、目標線量分布を表す。これらの値は、ユーザー（医療者等）が設定する。
Ｑ＿Ｐ＾ｏ、Ｑ＿Ｐ＾ｕ、Ｑ＿Ｏは適当な重み定数である。これらの値は、ユーザー（医療者等）が設定する。さらなる詳細は特許文献１に記載されているため省略する。

図１１は、本実施形態に係る情報処理装置の第１の変形例を示す。探索モード決定部２１と、線種データ記憶部２２と、探索モード切替えルール記憶部２３が追加されている。線種データ記憶部２２は、使用する線種のデータを記憶している。ユーザーは線種データ記憶部２２に、使用する線種を登録しておく。探索モード切替えルール記憶部２３は、線種ごとに、Ｌ−ＢＦＧＳ−Ｂ法を用いるか、本実施形態に係る改変Ｌ−ＢＦＧＳ−Ｂ法のいずれを用いるかを定めた探索モード切替えルールを記憶している。具体的に、探索モード切替えルールは、アクティブ集合生成ステップにおいて、ニュートン方向または最急降下方向のいずれが用いられるかを表している。

探索モード切替えルールは、あらかじめ、たとえばファントム（標的模型）を用いて行われた治療計画の試計算の結果などから定めておけばよい。アクティブ集合生成ステップでニュートン方向または最急降下方向のいずれを用いた方が良いかは、目的関数の関数形に依存することが想定される。治療計画の目的関数は生物効果モデルに依存しており、さらに生物効果モデルのモデルパラメーターは線種ごとに決まることが想定される。

探索モード切替えルールは、線種ごとに、常にニュートン方向または最急降下方向の一方のみを用いるものと定めてもよい。ニュートン方向の場合、ニュートン方向モード、最急降下方向の場合は、最急降下方向モードを決定する。または、解の更新の各回で取得される目的関数値または更新回数などから、機械学習によってニュートン方向または最急降下方向の一方を選ぶことを定めてもよい。前述した実施形態で示した重粒子線の場合は、ニュートン方向が好ましいため、線種が重粒子線の場合は、ニュートン方向モードを選択すればよい。探索モード決定部２１は、線種データ記憶部２２から線種のデータを読み出し、探索モード切替えルール記憶部から探索モード切替えルールを読み出し、探索モードを決定する。

探索方向算出部１１は、探索モード決定部２１で決定された探索モードに従って、アクティブ集合生成ステップにおいて、ニュートン方向または最急降下方向のいずれが用いられるかを決定する。それ以外の処理は、前述した実施形態と同じであるため、説明を省略する。

以上により、線種ごとに探索モードを決定することで、アクティブ集合生成ステップでニュートン方向または最急降下方向のうち適切な方向を選択できる。よって、少ない更新回数で、目的関数の良好な解を算出することができる。

図１２は、本実施形態に係る情報処理装置の第２の変形例を示す。第１の変形例に対して、パラメーター記憶部が除去され、パラメーター算出部３１と、照射設定データ記憶部３２が追加されている。照射設定データ記憶部３２は、照射設定データを記憶する。照射設定データは、放射線の照射方法を表すものであり、たとえば重粒子線による３次元スポットスキャニング法では照射方向や加速器パラメーター、スポットサイズ、線量分布の計算メッシュサイズなどを含む。パラメーター算出部３１は、照射設定データと、線種データ記憶部２２に記憶された線種のデータから、パラメーターを算出する。パラメーターとして具体的に、線量付与行列、および、ＲＢＥ関数ｒ（ｘ）の定数等を算出する。パラメーターの算出にはペンシルビーム法やモンテカルロ法などの公知の算出方法を用いることができる。

（ハードウェア構成図）
図１３は、本発明の実施形態に係る情報処理装置のハードウェア構成例を示す。図１３のハードウェア構成は、ＣＰＵ６１、入力回路６２、出力回路６３、通信回路６４、主記憶装置６５、外部記憶装置６６を備え、これらはバス６７により相互に通信可能に接続される。

外部記憶装置６６は、例えば、ＨＤＤ、ＳＳＤ、メモリ装置、ＣＤ−Ｒ、ＣＤ−ＲＷ、ＤＶＤ−ＲＡＭ、ＤＶＤ−Ｒ等の記憶媒体等を含む。外部記憶装置６６は、目的関数生成部３、探索処理部４、探索モード決定部２１、パラメーター算出部３１の機能を、プロセッサであるＣＰＵ６１に実行させるためのプログラムを記憶している。また、線量処方データ記憶部１、照射量データ記憶部５、探索モード切替えルール記憶部２３、線種データ記憶部２２、照射設定データ記憶部３２は、外部記憶装置６６または主記憶装置６５によって実現されている。ここでは、外部記憶装置６６を１つのみ示しているが、複数存在しても構わない。

主記憶装置６５は、ＣＰＵ６１による制御の下で、外部記憶装置６６に記憶された制御プログラムまたは各記憶部のデータを展開し、当該プログラムの実行時に必要なデータ、当該プログラムの実行により生じたデータ等を記憶する。主記憶装置６５は、例えば揮発性メモリ（ＤＲＡＭ、ＳＲＡＭ等）または不揮発性メモリ（ＮＡＮＤフラッシュメモリ、ＭＲＡＭ等）など、任意のメモリを含む。主記憶装置６５に展開された制御プログラムがＣＰＵ６１により実行されることで、目的関数生成部３、探索処理部４、探索モード決定部２１、パラメーター算出部３１の機能が実行される。主記憶装置６５および外部記憶装置６６は、ハードウェアストレージの一例である。

入力部である入力回路６２は、ユーザインタフェースからデータを取得する。出力部である出力回路６３は、線量処方データ記憶部１、照射量データ記憶部５、探索モード切替えルール記憶部２３、線種データ記憶部２２、照射設定データ記憶部３２内のデータを、表示装置に出力する。通信部である通信回路６４は、無線または有線の通信手段を含み、外部のサーバ等の通信装置と所定の通信方式で通信を行う。出力回路６３または通信回路６４は、本実施形態により算出した照射量を照射システムに出力または送信してもよい。通信回路６４は、主記憶装置６５または外部記憶装置６６内のデータを、通信装置に送信してもよい。入力回路６２、出力回路６３および通信回路６４は、それぞれ別個の集積回路等の回路で構成されていてもよいし、単一の集積回路等の回路で構成されてもよい。

なお、本発明は上記各実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記各実施形態に開示されている複数の構成要素を適宜組み合わせることによって種々の発明を形成できる。また例えば、各実施形態に示される全構成要素からいくつかの構成要素を削除した構成も考えられる。さらに、異なる実施形態に記載した構成要素を適宜組み合わせてもよい。

１ 線量処方データ記憶部
２ パラメーター記憶部
３ 目的関数生成部
４ 探索処理部
５ 照射量データ記憶部
１１ 探索方向算出部
１２ 改善解算出部
２１ 探索モード決定部
２２ 線種データ記憶部
２３ 探索モード切替えルール記憶部
３１ パラメーター算出部
３２ 照射設定データ記憶部
６１ ＣＰＵ
６２ 入力回路
６３ 出力回路
６４ 通信回路
６５ 主記憶装置
６６ 外部記憶装置
６７ バス

Claims (9)

1. 生体に照射する複数の放射線の照射量を表す複数の変数を含む目的関数を、前記複数の変数のそれぞれの制約条件を満たすように最適化することにより、前記複数の変数の値を決定する探索処理部を備え、
   前記目的関数は、前記生体における複数の位置のそれぞれに付与される線量と、前記複数の位置のそれぞれの目標線量との差異を定量化した値を定めており、
   前記複数の位置にそれぞれ付与される線量は、
   前記位置ごとに定義された複数の係数と、前記位置ごとに定義された、前記複数の変数を含む生物学的効果比とを、前記複数の変数に乗じて定義され、
   前記探索処理部は、現在の暫定解において前記目的関数のニュートン方向へ前記複数の変数の値を変更していくことにより、前記制約条件に違反する変数を含むアクティブ集合を生成し、前記アクティブ集合に基づいて、解の探索方向を決定し、決定した探索方向に向けて次の改善解を算出する
   情報処理装置。
2. 前記探索処理部は、前記アクティブ集合に含まれる変数を前記変数に対応する前記制約条件に違反しない値に固定して、前記現在の暫定解における前記目的関数の二次近似関数の最小解を算出し、前記最小解の方向へ探索することで、前記次の改善解を算出する
   請求項１に記載の情報処理装置。
3. 前記位置ごとに定義された前記複数の係数は、前記生体における前記複数の放射線の照射量のそれぞれに対する線量分布を表す
   請求項１または２に記載の情報処理装置。
4. 前記放射線は、粒子線、またはエックス線以外の光線である
   請求項１ないし３のいずれか一項に記載の情報処理装置。
5. 前記粒子線の前記生体への照射は、３次元スポットスキャニング法により行われる
   請求項４に記載の情報処理装置。
6. 前記粒子線は、重粒子線である
   請求項５に記載の情報処理装置。
7. 前記重粒子線は、炭素線である
   請求項６に記載の情報処理装置。
8. 生体に照射する複数の放射線の照射量を表す複数の変数を含む目的関数を、前記複数の変数のそれぞれの制約条件を満たすように最適化することにより、前記複数の変数の値を決定する探索処理ステップを備え、
   前記目的関数は、前記生体における複数の位置のそれぞれに付与される線量と、前記複数の位置のそれぞれの目標線量との差異を定量化した値を定めており、
   前記複数の位置にそれぞれ付与される線量は、
   前記位置ごとに定義された複数の係数と、前記位置ごとに定義された、前記複数の変数を含む生物学的効果比とを、前記複数の変数に乗じて定義され、
   前記探索処理ステップは、現在の暫定解において前記目的関数のニュートン方向へ前記複数の変数の値を変更していくことにより、前記制約条件に違反する変数を含むアクティブ集合を生成し、前記アクティブ集合に基づいて、解の探索方向を決定し、決定した探索方向に向けて次の改善解を算出する
   情報処理方法。
9. 生体に照射する複数の放射線の照射量を表す複数の変数を含む目的関数を、前記複数の変数のそれぞれの制約条件を満たすように最適化することにより、前記複数の変数の値を決定する探索処理ステップをコンピューターに実行させ、
   前記目的関数は、前記生体における複数の位置のそれぞれに付与される線量と、前記複数の位置のそれぞれの目標線量との差異を定量化した値を定めており、
   前記複数の位置にそれぞれ付与される線量は、
   前記位置ごとに定義された複数の係数と、前記位置ごとに定義された、前記複数の変数を含む生物学的効果比とを、前記複数の変数に乗じて定義され、
   前記探索処理ステップは、現在の暫定解において前記目的関数のニュートン方向へ前記複数の変数の値を変更していくことにより、前記制約条件に違反する変数を含むアクティブ集合を生成し、前記アクティブ集合に基づいて、解の探索方向を決定し、決定した探索方向に向けて次の改善解を算出する
   コンピュータプログラム。