(57)【特許請求の範囲】
【請求項１】
各細胞ｋの活動を、非負の値で表現される各位置ｌの細胞形状aklと非負の値で表現されるスパイク時系列uktより導出されるカルシウム濃度変化vktとの積により表現し、観測データFltを複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データにおける細胞の位置と活動を求める情報処理装置であって、
細胞のスパイク生成に起因する細胞内のカルシウム濃度の上昇に伴って起こる蛍光タンパク質からの発光を位置ｌにおける時刻tでの蛍光強度Fltとして読み取るデータ読み取り部と、
スパイク時系列uktとカルシウム濃度変化vktとの関係式のパラメータを設定するパラメータ設定部と、
前記データ読み取り部により読み取られた位置lにおける時刻tでの蛍光強度Fltに基づいて細胞形状aklの細胞候補を配置する初期設定部と、
予め定められた最小化問題の目的関数であって、かつ、細胞形状aklに関するペナルティ項、及びスパイク時系列uktに関するペナルティ項を含む目的関数を減少させるように、細胞形状akl、カルシウム濃度変化vktおよびベースラインをより真の解に近い推定値に更新する演算部と、
前記演算部における更新で目的関数の変化が十分小さくなった場合に目的関数が収束したと判定して処理を終了する収束判定部と、
を含み、
前記演算部は、
前記スパイク時系列及び前記ベースラインを更新するスパイク時系列ベースライン演算部と、
前記細胞形状及び前記ベースラインを更新する細胞形状ベースライン演算部と、
前記スパイク時系列ベースライン演算部及び前記細胞形状ベースライン演算部の出力の少なくとも一方に対して、前記目的関数の値を最小化するように、細胞候補を統合する細胞統合部と、
を含み、
前記細胞統合部は、空間的に近い細胞候補のペア（ｋ、ｋ’）について、
0006787561.tif 000073

が最小となる、細胞候補のペア（ｋ、ｋ’）を統合した細胞候補ｋ”のｓk”ａkl”ｖkt”を求め、細胞候補のペア（ｋ、ｋ’）を除去して細胞候補ｋ”を追加したときの前記目的関数の値を評価して、前記細胞候補ｋ”への統合を採用するか否かを決定し、
前記モデルは、以下の式で表わされる情報処理装置。
【数１】
0006787561.tif 000074
ただし、ｓkは、前記演算部により推定される、各細胞ｋの活動の強さを表す変数であり、ｂltは、時刻ｔでの各位置ｌのベースラインを表し、γτは、前記関係式のパラメータであるスパイク信号の線形インパルス応答を表し、εltは、時刻ｔでの各位置ｌのノイズを表し、Ｅは、前記目的関数に含まれる二乗誤差項であり、ＰAは、前記細胞形状に関するペナルティ項であり、ＰUは、前記スパイク時系列に関するペナルティ項であり、Ｋは、細胞の数を表し、Ｌは、位置ｌの総数を表し、Ｔは、時刻ｔの総数を表す。
【請求項２】
前記演算部は、ノイズレベルを算出し、
前記目的関数に含まれる前記二乗誤差項の強さが、前記ノイズレベルを用いて調整される請求項１記載の情報処理装置。
【請求項３】
前記ベースラインは、空間ベースライン及び時間ベースラインであり、
前記目的関数は、空間ベースラインに関するペナルティ項、及び時間ベースラインに関するペナルティ項を更に含む請求項１又は２記載の情報処理装置。
【請求項４】
前記演算部は、ノイズレベルを算出し、
前記目的関数に含まれる前記空間ベースラインに関するペナルティ項、及び前記時間ベースラインに関するペナルティ項の強さが、前記ノイズレベルを用いて調整される請求項３記載の情報処理装置。
【請求項５】
前記目的関数は、以下の式で表わされる請求項４記載の情報処理装置。
【数２】
0006787561.tif 000075
0006787561.tif 000076
0006787561.tif 000077
ただし、λLは、前記空間ベースラインに関するペナルティ項の強さを表し、λTは、前記時間ベースラインに関するペナルティ項の強さを表し、λAは、前記細胞形状に関するに関するペナルティ項の強さを表し、λUは、前記スパイク時系列に関するペナルティ項の強さを表し、￣ｂは、ベースライン全体の平均を表し、^ｂlは、位置ｌの空間ベースラインを表し、~ｂtは、時刻ｔの時間ベースラインを表す。
【請求項６】
前記目的関数は、以下の式で表わされる請求項５記載の情報処理装置。
【数３】
0006787561.tif 000078
ただし、σ2は、ノイズレベルを表す。
【請求項７】
前記初期設定部は、想定される細胞のサイズより大きめの特定の領域を覆うような細胞候補を十分な重なりを持たせて画像全体を埋め尽くすように真の細胞の位置に無関係に多数配置することを特徴とする請求項１〜請求項６のいずれか１項に記載の情報処理装置。
【請求項８】
前記演算部は、
前記スパイク時系列及び前記ベースラインを更新するスパイク時系列ベースライン演算部と、
前記細胞形状及び前記ベースラインを更新する細胞形状ベースライン演算部と、
前記スパイク時系列ベースライン演算部及び前記細胞形状ベースライン演算部の出力の少なくとも一方に対して、前記目的関数の値を最小化するように、細胞候補の除去を行う細胞候補除去部と、
を含む請求項１〜請求項７の何れか１項に記載の情報処理装置。
【請求項９】
前記演算部は、
細胞形状ａklを固定し、スパイク時系列の変数値が非負という拘束条件のもとで、MAP推定の目的関数を最小化するカルシウム濃度変化ｖktと空間ベースラインｂlと時間ベースラインｂtとを求める二次計画問題を解くスパイク時系列ベースライン演算部と、
カルシウム濃度変化ｖktを固定し、細胞形状の変数値が非負という拘束条件のもとで、MAP推定の目的関数を最小化する細胞形状ａklと空間ベースラインｂlと時間ベースラインｂtとを求める二次計画問題を解く細胞形状ベースライン演算部と、
を含む請求項１〜請求項８のいずれか１項に記載の情報処理装置。
【請求項１０】
前記スパイク時系列ベースライン演算部において、
当該二次計画問題を逐次二次計画法と主双対内点法とを組み合わせて解くことを特徴とする
請求項９に記載の情報処理装置。
【請求項１１】
前記スパイク時系列ベースライン演算部において、
前記逐次二次計画法の局所問題で更新される係数の計算を、主双対内点法の反復ステップ内で行うことを特徴とする請求項１０に記載の情報処理装置。
【請求項１２】
前記細胞形状ベースライン演算部において、
当該二次計画問題を逐次二次計画法と主双対内点法とを組み合わせて解くことを特徴とする
請求項９〜請求項１１の何れか１項に記載の情報処理装置。
【請求項１３】
前記細胞形状ベースライン演算部において、
前記逐次二次計画法の局所問題で更新される係数の計算を、主双対内点法の反復ステップ内で行うことを特徴とする請求項１２に記載の情報処理装置。
【請求項１４】
各細胞ｋの活動を、非負の値で表現される各位置ｌの細胞形状aklと非負の値で表現されるスパイク時系列uktより導出されるカルシウム濃度変化vktとの積により表現し、観測データFltを複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データにおける細胞の位置と活動を求める情報処理装置における情報処理方法であって、
データ読み取り部が、細胞のスパイク生成に起因する細胞内のカルシウム濃度の上昇に伴って起こる蛍光タンパク質からの発光を位置ｌにおける時刻tでの蛍光強度Fltとして読み取り、
パラメータ設定部が、スパイク時系列uktとカルシウム濃度変化vktとの関係式のパラメータを設定し、
初期設定部が、前記データ読み取り部により読み取られた位置lにおける時刻tでの蛍光強度Fltに基づいて細胞形状aklの細胞候補を配置し、
演算部が、予め定められた最小化問題の目的関数であって、かつ、細胞形状aklに関するペナルティ項、及びスパイク時系列uktに関するペナルティ項を含む目的関数を減少させるように、細胞形状akl、カルシウム濃度変化vktおよびベースラインをより真の解に近い推定値に更新し、
収束判定部が、前記演算部における更新で目的関数の変化が十分小さくなった場合に目的関数が収束したと判定して処理を終了する
ことを含み、
前記演算部が更新することは、
スパイク時系列ベースライン演算部が、前記スパイク時系列及び前記ベースラインを更新することと、
細胞形状ベースライン演算部が、前記細胞形状及び前記ベースラインを更新することと、
細胞統合部が、前記スパイク時系列ベースライン演算部及び前記細胞形状ベースライン演算部の出力の少なくとも一方に対して、前記目的関数の値を最小化するように、細胞候補を統合することと、
を含み、
前記細胞統合部が統合することでは、空間的に近い細胞候補のペア（ｋ、ｋ’）について、
0006787561.tif 000079

が最小となる、細胞候補のペア（ｋ、ｋ’）を統合した細胞候補ｋ”のｓk”ａkl”ｖkt”を求め、細胞候補のペア（ｋ、ｋ’）を除去して細胞候補ｋ”を追加したときの前記目的関数の値を評価して、前記細胞候補ｋ”への統合を採用するか否かを決定し、
前記モデルは、以下の式で表わされる情報処理方法。
【数４】
0006787561.tif 000080
ただし、ｓkは、前記演算部により推定される、各細胞ｋの活動の強さを表す変数であり、ｂltは、時刻ｔでの各位置ｌのベースラインを表し、γτは、前記関係式のパラメータであるスパイク信号の線形インパルス応答を表し、εltは、時刻ｔでの各位置ｌのノイズを表し、Ｅは、前記目的関数に含まれる二乗誤差項であり、ＰAは、前記細胞形状に関するペナルティ項であり、ＰUは、前記スパイク時系列に関するペナルティ項であり、Ｋは、細胞の数を表し、Ｌは、位置ｌの総数を表し、Ｔは、時刻ｔの総数を表す。
【請求項１５】
コンピュータを、請求項１〜請求項１３の何れか１項に記載の情報処理装置の各部として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、情報処理装置、方法、及びプログラムに係り、特に、神経細胞の活動電位を画像処理等により推測する情報処理装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
従来より、カルシウムイメージングデータからの細胞活動の推定手法の多くは、ROIと呼ばれる細胞形状を求めてからROIの平均カルシウム濃度の変化を計算し、スパイク時刻を推定するものである。しかし、ROIの重なりがある場合などには適切な方法と言えない。
【０００３】
それを改善する手法としては非特許文献1による比較的簡便なヒューリスティックな手法と、非特許文献２〜４、特許文献１のようなモデルベースの手法とがある。モデルベースの手法の多くは非負拘束条件を積極的に利用している。しかし、カルシウムイメージングデータの観測値が相対値であることから、非負拘束条件を正しく適用するためには細胞の活動と無関係なベースラインを正確に推定する必要があり、実験状況によってベースラインも時間的あるいは空間的に異なるので注意が必要である。
【０００４】
非特許文献３は非負拘束条件を用いた最も直接的な手法であるが、ベースラインの推定を別に行う必要がある。特許文献１、非特許文献４は非特許文献２を発展させた手法でありベースラインと非負拘束条件での細胞活動の推定を同時に行う手法である。特許文献１、非特許文献４には類似点も多いがそれぞれに欠点がある。
【０００５】
以下、特許文献１、非特許文献４について説明する。
【０００６】
（確率モデル）
【０００７】
観測値Ｆltは、ベースラインｂlt、細胞の数Ｋ、細胞の活動ａkl、ｖkt、ノイズεltにより構成される。各細胞の活動は空間と時間の成分ａkl、ｖktに分離できる。時間成分は、スパイク時系列ｕktとAutoregressive(AR)modelで表わされる。ａkl、ｖktはそれぞれ非負である。
【０００８】
【数１】
0006787561.tif 000002
・・・（１）
【０００９】
なお、各細胞の活動が同じでも時間成分と空間成分は、定数ｓに対して、（ａl、ｖt）と（ｓａl、ｓ-1ｖt）が同じａklｖktとなり一意に決まらない。
【００１０】
（二乗誤差）
パラメータの評価に二乗誤差を用いる。
【００１１】
【数２】
0006787561.tif 000003
・・・（２）
【００１２】
（スパース性（ペナルティ項）の導入）
二乗誤差のみの評価では、細胞数を増やせばよいことになってしまうので、適切なペナルティ項を導入する必要がある。また、細胞は局所的でスパイクも全か無かの性質を持つためスパースであることを仮定するのが自然である。よって、L1ノルムによるペナルティ項を設定する。
【００１３】
【数３】
0006787561.tif 000004
・・・（３）
【００１４】
（問題設定）
二乗誤差とペナルティ項を共に小さくするパラメータを求める。
【００１５】
（特許文献１の方法）
二乗誤差とペナルティ項の和について最小化を行う。
【００１６】
【数４】
0006787561.tif 000005
・・・（４）
【００１７】
細胞活動は自動的に以下を満たす必要があり一意に決定される。
【００１８】
【数５】
0006787561.tif 000006
・・・（５）
【００１９】
特許文献１では、時間と空間に関する二次計画問題を交互に反復して解き全体の問題を解いている。しかし、反復の途中では上記の条件を加えた形で解くことができないため収束が保障されない。
【００２０】
（非特許文献４の方法）
ピクセル毎の二乗誤差を拘束条件として、ペナルティ項を最小化する。拘束条件ではAR modelにより推定したピクセルｌごとのノイズレベルσl2を利用する。
【００２１】
【数６】
0006787561.tif 000007
・・・（６）
【００２２】
最適化の過程で時間方向と空間方向が独立して計算可能なようにｂlt＝ｂlｂtで非負と仮定する。このとき、（ａl、ｖt）は一意にならないが、ａlｖtは一意に決まる。
【００２３】
ベースラインと問題設定の単純化によって時間と空間の問題を交互に解いて収束することを保障している。また、一般に二次計画問題を解くよりも高速に数値解を得ることができる。ただし、計算途中で推定値を評価することは難しい。
【００２４】
（ベースライン）
特許文献１ではベースラインとして、空間成分と時間成分の和を用いている。この場合、全体の最適化問題を空間と時間にうまく分離できないことに起因する困難さがある。
【００２５】
【数７】
0006787561.tif 000008
・・・（７）
【先行技術文献】
【特許文献】
【００２６】
【特許文献１】特開２０１４−０９５６７８
【非特許文献】
【００２７】
【非特許文献１】E.A. Mukamel, A. Nimmerjahn and M.J. Schnitzer, “Automated analysis of cellular signals from large-scale calcium imaging data”, Neuron 63, 747-760, 2009.
【非特許文献２】J.T. Vogelstein, at. el., “Fast nonnegative doconvolution for spike train interface from population calcium imaging, J. Neurophysiol 104, 3691-3704, 2010.
【非特許文献３】R. Maruyama, at. el., “Detecting cells using non-negative matrix factorization on calcium imaging data”, Neural Network 55, 11-19, 2014.
【非特許文献４】E.A. Pnevmatikakis, at. el., “Simultaneous denoising, deconvolution, and demixing of calcium imaging data”, Neuron 89, 285-299, 2016.
【発明の概要】
【発明が解決しようとする課題】
【００２８】
非特許文献４では多くの近似や単純化を行っているため、精度が十分でない。特許文献１の方法では数値計算の負荷が高く、複数のモデルパラメータの調整が必要である。また、いずれの場合も細胞ごとのSN比に大きな差がある場合に不適切な細胞数を出力する。また、ベースラインは誤差の計算上は細胞成分とほぼ同じ扱いであるにも関わらずペナルティが科せられないため、本来細胞として検出されるべき成分の一部がベースラインとして扱われ検出されない場合がある。
【００２９】
本発明は、上記の事情を鑑みてなされたもので、観測データにおける細胞の位置と活動を精度よく推定することができる情報処理装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００３０】
上記の目的を達成するために本発明に係る情報処理装置は、各細胞ｋの活動を、非負の値で表現される各位置ｌの細胞形状aklと非負の値で表現されるスパイク時系列uktより導出されるカルシウム濃度変化vktとの積により表現し、観測データFltを複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データにおける細胞の位置と活動を求める情報処理装置であって、細胞のスパイク生成に起因する細胞内のカルシウム濃度の上昇に伴って起こる蛍光タンパク質からの発光を位置ｌにおける時刻tでの蛍光強度Fltとして読み取るデータ読み取り部と、スパイク時系列uktとカルシウム濃度変化vktとの関係式のパラメータを設定するパラメータ設定部と、前記データ読み取り部により読み取られた位置lにおける時刻tでの蛍光強度Fltに基づいて細胞形状aklの細胞候補を配置する初期設定部と、予め定められた最小化問題の目的関数であって、かつ、細胞形状aklに関するペナルティ項、及びスパイク時系列uktに関するペナルティ項を含む目的関数を減少させるように、細胞形状akl、カルシウム濃度変化vktおよびベースラインをより真の解に近い推定値に更新する演算部と、前記演算部における更新で目的関数の変化が十分小さくなった場合に目的関数が収束したと判定して処理を終了する収束判定部と、を含み、前記モデルは、以下の式で表わされる。
【００３１】
【数８】
0006787561.tif 000009
【００３２】
ただし、ｓkは、各細胞ｋの強さを表す変数であり、ｂltは、時刻ｔでの各位置ｌのベースラインを表し、γτは、前記関係式のパラメータであるスパイク信号の線形インパルス応答を表し、εltは、時刻ｔでの各位置ｌのノイズを表し、Ｅは、前記目的関数に含まれる二乗誤差項であり、ＰAは、前記細胞形状に関するペナルティ項であり、ＰUは、前記スパイク時系列に関するペナルティ項であり、Ｋは、細胞の数を表し、Ｌは、位置ｌの総数を表し、Ｔは、時刻ｔの総数を表す。
【００３３】
本発明に係る情報処理方法は、各細胞ｋの活動を、非負の値で表現される各位置ｌの細胞形状aklと非負の値で表現されるスパイク時系列uktより導出されるカルシウム濃度変化vktとの積により表現し、観測データFltを複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データにおける細胞の位置と活動を求める情報処理装置における情報処理方法であって、データ読み取り部が、細胞のスパイク生成に起因する細胞内のカルシウム濃度の上昇に伴って起こる蛍光タンパク質からの発光を位置ｌにおける時刻tでの蛍光強度Fltとして読み取り、パラメータ設定部が、スパイク時系列uktとカルシウム濃度変化vktとの関係式のパラメータを設定し、初期設定部が、前記データ読み取り部により読み取られた位置lにおける時刻tでの蛍光強度Fltに基づいて細胞形状aklの細胞候補を配置し、演算部が、予め定められた最小化問題の目的関数であって、かつ、細胞形状aklに関するペナルティ項、及びスパイク時系列uktに関するペナルティ項を含む目的関数を減少させるように、細胞形状akl、カルシウム濃度変化vktおよびベースラインをより真の解に近い推定値に更新し、収束判定部が、前記演算部における更新で目的関数の変化が十分小さくなった場合に目的関数が収束したと判定して処理を終了することを含み、前記モデルは、上記の式で表わされる。
【００３４】
本発明に係るプログラムは、コンピュータを、上記の情報処理装置の各部として機能させるためのプログラムである。
【発明の効果】
【００３５】
以上説明したように、本発明の情報処理装置、方法、及びプログラムによれば、細胞形状aklに関するペナルティ項、及びスパイク時系列uktに関するペナルティ項を含む目的関数を減少させるように、細胞形状akl、カルシウム濃度変化vktおよびベースラインの推定値を更新することにより、観測データにおける細胞の位置と活動を精度よく推定することができる、という効果が得られる。
【図面の簡単な説明】
【００３６】
【図１】本発明の実施の形態に係る情報処理装置を示すブロック図である。
【図２】本発明の実施の形態に係る情報処理装置を示す機能ブロック図である。
【図３】逐次２次計画法と内点法とを組み合わせた場合の通常の非線形問題の解法の処理の流れを示すフローチャートである。
【図４】本発明の実施の形態に係る情報処理装置の解析処理ルーチンの内容を示すフローチャートである。
【図５】細胞形状を固定してスパイク時系列及びベースラインを求める処理、及びスパイク時系列を固定して細胞形状及びベースラインを求める処理の流れを示すフローチャートである。
【発明を実施するための形態】
【００３７】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００３８】
＜本発明の実施の形態の概要＞
（ペナルティ項の変更（細胞毎の正規化））
以下の式に示すように、細胞の活動の強さを表す変数ｓkを導入することで、細胞ごとのSN比に大きな差がある場合に不適切な細胞数を出力することへの対応を行う。また、特許文献１では、時間方向と空間方向の問題を交互に反復して解くことで全体の解に収束させる方法を採用していたが、時間成分と空間成分の依存関係により実際には収束しない場合があった。細胞形状に関するペナルティ項及びスパイク時系列に関するペナルティ項を導入することにより、独立性が保障され必ず収束する。なお、細胞形状に関するペナルティ項ＰA及びスパイク時系列に関するペナルティ項ＰUの評価が変わるので、注意する。
【００３９】
【数９】
0006787561.tif 000010
・・・（８）
【００４０】
（パラメータの自動チューニング（最尤推定、MAP推定））
誤差εltが正規分布に従うことを仮定して、以下の式に示すように、ノイズレベルσ2を変数として評価関数に導入する。最尤法を用いるとペナルティの強さが自動的に調整されるのと同じ効果となる。また、a0は細胞の形状を表す画像の平均値、u0は平均発火率を表すパラメータで、a,uはそれぞれ指数分布に従うとする。
【００４１】
【数１０】
0006787561.tif 000011
・・・（９）
【００４２】
（ベースラインへのペナルティ項の導入）
ベースラインの変化は例外的であり基本的に0であることが自然であることを踏まえて、ベースラインの分散にペナルティを導入する。このことにより、ベースライン成分が過剰になることを抑えて細胞がより正しく検出されるようになる。
【００４３】
この際、ベースラインは撮影状況による空間的なベースラインのばらつきと細胞の代謝などの影響による時間的なベースラインの変化の和となるとする。カルシウムイメージングデータでは観測値がカルシウム濃度の相対値として得られるため、非特許文献４の積のモデルよりも自然であり、ベースラインへのペナルティの導入も容易である。
【００４４】
【数１１】
0006787561.tif 000012
・・・（１０）
【００４５】
（パラメータの自動チューニングとベースラインのペナルティ項との組み合わせ）
ベースラインの変動はノイズレベルに比例し、正規分布に従うという事前分布を導入すると以下のようになる。
【００４６】
【数１２】
0006787561.tif 000013
・・・（１１）
【００４７】
（安定で高速な解法）
また、ベースラインのペナルティ項の導入は、単にモデルを改良するだけでなく、内点法で計算する連立一次方程式の対角成分が安定化するため、内点法の収束が加速し解が安定化する副作用ももたらす。
【００４８】
細胞形状及びスパイク時系列に関するペナルティ項と、ベースラインのペナルティ項とを導入し、時間ステップ、空間ステップそれぞれにおいて、ベースライン全体を考慮することによって、特許文献１における全体の収束に関する問題を解決する。パラメータの自動チューニングとベースラインのペナルティ項との組み合わせの導入により、ペナルティ項の係数がノイズレベルに応じて自動的に決定されることで係数のチューニングが不要になる。
【００４９】
計算上の工夫（aではなく，saを変数と見て計算する）ことにより、非線形の効果が限定的にしか現れず、計算量を増やすことなくモデルの改良の効果を得ることができる。
【００５０】
しかし、一方で各ステップの最適化問題に非線形成分が含まれることになり、一般には解くことが難しくなる。そこで、逐次二次計画法と主双対内点法を組み合わせて工夫を行うことで、特許文献１の方法とほぼ同じ計算量で解を求めることができる方法を開発した。具体的には逐次二次計画法の局所問題で更新される計算を内点法の反復ステップ内で行ってしまう。この方法は、一般に利用できるわけではないが、式変形により逐次2次計画法の問題の更新で更新される係数が限定されていることを利用して可能にしている。
【００５１】
（細胞数の決定）
上記（１１）式で一元的に適切な細胞数を評価できる。上記（１１）式を直接評価することは計算コストが高いが、中間の計算結果を再利用して効率的に評価が可能である。
【００５２】
得られた推定値のセットから、特定の細胞候補を除去してベースラインのみを再計算する。上記（１１）式により除去前と除去後の評価を行い、細胞を除去するかを決定する。
【００５３】
また、空間的に近い細胞候補ｋ、ｋ´について
【数１３】
0006787561.tif 000014
が最小となるｓｋ'', ｓｋl''，ｖｋt''を求める。ｋ、ｋ'を除去してｋ''を追加し、ベースラインを再計算する（統合）。上記（１１）式を評価し統合を採用するかを決定する。
【００５４】
また、ａklに複数のクラスタが存在するときには，細胞を分割する。
【００５５】
＜システム構成＞
図１に示すように、本発明の実施の形態に係る情報処理装置１０は、ＣＰＵ１２、ＲＯＭ１４、ＲＡＭ１６、ＨＤＤ１８、通信インタフェース２０、及びこれらを相互に接続するためのバス２２を備えている。
【００５６】
ＣＰＵ１２は、各種プログラムを実行する。ＲＯＭ１４には、各種プログラムやパラメータ等が記憶されている。ＲＡＭ１６は、ＣＰＵ１２による各種プログラムの実行時におけるワークエリア等として用いられる。記録媒体としてのＨＤＤ１８には、後述する解析処理ルーチンを実行するためのプログラムを含む各種プログラムや各種データが記憶されている。
【００５７】
本実施の形態における情報処理装置１０は、各細胞ｋの活動を、非負の値で表現される各位置ｌの細胞形状aklと非負の値で表現されるスパイク時系列uktより導出されるカルシウム濃度変化vktとの積により表現し、観測データFltを複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データにおける細胞の位置と活動を求める。
【００５８】
本実施の形態における情報処理装置１０は、機能的には、図２に示すように、データ読み込み部３０、初期設定部３２、スパイク時系列ベースライン演算部３４、細胞処理部３６、細胞形状ベースライン演算部３８、及び出力部４２を備えている。なお、初期設定部３２が、パラメータ設定部及び初期設定部の一例である。
【００５９】
また、細胞処理部３６は、縮退細胞候補除去部５０、類似細胞統合部５２、及び反復・収束判定部５４を備えている。
【００６０】
＜本発明の実施の形態の原理＞
次に、本実施の形態のモデルと数値解法の原理について説明する。
【００６１】
本実施の形態では、以下の問題を解く。
【００６２】
（確率モデル）
【数１４】
0006787561.tif 000015
・・・（１２）
（等式拘束条件）
【数１５】
0006787561.tif 000016
・・・（１３）
（不等式拘束条件）
【数１６】
0006787561.tif 000017
・・・（１４）
（最適化問題）
【数１７】
0006787561.tif 000018
・・・（１５）
【００６３】
ただし、ｓkは、各細胞ｋの活動の強さを表す変数であり、γτは、スパイク信号の線形インパルス応答を表し、εltは、時刻ｔでの各位置ｌのノイズを表し、Ｅは、目的関数に含まれる二乗誤差項であり、ＰAは、細胞形状に関するペナルティ項であり、ＰUは、スパイク時系列に関するペナルティ項であり、Ｋは、細胞の数を表し、Ｌは、位置ｌの総数を表し、Ｔは、時刻ｔの総数を表す。また、￣ｂは、ベースライン全体の平均を表し、^ｂlは、位置ｌの空間ベースラインを表し、~ｂtは、時刻ｔの時間ベースラインを表す。また、ηLは、データのノイズレベルσに対して、想定される空間ベースラインの標準偏差がどの程度の大きさかを比で表わしたものであり、ηTは、データのノイズレベルσに対して、想定される時間ベースラインの標準偏差がどの程度の大きさかを比で表わしたものである。実験状況からベースラインが大きく変動しているようなデータでは、ηL、ηTの値を大きく設定する。また、a0は細胞の形状を表す画像の平均値、u0は平均発火率を表すパラメータで、a,uはそれぞれ指数分布に従うとする。
【００６４】
（モデルパラメータの決定）
まず、実験条件より、a0, u0, ηL2, ηT2 を設定する。
【００６５】
また、γτは非特許文献４と同様にAR modelを用いてデータから決定できる。
【００６６】
（手順）
手順１．  a0, u0, ηL2, ηT2を設定する。
【００６７】
手順２．  AR model を利用して γτを求める。
【００６８】
手順３．  初期細胞候補数K{0}を決め初期値
0006787561.tif 000019
をランダムに作成する。このとき、Fltの性質を利用してなるべく良い性質の初期値を設定する。
【００６９】
手順４．
0006787561.tif 000020
を固定した問題（詳細を後述する）を解く。
【００７０】
0006787561.tif 000021
が求まる。
【００７１】
手順５．
0006787561.tif 000022
を固定した問題（詳細を後述する）を解く。
【００７２】
0006787561.tif 000023
が求まる．
【００７３】
手順６．
0006787561.tif 000024
をもとに細胞候補の削除、統合、分割を行う。
【００７４】
0006787561.tif 000025
が求まる。
【００７５】
手順７．  収束条件を確認し、収束したら終了する。収束していなければ手順４に戻る。
【００７６】
（v{i}を固定した問題を解く）
（定式化）
データにのみ依存する係数部分を以下のように整理する。
【００７７】
【数１８】
0006787561.tif 000026
【００７８】
また、ｓkａklに関する形式的な変形により、sをaに従属させる。
【００７９】
【数１９】
0006787561.tif 000027
・・・（１６）
【００８０】
また、以下のようにｖに依存している部分の係数を整理する。
【００８１】
【数２０】
0006787561.tif 000028
・・・（１７）
【００８２】
解くべき問題は、以下のように表される。
【００８３】
【数２１】
0006787561.tif 000029
・・・（１８）
0006787561.tif 000030
・・・（１９）
【００８４】
上記（１８）式では、Qの一部の対角成分の項の効果により解が安定化する（対角成分が優位になる）。
【００８５】
上記（１９）式では、第１項目（ｌｏｇσ2）と、第３項目の前半部分（ｒaの部分）とが非線形である。ｒaに含まれるρが、ｘ（ａkl*）とそこから計算されるｓkに依存する。
【００８６】
（逐次二次計画法）
非線形項の扱い方で幾つかの方法が考えられる。いずれの方法であっても、ｘに依存する変数を計算し、一般の二次計画法となるように問題を変形し、変形した問題の解で再度依存する変数を計算する。上記（１９）式の特性により、ごく一部の変数を計算すればほとんど再計算なく変形された一般の二次方程式が得られる。
【００８７】
（逐次二次計画法の第１の例の手順）
簡単のために上述した２つの非線形部分について、それぞれσ2とｒa(ｓk)とが定数であるとみなすと、以下の2つのステップの繰り返しとなる。
【００８８】
ステップ１．  上記（１９）式でｘからσ2とｒa(ｓk)を求める。
【００８９】
【数２２】
0006787561.tif 000031
【００９０】
ステップ２．  σ2とｒaを定数とみなして上記（１９）式を解きｘを求める。
【００９１】
【数２３】
0006787561.tif 000032
【００９２】
（逐次二次計画法の第２の例の手順）
非線形部分のうち、ｒaの部分については微分が連続な関数とならないため、線形化が難しいが、Qの一部の対角成分の項については、（２０）式のように線形化が可能である。ここでもｒaの部分については、ｘを更新した後に、上記（１６）式により、ｓkが更新され、（１８）式によりｒa(ｓk)が更新される。
【００９３】
【数２４】
0006787561.tif 000033
・・・（２０）
0006787561.tif 000034
・・・（２１）
【００９４】
具体的には、以下の４つのステップで実現する。
【００９５】
ステップ１．  初期値ｘ{0}*を定める。
【００９６】
ステップ２．  ｘ{i}*を用いて上記（２１）式に従ってｓ{i}を求める。
【００９７】
ステップ３．  ｘ{i}*、ｓ{i}を用いて、Ｑ{i}*、ｒ{i}*を計算する。
【００９８】
ステップ４．  Ｑ{i}*、ｒ{i}*を用いた、以下の（２２）式で表わされる問題を解き、ｘ{i+1}を求める。
【００９９】
（主双対問題と内点法）
（主双対問題）
局所二次計画問題（上記（２１）式）の解は以下の問題の解と等しい。
【０１００】
【数２５】
0006787561.tif 000035
・・・（２２）
【０１０１】
（内点法）
予測子修正子法を導入した内点法による解法を、以下のステップ１〜ステップ９に示す。本発明の内容は内点法の詳細には依らない。
【０１０２】
ステップ１．  初期値
0006787561.tif 000036
を設定する。
【０１０３】
ステップ２．  以下の式を計算する。
【０１０４】
【数２６】
0006787561.tif 000037
【０１０５】
ステップ３．
0006787561.tif 000038
から、（２３）式により
0006787561.tif 000039
を求める。大規模データでも行列が疎で構造があるため共役勾配法などの反復法で効率よく解が求まる。
【０１０６】
【数２７】
0006787561.tif 000040
・・・（２３）
【０１０７】
ステップ４．  定数
0006787561.tif 000041
を求める。
【０１０８】
ステップ５．
0006787561.tif 000042
から（２４）式により、
0006787561.tif 000043
を求める。
【０１０９】
【数２８】
0006787561.tif 000044
・・・（２４）
【０１１０】
上記（２４）式は、（２３）式と比べると、右辺がわずかに異なっている。
【０１１１】
ステップ６．  非負条件を満たすステップ幅αを求める。
【０１１２】
ステップ７．  評価値
0006787561.tif 000045
で収束を確認する。
【０１１３】
ステップ８．α’＝αsafeとして、以下の式で、ｘを更新する。
【０１１４】
【数２９】
0006787561.tif 000046
【０１１５】
ステップ９．  収束していない場合には、ステップ２へ戻る。
【０１１６】
（計算上の工夫）
上述したパラメータの自動チューニングを導入しても逐次2次計画問題における更新は、結局元の問題のペナルティ項の大きさを自動調整する部分しか変更されない。このため、内点法の内部で問題の更新を行っても十分安定に解くことができる。このような特性を考慮して逐次二次計画法と二次計画問題の反復解法（内点法）を組み合わせると、通常、図３のように２重ループとなるところを、1ステップごとに非線形の効果部分の係数を計算し直すことで大幅にトータルの計算量を削減することができる。
【０１１７】
つまり、ａを求める処理の手順は、以下のステップ１〜ステップ１１のようになる。
【０１１８】
ステップ１．  （１８）式に従って、Ｑ、ｒを計算する。
【０１１９】
ステップ２．  初期値
0006787561.tif 000047
を設定する。
【０１２０】
例えば、
0006787561.tif 000048
と設定すればよい。
【０１２１】
ステップ３．  （２０）式に従って、ｘ{i}*,ｓ{i}を用いて Ｑ{i}*,ｒ{i}*を計算する。
【０１２２】
ステップ４．  以下の式を計算する。
【０１２３】
【数３０】
0006787561.tif 000049
【０１２４】
ステップ５．
0006787561.tif 000050
から、（２５）式により
0006787561.tif 000051
を求める。大規模データでも行列が疎で構造があるため共役勾配法などの反復法で効率よく解が求まる。
【０１２５】
【数３１】
0006787561.tif 000052
・・・（２５）
【０１２６】
ステップ６．  定数
0006787561.tif 000053
を求める。
【０１２７】
ステップ７．
0006787561.tif 000054
から（２６）式により、
0006787561.tif 000055
を求める。
【０１２８】
【数３２】
0006787561.tif 000056
・・・（２６）
【０１２９】
上記（２６）式は、（２５）式と比べると、右辺がわずかに異なっている。
【０１３０】
ステップ８．  非負条件を満たし、かつ、ｓの変動が一定以下となるように、ステップ幅αを求める。
【０１３１】
ステップ９．  評価値
0006787561.tif 000057
で収束を確認する。
【０１３２】
ステップ１０．  α’＝αsafeとして、以下の式で、ｘを更新する。
【０１３３】
【数３３】
0006787561.tif 000058
【０１３４】
ステップ１１．  収束していない場合には、ステップ２へ戻る。
【０１３５】
ペナルティ項の導入は、最大値の選択という比較的一般的な解法への適応が難しい問題である。ただし、今回の問題では内点法での更新ステップ幅を選択する際に変数が大きく変化しすぎないように修正することで対応できる。実際には通常の内点法よりもステップ幅が小さく修正されることはほとんどなく実行速度に影響はない。
【０１３６】
上述した原理に従って、本実施の形態に係るスパイク時系列ベースライン演算部３４は、上記ステップ１〜ステップ１１を実行することにより、スパイク時系列及びベースラインを更新する。
【０１３７】
（v{i}を固定した問題を解き評価する）
ｘにｕとｖを含む。Ｃにｕとｖの関係を拘束条件として入れる。Ｄはｘからｕ部分を取り出す変換となる。Ｑ,ｒ,Ｃ,Ｄを定式化した後の解法は、上述した、v{i}を固定した問題を解く場合と同じである。
【０１３８】
（定式化）
以下のように、ｓkｕktに関する変形と係数の整理を行う。
【０１３９】
【数３４】
0006787561.tif 000059
・・・（２７）
0006787561.tif 000060
・・・（２８）
【０１４０】
解くべき問題は、以下である。
【０１４１】
【数３５】
0006787561.tif 000061
0006787561.tif 000062
・・・（２９）
0006787561.tif 000063
・・・（３０）
【０１４２】
（計算上の工夫）
ａを求める問題を解く場合と同様に、逐次二次計画法と主双対内点法の組み合わせをベースに、主双対内点法のループないで逐次2次計画法に相当する部分の計算も行い、2重ループを1重ループにして大幅に計算量を減らす。
【０１４３】
以上説明した原理に従って、本実施の形態に係る細胞形状ベースライン演算部３８は、スパイク時系列ベースライン演算部３４の上記ステップ１〜ステップ１１と同様のステップを実行することにより、細胞形状及びベースラインを更新する。
【０１４４】
（細胞候補の増減）
細胞候補を増減させたときに上記（１５）式の評価関数が小さくなれば変更を採用する。
【０１４５】
本実施の形態に係る縮退細胞候補除去部５０は、細胞形状ベースライン演算部３８の出力に対して、上記（１５）式の評価関数の値を最小化するように、細胞候補の除去を行う。
【０１４６】
類似細胞統合部５２は、細胞形状ベースライン演算部３８の出力に対して、上記（１５）式の評価関数の値を最小化するように、細胞候補を統合する。
【０１４７】
細胞候補を増減させた際にはベースラインの再計算が必要となる。
【０１４８】
（ａ，ｖを固定した問題を解き評価する）
ベースラインの再計算では、ベースラインを単純な二次計画問題として解くことができる。
【０１４９】
【数３６】
0006787561.tif 000064
・・・（３１）
【０１５０】
＜情報処理装置の動作＞
次に、本実施の形態に係る情報処理装置１０の作用について説明する。
【０１５１】
まず、情報処理装置１０に、観測データが入力されると、情報処理装置１０によって、図４に示す解析処理ルーチンが実行される。
【０１５２】
ステップＳ１００において、データ読み出し部３０が、ある時間ｔとその時間ｔにおける活動電位を発生する際における細胞内のカルシウム濃度が上昇に起因する蛍光タンパク質からの発光を二光子顕微鏡で観察した動画像データを読み込み、パラメータを設定する。読み込むデータとしては、位置ｌにおける時刻tでの蛍光強度Ｆltがデータとして与えられる。ただし、ｌ,ｔはそれぞれＬ，Ｔ以下の自然数で、サンプリングレートおよび空間解像度から決まるスケールになっているとする。尚、サンプリングレートＳ[Hz]、記録時間Ｒ[sec]とすると、Ｔ＝ＳＲである。また、Ｌは画像のピクセル数である。蛍光強度Ｆは相対的な値であるため、何らかのベースラインを規定する必要がある。
【０１５３】
次いで、ステップＳ１０２において、初期設定部３２が、ベースライン、細胞の個数、細胞の形状、及びスパイク時系列を初期設定する。例えば、細胞が存在しないと仮定した場合の最適なベースラインを初期条件として設定する。また、初期設定部３２は、ＡＲ modelを利用して、γτを求める。
【０１５４】
ここで、想定される細胞のサイズより大きめの特定の領域を覆うような細胞候補を十分な重なりを持たせて画像全体を埋め尽くすように真の細胞の位置に無関係に多数配置する。疎な解を得られるような条件でMAP推定を行うことで、不要な細胞候補は自動的に縮退し（信号が小さくなり）細胞数が自動的に決定される。
【０１５５】
次いで、ステップＳ１０４において、スパイク時系列ベースライン演算部３４が、スパイク時系列およびベースラインを求める。
【０１５６】
スパイク時系列ベースライン演算部３４は、上記（１３）式、（１４）式に示す拘束条件のもとで、上記（１５）式に示す評価関数を最小化するスパイク時系列ｕkt、カルシウム濃度変化ｖkt、空間ベースライン^ｂl、時間ベースライン~ｂtを求める二次計画問題を解く。
【０１５７】
細胞形状aklを固定すると、上記（１３）式、（１４）式に示す拘束条件のもとで、上記（１５）式に示す評価関数を最小化するスパイク時系列ｕkt、カルシウム濃度変化ｖkt、空間ベースライン^ｂl、時間ベースライン~ｂtを求める二次計画問題となる。この問題を逐次二次計画法と主双対内点法とを組み合わせて解く。
【０１５８】
そして、ステップＳ１０６において、細胞形状ベースライン演算部３８が、細胞形状およびベースラインを求める。
【０１５９】
スパイク時系列ｕktを固定すると、上記（１３）式、（１４）式の拘束条件のもとで目的関数(19)を最小化する細胞形状akl、空間ベースライン^ｂl、時間ベースライン~ｂtを求める二次計画問題となる。この問題を逐次二次計画法と主双対内点法とを組み合わせて解く。
【０１６０】
次いで、ステップＳ１０８において、縮退細胞候補除去部５０が、細胞候補の除去を行う。細胞候補の各々について、当該細胞候補を除去する前と除去した後とで、上記（１５）式の評価関数の値を算出し、当該細胞候補を除去した後の評価関数が小さくなれば、当該細胞候補を除去する。
【０１６１】
ステップＳ１１０では、類似細胞統合部５２が、細胞候補の統合を行う。具体的には、空間的に近い細胞候補のペア（ｋ、ｋ’）の各々について、
0006787561.tif 000065
が最小となる、細胞候補のペア（ｋ、ｋ’）を統合した細胞候補ｋ”のｓk”ａkl”ｖkt”を求め、細胞候補のペア（ｋ、ｋ’）を除去する前の評価関数の値と、細胞候補のペア（ｋ、ｋ’）を除去して細胞候補ｋ”を追加したときの評価関数が小さくなれば、当該細胞候補のペア（ｋ、ｋ’）を除去して、統合した細胞候補ｋ”を追加する。なお、本ステップＳ１１０は、省略してもよい。
【０１６２】
次いで、ステップＳ１１２において、反復・収束判定部５４が、収束判定を行う。目的関数の変化が十分小さくかつ除去された細胞候補がない場合に収束したと判定して、出力部４２により細胞の推定結果を出力して、終了する（ステップＳ１１４）。ステップＳ１１２で収束していないと判定されると、ステップＳ１０４に戻る。
【０１６３】
尚、収束の判定は、例えば、目的関数の変化が、閾値以下であれば収束したとして終了することができる。この閾値の値は、経験則で決めることができる。
【０１６４】
以上の処理を繰り返すことで収束すると、細胞数Ｋ、細胞形状akl、スパイク時系列ｕkt、カルシウム濃度変化ｖkt、空間ベースライン^ｂl、時間ベースライン~ｂtが全て求まり、これらの値を出力部４２が結果を出力することができる。
【０１６５】
上記ステップＳ１０４は、図５に示す処理ルーチンによって実現される。
【０１６６】
ステップＳ１２０において、（１８）式に従って、Ｑ、ｒを計算する。
【０１６７】
そして、ステップＳ１２２において、初期値
0006787561.tif 000066
を設定する。
【０１６８】
次のステップＳ１２４では、上記（２０）式に従って、ｘ{i}*,ｓ{i}を用いて Ｑ{i}*,ｒ{i}*を計算する。
【０１６９】
そして、ステップＳ１２６において、ω{ij}を計算し、
0006787561.tif 000067
から、上記（２５）式により
0006787561.tif 000068
を求める。
【０１７０】
そして、定数ω{ij}P、ρ{ij}を求める。また、
0006787561.tif 000069
から上記（２６）式により、
0006787561.tif 000070
を求める。
【０１７１】
そして、ステップＳ１２８において、非負条件を満たし、かつ、ｓの変動が一定以下となるように、ステップ幅αを求める。
【０１７２】
ステップＳ１３０では、評価値
0006787561.tif 000071
で収束を確認し、収束している場合には、処理ルーチンを終了する。一方、収束していない場合には、ステップＳ１３２において、α’＝αsafeとして、以下の式のように更新し、上記ステップＳ１２４へ戻る。
【０１７３】
【数３７】
0006787561.tif 000072
【０１７４】
ステップＳ１０６も、上記図５に示す処理ルーチンと同様の処理ルーチンにより実現される。
【０１７５】
以上説明したように、本発明の実施の形態に係る情報処理装置によれば、細胞形状aklに関するペナルティ項、スパイク時系列uktに関するペナルティ項、空間ベースラインに関するペナルティ項、及び時間ベースラインに関するペナルティ項を含む目的関数を減少させるように、細胞形状akl、カルシウム濃度変化vktおよびベースラインの推定値を更新することにより、観測データにおける細胞の位置と活動を精度よく推定することができる。
【０１７６】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１７７】
本発明のプログラムは、記憶媒体に格納して提供するようにしてもよい。
【符号の説明】
【０１７８】
１０   情報処理装置
１２   ＣＰＵ
１４   ＲＯＭ
１６   ＲＡＭ
１８   ＨＤＤ
３０   データ読み取り部
３２   初期設定部
３４   スパイク時系列ベースライン演算部
３６   細胞処理部
３８   細胞形状ベースライン演算部
４０   解析モデル記憶部
４２   出力部
５０   縮退細胞候補除去部
５２   類似細胞統合部
５４   反復・収束判定部
akl    細胞形状
ukt    スパイク時系列
vkt    カルシウム濃度変化
【図１】
0006787561.tif 000081
【図２】
0006787561.tif 000082
【図３】
0006787561.tif 000083
【図４】
0006787561.tif 000084
【図５】
0006787561.tif 000085