(57)【要約】【課題】画像に含まれる細胞数が不明な状態から初期条件を個別に設定することなく高速に十分な精度で安定して細胞数およびそれぞれの細胞の形状、活動時系列を推定する技術を提供する。【解決手段】細胞の活動を非負の値で表現される細胞形状ac(x,y)と非負の値で表現されるスパイク時系列uc(t)より導出されるカルシウム濃度変化vc(t)の積により表現し、観測データf(t,x,y)を複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データとモデルと事前分布を元にしたベイズ理論に基づく推定により想定される細胞の位置と活動を求める情報処理装置。【選択図】図１
2014095678.tif 000002
【特許請求の範囲】
【請求項１】
細胞の活動を非負の値で表現される細胞形状ac(x,y)と非負の値で表現されるスパイク時系列uc(t)より導出されるカルシウム濃度変化vc(t)の積により表現し、観測データf(t,x,y)を複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データとモデルと事前分布を元にしたベイズ理論に基づく推定により想定される細胞の位置と活動を求める情報処理装置であって、
細胞のスパイク生成に起因する細胞内のカルシウム濃度の上昇に伴って起こる蛍光タンパク質からの発光を位置(x,y)における時刻tでの蛍光強度f(t,x,y)として読み取るデータ読み取り部と、
スパイク時系列uc(t)とカルシウム濃度変化vc(t)との関係式と細胞形状ac(x,y)とスパイク時系列uc(t)に関する事前分布を設定するパラメータ設定部と、
前記データ読み取り部により読みとられた位置(x,y)における時刻tでの蛍光強度f(t,x,y)に基づいて細胞形状ac(x,y)の細胞候補を配置し、ベースラインを蛍光強度f(t,x,y)の統計量に基づいて設定する初期設定部と、
前記ベイズ理論に基づく推定で導出される最小化問題の目的関数を減少させ、細胞形状ac(x,y)、カルシウム濃度変化vc(t)およびベースラインをより真の解に近い推定値に更新する演算部と、
前記演算部における更新で目的関数の変化が十分小さくなった場合に目的関数が収束したと判定して処理を終了する収束判定部とを有することを特徴とする情報処理装置。
【請求項２】
前記演算部において、
空間ベースラインabaseline(x,y)に加えて時間ベースラインvbaseline(t)を導入したモデルに基づいた推定値更新を行うことを特徴とする請求項１に記載の情報処理装置。
【請求項３】
前記パラメータ設定部において、
スパイク時系列に対するカルシウム濃度変化のインパルス応答に基づいた変換行列により、スパイク時系列とカルシウム濃度変化との関係式を設定し、
前記演算部において、
変換行列に基づいた推定値更新を行うことを特徴とする請求項１又は２に記載の情報処理装置。
【請求項４】
前記パラメータ設定部において、
2014095678.tif 000037
前記演算部において、
設定された事前分布に基づいた推定値更新を行うことを特徴とする請求項１から３までのいずれか１項に記載の情報処理装置。
【請求項５】
前記初期設定部において、
想定される細胞のサイズより大きめの特定の領域を覆うような細胞候補を十分な重なりを持たせて画像全体を埋め尽くすように真の細胞の位置に無関係に多数配置することを特徴とする請求項１から４までのいずれか１項に記載の情報処理装置。
【請求項６】
さらに、前記スパイク時系列／時間方向のベースライン演算部と、前記細胞形状／空間方向のベースライン演算部と、の出力の少なくともいずれか一方に対して縮退した細胞候補の除去を行う縮退細胞候補除去部を有することを特徴とする請求項１から５までのいずれか１項に記載の情報処理装置。
【請求項７】
さらに、前記スパイク時系列／時間方向のベースライン演算部と、前記細胞形状／空間方向のベースライン演算部と、の出力の少なくともいずれか一方に対して類似した細胞候補を統合する類似細胞統合部
を有することを特徴とする請求項１から６までのいずれか１項に記載の情報処理装置。
【請求項８】
前記演算部は、
細胞形状ac(x,y)と空間ベースラインabaseline(x,y)を固定し、スパイク時系列の変数値が非負という拘束条件のもとで、MAP推定の目的関数を最小化するカルシウム濃度変化vc(t)と時間ベースラインvbaseline(t)を求めるCT+T次元の二次計画問題を解くカルシウム濃度変化／時間方向ベースライン演算部と、
カルシウム濃度変化vc(t)と時間ベースラインvbaseline(t)を固定し、細胞形状の変数値が非負という拘束条件のもとで、MAP推定の目的関数を最小化する細胞形状ac(x,y)と空間ベースラインabaseline(x,y)を求めるXY個の(C+1)次元の二次計画問題を解く細胞形状／空間方向ベースライン演算部と
を有することを特徴とする請求項１から７までのいずれか１項に記載の情報処理装置。
【請求項９】
前記カルシウム濃度変化／時間方向ベースライン演算部において、
当該二次計画問題を主双対内点法と共役勾配法を用いて反復的に解くことを特徴とする請求項８に記載の情報処理装置。
【請求項１０】
細胞の活動を非負の値で表現される細胞形状ac(x,y)と非負の値で表現されるスパイク時系列uc(t)より導出されるカルシウム濃度変化vc(t)の積により表現し、観測データf(t,x,y)を複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データとモデルと事前分布を元にしたベイズ理論に基づく推定により想定される細胞の位置と活動を求める情報処理装置を用いた情報処理方法であって、
細胞のスパイク生成に起因する細胞内のカルシウム濃度の上昇に伴って起こる蛍光タンパク質からの発光位置(x,y)における時刻tでの蛍光強度f(t,x,y)を読み取るデータ読み取りステップと、
スパイク時系列uc(t)とカルシウム濃度変化vc(t)との関係式と細胞形状ac(x,y)とスパイク時系列uc(t)に関する事前分布を設定するパラメータ設定ステップと、
前記データ読み取り部により読みとられた位置(x,y)における時刻tでの蛍光強度f(t,x,y)に基づいて細胞形状ac(x,y)の細胞候補を配置し、ベースラインを蛍光強度f(t,x,y)の統計量に基づいて設定する初期設定ステップと、
前記ベイズ理論に基づく推定で導出される最小化問題の目的関数を減少させ、細胞形状ac(x,y)、カルシウム濃度変化vc(t)およびベースラインをより真の解に近い推定値に更新する演算ステップと、
前記演算ステップにおける更新で目的関数の変化が十分小さくなった場合に目的関数が収束したと判定して処理を終了する収束判定ステップと
を有することを特徴とする情報処理方法。
【請求項１１】
コンピュータに請求項１０に記載の情報処理方法を実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、情報処理技術に関し、特に、神経細胞の活動電位を画像処理等により推測する情報処理技術に関する。
【背景技術】
【０００２】
神経細胞は活動電位（スパイク）を発生した際に他の神経細胞に情報を伝達する。脳内の神経細胞集団がどのように情報処理を行っているかを理解するためには複数の神経細胞の活動電位を同時に記録する必要がある。活動電位を発生する際には細胞内のカルシウム濃度が上昇するため、カルシウム濃度に依存して光る蛍光タンパク質を導入した脳を二光子顕微鏡で動画として記録することにより、複数の細胞に活動電位が生じている様子を観察することができる。
【０００３】
しかしながら、実際に神経細胞集団のデータとして扱うためには動画データから、各細胞の位置と活動の時系列を抽出する必要がある。実際には、データにはノイズが多く含まれるため、得られるのは確定した情報ではなく、データを元に推定された情報である。大容量の記録データから、高速・高精度に細胞位置と活動時系列を推定することが重要である。
【０００４】
従来手法として主に用いられている手法は、まず主成分分析などを用いて類似した活動パターンを見せる領域を見つけ、いくつかの処理をすることでROI(Region of Interest)と呼ばれる個々の細胞の形状を確定させ、その領域内での蛍光タンパク質の反応の平均をその細胞の活動と見なし、必要であればスパイクタイミングを正確に求めるという手順を取るものである（下記非特許文献１参照）。当該非特許文献は、主成分分析/独立成分分析によるROI決定アルゴリズムに関する。
【０００５】
実際には細胞集団は３次元の構造を持ち２次元で射影した画像データ中では複数の細胞が重なっている領域が多数あるが、そのクロストークを考慮していないことを解決するために、本発明のベースとなるモデルに関する論文として、近年、観測画像のモデルとして複数細胞の活動の総和を用いる考え方が提唱され、そのモデルに基づいて初期条件を与えて反復法により収束させることで細胞の位置と活動を推定できることをデモンストレーションにより示させている（下記非特許文献２参照）。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】Automated anlalysis of cellular signals from large-scale calcium imaging data, E. A. Mukamel, A. Nimmerjahn and M. J. Schnitzer, Neuron 63, 747-760, (2009).
【非特許文献２】Fast nonnegative deconvolution for spike train inference from population calcium imaging, J. T. Vogelstein, A. M. Packer, T. A. Machado, T. Sippy, B. Babadi, R Yuste and L Paninski, J. Neurophysiol 104:3691-3704 (2010).
【発明の概要】
【発明が解決しようとする課題】
【０００７】
非特許文献１のような方法の問題点は、実際には細胞集団は３次元の構造を持ち２次元で射影した画像データ中では複数の細胞が重なっている領域が多数あるが、そのクロストークを考慮していないことである。
【０００８】
また、非特許文献２の方法は、現実の大容量データを処理するには多くの問題点が存在する。
【０００９】
以下に、非特許文献２に開示されている既存処理の内容とその問題点について図面を参照しながら説明を行う。
（既存処理）
図４は、非特許文献２に示される既存処理を行う情報処理装置の一構成例を示す機能ブロック図であり、図５は、非特許文献２に示される既存処理の流れを示すフローチャート図である。まず、ステップＳ１０１で情報処理を開始すると、ステップＳ１０２で、データ読み込み部１０１が、データの読み込みを行う。読み込むデータとして、位置(x,y)における時刻tでの蛍光強度f(t,x,y)がデータとして与えられる。ただし、t,x,yはそれぞれT,X,Y以下の自然数で、サンプリングレートおよび空間解像度から決まるスケールになっているとする。蛍光強度fは相対的な値であるため、何らかのベースラインを規定する必要がある。
【００１０】
2014095678.tif 000003
【００１１】
【数１】
2014095678.tif 000004
【００１２】
時間軸に対するカルシウム濃度はスパイク（図３参照）によりデルタ関数的に増加し、その後、減衰する。uc(0)=vc(0)よりuc,vcは一方の時系列からもう一方が線形変換で導出可能である。
【００１３】
【数２】
2014095678.tif 000005
【００１４】
2014095678.tif 000006
【００１５】
【数３】
2014095678.tif 000007
【００１６】
2014095678.tif 000008
【００１７】
【数４】
2014095678.tif 000009
【００１８】
2014095678.tif 000010
ベースラインは平均あるいはmedianで初期設定する。
【００１９】
【数５】
2014095678.tif 000011
vc(t)=0として、式（４）よりσ2を設定する。
2014095678.tif 000012
【００２０】
但し、細胞の個数Cおよび細胞の形状ac(x,y)の求め方は規定されておらず、これらの初期条件を得る何らかの手順を別途設定する必要があるが、解が適切に収束するためには初期条件が解に十分近い必要がある。
【００２１】
次いで、ステップＳ１０４において、カルシウム濃度／スパイク時系列演算部１０５が、カルシウム濃度vおよびスパイク時系列uを求める。上記の式(1),(2),(3)において、ac,abaselineおよびパラメータσ,λcを固定して、u,vに関するMAP(事後分布最大化)推定を行うと、式（２）の条件下で次の二次計画問題を解くことになる。このとき、ベースラインが固定されており、かつ非負拘束条件が設定されているため、固定されたベースラインの影響を強く受け系全体としての適切な解が得られない。
【００２２】
【数６】
2014095678.tif 000013
【００２３】
Log障壁法とニュートン法を用いてこの最小化問題を解く。Ｃ＝１ならば３重対角行列になり比較的容易に解けるが、Ｃ＞１では３重ブロック対角行列を考える必要があり、精度よく解くことは困難である。
【００２４】
次いで、ステップＳ１０５において、正規化部１０７が、各細胞候補のスパイク時系列と細胞形状を正規化する処理を行う。
【００２５】
各細胞の信号はfc(t,x,y)=ac(x,y)vc(t)と表されるが、
【数７】
2014095678.tif 000014
としても信号全体は変わらない。
【００２６】
しかし、式（７）は最小二乗の項とuに対するペナルティ項の和になっているため、式（２）より、scを大きくすればするほどペナルティ項が小さくなり目的関数の値は小さくなってしまう。これではペナルティ項の意味がないため、何らかの形で正規化する必要がある。そこで、sc=maxvcとする。しかし、拘束条件sc=maxvcをいれて問題を解いたわけではないため、正規化後は式（７）の解になっていない。
【００２７】
ステップＳ１０６で、細胞形状／ベースライン演算部が、細胞形状acおよびベースラインabaselineを求める。ここではu,vを固定して、式（７）が最小となるac,abaselineを求める。最尤法を行うと、ピクセルごとに独立して考えることができて、（９）式をそれぞれ解けばよい。
【００２８】
【数８】
2014095678.tif 000015
【００２９】
これは、拘束条件がないので簡単に解けるが、ac(x,y)には負の値も含まれるため細胞の形状はノイズの影響を強く受け、安定しない。
【００３０】
ステップＳ１０７で、パラメータ更新部１１１が、パラメータの更新を行う。式（４）、（５）に従いパラメータを更新する。しかしながら、真の解で規定されるべきパラメータを現在の仮の値で近似しており、事前分布が変数から定義されることになるため、目的関数を正しく設定できない。
【００３１】
ステップＳ１０４とＳ１０７の処理を、反復・収束判定部１１７が、反復・収束判定する。ステップＳ１０８で処理が終了すると、出力部１２１が結果を出力する。
【００３２】
正規化処理の説明においても述べたように、そもそも問題設定が不良のため最終的な収束判定が不能であるため、各ステップを適当な回数繰り返したところで打ち切ることになる。
【００３３】
このように、非特許文献２の問題点としては、アルゴリズム中で細胞数が固定であること、初期条件として細胞の形状を正しい形状に十分近くしなければならないこと、示された解法では速度および精度が不十分であること、多数の細胞が存在する際に適切な収束が得られないこと、などである。
【００３４】
本発明は、画像に含まれる細胞数が不明な状態から初期条件を個別に設定することなく高速に十分な精度で安定して細胞数およびそれぞれの細胞の形状、活動時系列を推定する技術を提供することを目的とする。
【課題を解決するための手段】
【００３５】
本発明では、パラメータの非負拘束条件と適切な事前分布の設定による正規化項（ペナルティ項）をモデルに導入し得られる解の精度を向上させた上で、多数の細胞候補を画像中に重なりを持たせて密に埋め尽くす様な初期条件設定と、細胞候補の取捨選択の工夫により、細胞数とパラメータを正しく推定できる装置（プログラム）を開発した。主たるアルゴリズムとしては、形状と時間変化を交互に反復して改善する反復法を用い、途中でノイズレベル以下の小さな信号しか持たない細胞候補を取り除き、類似した活動状態を持つ細胞候補を統合することで計算量を削減しながら画像中に含まれる活動した細胞をすべて正しく検出する。各ステップで必要な非線形計画の解法には、主双対内点法および大規模一次方程式の高速な解法としての共役勾配法を用いた独自の並列化実装を用いる。
【００３６】
最終的に、画像に含まれる細胞数が不明な状態から初期条件を個別に設定することなく高速に十分な精度で安定して細胞数およびそれぞれの細胞の形状、活動時系列を推定することができるプログラムが得られた。
【００３７】
従来法２では、細胞数が固定かつ初期条件が厳しいという問題点がある。そこで、多数の細胞候補を画像中に重なりを持たせて密に埋め尽くす様に初期条件を設定し、解が疎になるような条件をモデルに付け加える。前者の工夫により、真の細胞はいずれかの細胞候補の十分近くにあることになるため取りこぼすことなく細胞を検出でき、かつ後者の工夫とクロストークをきちんと考慮したモデルを用いている効果により多数の細胞候補があっても必要以上に真の細胞が複数の細胞候補に分割されてしまうことを防ぐことができる。実際、アルゴリズムの反復を繰り返すと関連度自動決定と呼ばれる効果により大きな信号を持つ細胞候補と、ノイズレベル以下の小さな信号しか持たない細胞候補が生じるようになるため、小さな信号の細胞候補を取り除き、類似した活動状態を持つ細胞候補を統合することで最終的に画像中に含まれる活動した細胞をすべて正しく検出できる。ここで新しく提案した手法は、観測画像が各細胞の活動の和であるというモデルに対して適切な拘束条件を入れることで、反復法が正しく収束することをサポートしている。
【００３８】
この手法では、多数の細胞候補を用意するという性質上、多くのメモリおよび計算量を必要とすることになるが、画像を適切なサイズに分割して前処理しそれらの結果を統合する方針をとることで計算量を大幅に削減することができ、個々のステップでの適切な解法と並列化による高速化した実装により現実的な実行時間で目的を達成することができるようになっている。
【００３９】
最終的に、画像に含まれる細胞数が不明な状態から初期条件を個別に設定することなく高速に十分な精度で安定して細胞数およびそれぞれの細胞の形状、活動時系列を推定することができる様になった。
【００４０】
本発明の一観点によれば、細胞の活動を非負の値で表現される細胞形状ac(x,y)と非負の値で表現されるスパイク時系列uc(t)より導出されるカルシウム濃度変化vc(t)の積により表現し、観測データf(t,x,y)を複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データとモデルと事前分布を元にしたベイズ理論に基づく推定により想定される細胞の位置と活動を求める情報処理装置であって、細胞のスパイク生成に起因する細胞内のカルシウム濃度の上昇に伴って起こる蛍光タンパク質からの発光を位置(x,y)における時刻tでの蛍光強度f(t,x,y)として読み取るデータ読み取り部と、スパイク時系列uc(t)とカルシウム濃度変化vc(t)との関係式と細胞形状ac(x,y)とスパイク時系列uc(t)に関する事前分布を設定するパラメータ設定部と、前記データ読み取り部により読みとられた位置(x,y)における時刻tでの蛍光強度f(t,x,y)に基づいて細胞形状ac(x,y)の細胞候補を配置し、ベースラインを蛍光強度f(t,x,y)の統計量に基づいて設定する初期設定部と、前記ベイズ理論に基づく推定で導出される最小化問題の目的関数を減少させ、細胞形状ac(x,y)、カルシウム濃度変化vc(t)およびベースラインをより真の解に近い推定値に更新する演算部と、前記演算部における更新で目的関数の変化が十分小さくなった場合に目的関数が収束したと判定して処理を終了する収束判定部と、を有することを特徴とする情報処理装置が提供される。
【００４１】
前記演算部において、空間ベースラインabaseline(x,y)に加えて時間ベースラインvbaseline(t)を導入したモデルに基づいた推定値更新を行うことを特徴とする。
【００４２】
また、前記パラメータ設定部において、スパイク時系列に対するカルシウム濃度変化のインパルス応答に基づいた変換行列により、スパイク時系列とカルシウム濃度変化との関係式を設定し、前記演算部において、変換行列に基づいた推定値更新を行うことを特徴とする。
【００４３】
また、前記パラメータ設定部において、細胞cの細胞形状ac(x,y)とスパイク時系列uc(t)の事前分布を細胞の平均サイズ、スパイクの発火率、それらの積である信号強度を考慮した指数分布として設定し、前記演算部において、設定された事前分布に基づいた推定値更新を行うことを特徴とする。
【００４４】
また、前記演算部においてスパイク時系列uc(t)に加えて細胞形状ac(x,y)に対する非負拘束条件を導入したモデルに基づいた推定値更新を行うことを特徴とする。
【００４５】
これにより、疎な解を得られるような条件で推定を行うことで関連度自動決定の効果が働き、不要な細胞候補は自動的に縮退し（信号が小さくなり）細胞数が自動的に決定される。
【００４６】
また、前記初期設定部において、想定される細胞のサイズより大きめの特定の領域を覆うような細胞候補を十分な重なりを持たせて画像全体を埋め尽くすように真の細胞の位置に無関係に多数配置することを特徴とする。
【００４７】
さらに、前記スパイク時系列／時間方向のベースライン演算部と、前記細胞形状／空間方向のベースライン演算部と、の出力の少なくともいずれか一方に対して縮退した細胞候補の除去を行う縮退細胞候補除去部を有することを特徴とする。
【００４８】
さらに、前記スパイク時系列／時間方向のベースライン演算部と、前記細胞形状／空間方向のベースライン演算部と、の出力の少なくともいずれか一方に対して類似した細胞候補を統合する類似細胞統合部を有することを特徴とする。
【００４９】
また、前記演算部は、定められたモデル、拘束条件、事前分布のもとでMAP推定の目的関数を減少させる演算を行う。尚、MAP推定はベイズ推定という枠組みの中の比較的簡易な推定方法の例示であり、広く「ベイズ理論に基づく推定」を用いることができる。具体的には、定められたモデル、拘束条件、事前分布のもとで事後分布を近似したテスト関数が真の事後分布に近づくような演算を行う。
【００５０】
このベイズ推定法による定量的カルシウム動態推定については、以下発表を参照することができる。
角田敬正・織田善晃・大森敏明・井上雅司・宮川博義・岡田真人・青西亨、信学技報, vol. 111, no. 483, NC2011-175, pp. 317-322, 2012年3月.
http://www.ieice.org/ken/paper/20120316a0pt/
【００５１】
さらに、前記演算部は、細胞形状ac(x,y)と空間ベースラインabaseline(x,y)を固定し、スパイク時系列の変数値が非負という拘束条件のもとで、MAP推定の目的関数を最小化するカルシウム濃度変化vc(t)と時間ベースラインvbaseline(t)を求める(C+1)T次元の二次計画問題を解くカルシウム濃度変化／時間方向ベースライン演算部と、カルシウム濃度変化vc(t)と時間ベースラインvbaseline(t)を固定し、細胞形状の変数値が非負という拘束条件のもとで、MAP推定の目的関数を最小化する細胞形状ac(x,y)と空間ベースラインabaseline(x,y)を求めるXY個の(C+1)次元の二次計画問題を解く細胞形状／空間方向ベースライン演算部とを有することを特徴とする。
【００５２】
前記カルシウム濃度変化／時間方向ベースライン演算部において、当該二次計画問題を主双対内点法と（前処理付き）共役勾配法を用いて反復的に解くことを特徴とする。
【００５３】
本発明の他の観点によれば、細胞の活動を非負の値で表現される細胞形状ac(x,y)と非負の値で表現されるスパイク時系列uc(t)より導出されるカルシウム濃度変化vc(t)の積により表現し、観測データf(t,x,y)を複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データとモデルと事前分布を元にしたベイズ理論に基づく推定により想定される細胞の位置と活動を求める情報処理装置を用いた情報処理方法であって、細胞のスパイク生成に起因する細胞内のカルシウム濃度の上昇に伴って起こる蛍光タンパク質からの発光位置(x,y)における時刻tでの蛍光強度f(t,x,y)を読み取るデータ読み取りステップと、スパイク時系列uc(t)とカルシウム濃度変化vc(t)との関係式と細胞形状ac(x,y)とスパイク時系列uc(t)に関する事前分布を設定するパラメータ設定ステップと、前記データ読み取り部により読みとられた位置(x,y)における時刻tでの蛍光強度f(t,x,y)に基づいて細胞形状ac(x,y)の細胞候補を配置し、ベースラインを蛍光強度f(t,x,y)の統計量に基づいて設定する初期設定ステップと、前記ベイズ理論に基づく推定で導出される最小化問題の目的関数を減少させ、細胞形状ac(x,y)、カルシウム濃度変化vc(t)およびベースラインをより真の解に近い推定値に更新する演算ステップと、前記演算ステップにおける更新で目的関数の変化が十分小さくなった場合に目的関数が収束したと判定して処理を終了する収束判定ステップとを有することを特徴とする情報処理方法が提供される。
【００５４】
本発明は、コンピュータに上記に記載の情報処理方法を実行させるためのプログラムであっても良く、当該プログラムを記録するコンピュータ読み取り可能な記憶媒体であっても良い。
【発明の効果】
【００５５】
本発明によれば、画像に含まれる細胞数が不明な状態から初期条件を個別に設定することなく高速に十分な精度で安定して細胞数およびそれぞれの細胞の形状、活動時系列を推定することができる。
【図面の簡単な説明】
【００５６】
【図１】本発明の一実施の形態による情報処理装置の一構成例を示す機能ブロック図である。
【図２Ａ】本発明の一実施の形態による情報処理装置による情報処理の流れを示すフローチャート図である。
【図２Ｂ】本発明の一実施の形態による情報処理装置による情報処理の流れを示すフローチャート図であり、図２Ａに続く図である。
【図３】本発明の一実施の形態による情報処理装置による情報処理の概要を示す模式的な図である。
【図４】非特許文献２に記載の情報処理装置の一構成例を示す機能ブロック図である。
【図５】非特許文献２に記載の情報処理装置による情報処理の流れを示すフローチャート図である。
【発明を実施するための形態】
【００５７】
以下、本発明の一実施の形態による情報処理技術について図面を参照しながら説明を行う。
本発明は、パラメータの非負拘束条件と適切な事前分布の設定による正規化項（ペナルティ項）をモデルに導入し得られる解の精度を向上させた上で、多数の細胞候補を画像中に重なりを持たせて密に埋め尽くす様な初期条件設定と、細胞候補の取捨選択の工夫により、細胞数とパラメータを正しく推定できる技術に関するものである。主たるアルゴリズムとしては、形状と時間変化を交互に反復して改善する反復法を用い、途中でノイズレベル以下の小さな信号しか持たない細胞候補を取り除き、類似した活動状態を持つ細胞候補を統合することで計算量を削減しながら画像中に含まれる活動した細胞をすべて正しく検出する。各ステップで必要な非線形計画の解法には、主双対内点法および大規模一次方程式の高速な解法としての共役勾配法を用いた独自の並列化実装を用いる。
【００５８】
図３は、本実施の形態による情報処理の概要を示す模式的な図である。多数の画像からなる動画像を観察し、ある時間ｔとその時間ｔにおける活動電位を発生する際における細胞内のカルシウム濃度の上昇に起因する蛍光タンパク質からの発光を二光子顕微鏡で観察した動画像を元に、蛍光強度f(t,x,y)を得ることができる。
【００５９】
この蛍光強度f(t,x,y)の観測値を読み込み、初期条件を設定して、細胞の形状とスパイク列の対応するカルシウム濃度との乗算値の総和を求めることで、細胞数及びそれぞれの細胞の形状、活動時系列を求めることができる。
【００６０】
本実施の形態では、既存手法の問題点を解消するため、以下のような工夫を行った。
１）既存手法のモデルに条件等を追加して問題設定を明確化し、解が一意に求まるようにする。
２）併せて、目的に合致した解の性質を持つようにモデルを改良する。
３）細胞数および細胞位置が未定である状態を想定し、多数の細胞候補から関連度自動決定により適切な細胞数および細胞位置を正しく推定できるようにする。
４）大規模なデータに対しても高精度に推定が行えるようなモデルに即した解法を示す。
【００６１】
図１は、本実施の形態による情報処理装置（画像処理装置）の一構成例を示す機能ブロック図であり、図２Ａ、図２Ｂは、情報処理の流れを示す機能ブロック図であり、図３は、情報処理の概要を示す模式的な図である。
【００６２】
まず、処理を開始すると（ステップＳ１）、ステップＳ２において、データ読み出し部１が、ある時間ｔとその時間ｔにおける活動電位を発生する際における細胞内のカルシウム濃度が上昇に起因する蛍光タンパク質からの発光を二光子顕微鏡で観察した動画像データを読み込み、パラメータを設定する。読み込むデータとしては、位置(x,y)における時刻tでの蛍光強度f(t,x,y)がデータとして与えられる。ただし、t,x,yはそれぞれT,X,Y以下の自然数で、サンプリングレートおよび空間解像度から決まるスケールになっているとする。尚、サンプリングレートＳ[Hz]、記録時間Ｒ[sec]とすると、Ｔ＝ＳＲである。また、Ｘ、Ｙは画像のピクセル数である。蛍光強度fは相対的な値であるため、何らかのベースラインを規定する必要がある。解析モデル記憶部２１に記憶される変数を更新する本実施の形態によるモデルは、既存手法のモデルに対して以下の変更を行ったモデルである。
【００６３】
１）時間方向のベースラインvbaseline(t)を導入
既存手法のモデルでは、式（１）のように空間方向のみについてベースラインを考慮したが、本実施の形態では、時系列推定時に時間ベースラインvbaseline(t)を同時に推定することで、解が正しく収束するようにする。
【００６４】
2014095678.tif 000016
【００６５】
【数９】
2014095678.tif 000017
【００６６】
２）
既存手法では、ucとvcの対応を式（2）の形に限定していたが、本実施の形態では、より一般的な形に拡張することで、現実のカルシウム応答に近づけより精度の高い解が得られるようにする。
【００６７】
拡張したモデルでは、カルシウム濃度vc=[vc(1) vc(2) ・・・ vc(T)]Tはインパルス応答 [g1,g2,・・・]とスパイク時系列uc=[uc(1) uc(2) ・・・ uc(T)]Tとのたたみ込みで記述される。
【００６８】
【数１０】
2014095678.tif 000018
【００６９】
2014095678.tif 000019
【００７０】
【数１１】
2014095678.tif 000020
【００７１】
たとえば、既存の方法ではインパルス応答を指数関数で{1,γ,γ2,・・・,γT}と定義していた。このとき、R=1である。
【００７２】
【数１２】
2014095678.tif 000021
【００７３】
現実的な拡張としては、カルシウム濃度の立ち上がりにデルタ関数ではなく有限の値の時定数γ1での立ち上がりを仮定した{1,(γ12-γ22)/(γ1-γ2),(γ13-γ23)/(γ1-γ2),・・・,(γ1T-γ2T)/(γ1-γ2)}という形のインパルス応答を定義して以下の式(14)を用いることが考えられる。このとき、R=2 である。
【００７４】
【数１３】
2014095678.tif 000022
【００７５】
このように、本実施の形態によれば、わずかな計算量の増加でモデルを拡張することができる。
【００７６】
３）acに対する非負拘束条件を導入
上記既存手法では、式(3)および(5)により結果的にペナルティ項が定数となってしまい、解が一意に定まらない。スパイク時系列uc(t)には非負条件が課されているが、細胞の形状ac(x,y)の拘束条件が規定されておらず負の値もとり得るため、細胞とは実際には無関係な領域の影響を受けてしまう。本実施の形態では、本来の物理的な条件を反映させ細胞形状ac(x,y)にも非負の拘束条件を導入する。対応する細胞に対する相関が十分強い領域のみが値を持ち、それ以外がすべて0となる（解が疎になる）ため、細胞の形状が明確になる（図３参照）。
【００７７】
【数１４】
2014095678.tif 000023
【００７８】
2014095678.tif 000024
既存手法では、式(3)の事前分布のハイパーパラメータが式(5)により変数自身により決定される自己一貫性が要求されており、解が一意に定まらなかった。また、スパイク時系列のみに事前分布が設定されているため、信号を形成する細胞形状とスパイク時系列の比も一意に定まらなかった。本実施の形態では、細胞形状とスパイク時系列に関して式(16)の形式の事前分布を導入する。
【００７９】
【数１５】
2014095678.tif 000025
であり、a0,u0,s0は事前分布のハイパーパラメータである。この事前分布はモデルにおいてスパイク時系列の発火率および細胞形状の領域サイズさらにそれらの積で表される細胞の信号強度に対してハイパーパラメータで表現される事前知識を用いることに対応し、MAP推定においてはac,uc,acucに関するL1ノルム正規化項を加えることに対応する。この事前分布の導入により、MAP推定における解が一意に定まり、式(12)(15)と併せてより疎な解を得ることができる。
【００８０】
目的関数とパラメータ設定処理においては、拘束条件（式(12),(15)）、尤度関数（式(10)）、事前分布（式(16)(17)）のMAP推定は、拘束条件（式(12),(15)）の元で目的関数式(18)を最小化する細胞数Cと変数ac,abaseline,vc,vbaselineを求める問題となる。
【００８１】
【数１６】
2014095678.tif 000026
【００８２】
σ2は真の解が得られたときの二乗誤差であり現時点では不明であるため、全体を定数σ2倍した上で、パラメータρa,ρu,ρを用いて次のように置き換える。
【００８３】
【数１７】
2014095678.tif 000027
【００８４】
ただし、任意の正の定数sに対して、
【数１８】
2014095678.tif 000028
としても目的関数の値も個々の細胞の信号も変わらずac,ucの比を変化させるだけであるため、ρa,ρuの比は解の実質的な形に影響しない。このため、パラメータγ,ρa,ρu,ρを定めて式(19)を最小化する解を求めておいてから、必要であれば正の定数sを用いて(20)の変換を行うことにより最終的なac,vcのスケールを決めればよい。(19)を最小化する解は、適当なac,abaselineを初期条件として設定した上で、ac,abaselineを固定することにより(19)をvc,vbaselineに関する二次計画問題として解く時間変数更新ステップと、vc,vbaselineを固定することにより式(19)をac,abaselineに関する二次計画問題として解く空間変数更新ステップとを交互に繰り返すことにより、目的関数を減少させ最終的に収束させることで求める。ベースラインの初期条件としては、最小二乗法により簡単に求めることのできる細胞が存在しない場合の最適値を設定し、細胞形状に関しては想定される細胞のサイズより大きめの特定の領域を覆うような細胞候補を十分な重なりを持たせて画像全体を埋め尽くすように多数配置する。不要な細胞候補は関連度自動決定により縮退するため、自動的に正しい細胞数が得られる。
【００８５】
パラメータは、予想される細胞の画像上での大きさAsize、発火率UfreqおよびSN比Sに基づいて設定する。具体的にはac,ucの解がそれぞれ２つの値{0,σ},{0,S}しかとらず、σの値をとるピクセル数がAsize、Sの値をとる確率がUfreqであることを想定して、ρa≒σAsize,ρu≒STUfreqと設定する。
【００８６】
上記のアルゴリズムを用いることで、ここで設定した問題の大域解を高精度に求めることができる。
【００８７】
次いで、ステップＳ３において、初期設定部３が、ベースライン、細胞の個数および細胞の形状を初期設定する。例えば、細胞が存在しないと仮定した場合の最適なベースラインを初期条件として設定する。
【００８８】
ここで、想定される細胞のサイズより大きめの特定の領域を覆うような細胞候補を十分な重なりを持たせて画像全体を埋め尽くすように真の細胞の位置に無関係に多数配置する。疎な解を得られるような条件でMAP推定を行うことで関連度自動決定の効果が働き、不要な細胞候補は自動的に縮退し（信号が小さくなり）細胞数が自動的に決定される。
【００８９】
次いで、ステップＳ４において、スパイク時系列／時間方向のベースライン演算部５が、スパイク時系列および時間方向のベースラインを求める。
【００９０】
スパイク時系列／時間方向のベースライン演算部５は、拘束条件（式(2)）のもとで、下記の目的関数（式（１７））を最小化するvc,vbaseline,ucを求めるCT+T次元の二次計画問題を解く。
【００９１】
ac,abaselineを固定すると、拘束条件(12)(15)のもとで目的関数(21)を最小化するvc,vbaselineを求めるCT+T次元の二次計画問題となる。この問題を主双対内点法と前処理付き共役勾配法を組み合わせて解く。
【００９２】
【数１９】
2014095678.tif 000029
ただし、etはt番目の要素が1となる規定ベクトルとし、パラメータは次の式で表される。
【００９３】
【数２０】
2014095678.tif 000030
【００９４】
・主双対内点法
二次計画問題(21)の基本的な解法として主双対内点法を導入する。主双対内点法では変数vc,vbaseline,ucおよびスラック変数scに対して初期条件を設定し、それらの値を目的の二次計画問題の解に近づくように更新していく。更新する方向は次の線形方程式の解により決定される。ここで、βは制御パラメータ、eは要素がすべて１のベクトルである。
【００９５】
【数２１】
2014095678.tif 000031
2014095678.tif 000032
【００９６】
【数２２】
2014095678.tif 000033
このとき、(21)からΔvbaseline,Δs,Δuを消去すると、CT次元の線形方程式(28)が得られる。
【００９７】
【数２３】
2014095678.tif 000034
ただし、
【数２４】
2014095678.tif 000035
である。
【００９８】
CT次元の線形方程式は通常の解法では膨大なメモリと計算量が必要となるが、共役勾配法では対象となる行列と任意のベクトルの積が計算できればよいので、(28)のA’およびQijのみ記憶し、その構造を利用して消費メモリO(C2+CTR)および計算量O(C2T+CTR)で行列とベクトルの積を計算することにより、（前処理付きの）共役勾配法を適用する。通常CT×CT次元の行列とCT次元のベクトルの積の計算量および消費メモリはO(C2T2)である。
【００９９】
一方、ステップＳ９において、細胞形状／空間ベースライン演算部５が、細胞形状および空間方向のベースラインを求める。
【０１００】
vc,vbaselineを固定すると、(15)の拘束条件のもとで目的関数(19)を最小化するac,abaselineを求めるXY個の(C+1)次元の二次計画問題となる。この問題を主双対内点法を用いて解く。ここでは、時間方向と違い問題の次元は比較的小さいので線形方程式の解法にはコレスキー分解を用いた直接法を用いる。このとき、各問題を並列処理で解くことで速度を向上させることができる。
【０１０１】
【数２５】
2014095678.tif 000036
【０１０２】
次いで、ステップＳ４、ステップＳ１０の後に、ステップＳ５又はステップＳ１０において、縮退細胞候補除去部１１が、縮退した細胞候補の除去を行う。時系列および細胞形状を求める各ステップについてステップ後に細胞候補の除去を試みる。この際のアルゴリズムとして、縮退した細胞候補を含んだ形で計算を進めても良いが、縮退した細胞候補に関しても通常の細胞候補と同様の計算量がかかるため、適切に早めに除去していくことが望ましい。このため、maxac×maxvc<10-4σの場合に細胞候補を除去すると良い。
【０１０３】
また必要に応じて、スパイク時系列／時間方向のベースライン演算部５と、細胞形状／空間方向のベースライン演算部７と、の出力の少なくともいずれか一方に対して、類似した細胞候補を統合する類似細胞統合部１３を任意に設けても良い（ステップＳ６，ステップＳ１１）。
【０１０４】
類似細胞統合部１３は、複数の細胞候補の時系列の比較を行い十分に類似性が高いと判定した場合に、細胞形状を足し合わせた細胞候補を新たに導入し、元の細胞候補を除去することにより、同一細胞が複数の細胞候補に分割されることを防ぐことができる。
【０１０５】
ただし、関連度自動決定の効果によりこのような統合処理を必要としない場合もあるため、関連度自動決定の効果に依存して、これを設けるかどうかを決めれば良い。
【０１０６】
次いで、ステップＳ７又はステップＳ１２において、反復・収束判定部１５が、収束判定を行う。目的関数の変化が十分小さくかつ除去された細胞候補がない場合に収束したと判定して終了する（ステップＳ８、Ｓ１３）。ステップＳ７で収束していないと判定されると、ステップＳ９に進み、ステップＳ１２で収束していないと判定されると、ステップＳ４に戻る。
【０１０７】
尚、収束の判定は、例えば、細胞がない場合の二乗誤差をσ02、目的関数を０．５σ02としたとき、各ステップでの目的関数の変化が、閾値０．００１σ02以下であれば収束したとして終了することができる。この閾値の値は、経験則で決めることができるが、この方法に限定されるものではない。
【０１０８】
以上の処理を繰り返すことで収束すると、細胞数Cと変数ac,abaseline,vc,vbaselineおよびucすべてが全て求まり、これらの値を出力部１７が結果を出力することができる。
【０１０９】
上記の実施の形態において、添付図面に図示されている構成等については、これらに限定されるものではなく、本発明の効果を発揮する範囲内で適宜変更することが可能である。その他、本発明の目的の範囲を逸脱しない限りにおいて適宜変更して実施することが可能である。また、本発明の各構成要素は、任意に取捨選択することができ、取捨選択した構成を具備する発明も本発明に含まれるものである。
【０１１０】
また、本実施の形態で説明した機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより各部の処理を行ってもよい。尚、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。
【産業上の利用可能性】
【０１１１】
本発明は、情報処理装置として利用可能である。
【符号の説明】
【０１１２】
Ａ…情報処理装置、１…データ読み込み部、３…初期設定部、５…スパイク時系列／時間ベースライン演算部、７…細胞形状／空間ベースライン演算部、１１…縮退細胞除去部、１５…反復・収束判定部、１７…出力部、２１…解析モデル記憶部。
【図１】
2014095678.tif 000038
【図２Ａ】
2014095678.tif 000039
【図２Ｂ】
2014095678.tif 000040
【図３】
2014095678.tif 000041
【図４】
2014095678.tif 000042
【図５】
2014095678.tif 000043