JP6029101B2

JP6029101B2 - 情報処理装置、情報処理プログラム

Info

Publication number: JP6029101B2
Application number: JP2012248839A
Authority: JP
Inventors: 高志竹川; 宏山川
Original assignee: Fujitsu Ltd; RIKEN Institute of Physical and Chemical Research
Current assignee: Fujitsu Ltd; RIKEN Institute of Physical and Chemical Research
Priority date: 2012-11-12
Filing date: 2012-11-12
Publication date: 2016-11-24
Anticipated expiration: 2032-11-12
Also published as: JP2014095678A

Description

本発明は、情報処理技術に関し、特に、神経細胞の活動電位を画像処理等により推測する情報処理技術に関する。

神経細胞は活動電位（スパイク）を発生した際に他の神経細胞に情報を伝達する。脳内の神経細胞集団がどのように情報処理を行っているかを理解するためには複数の神経細胞の活動電位を同時に記録する必要がある。活動電位を発生する際には細胞内のカルシウム濃度が上昇するため、カルシウム濃度に依存して光る蛍光タンパク質を導入した脳を二光子顕微鏡で動画として記録することにより、複数の細胞に活動電位が生じている様子を観察することができる。

しかしながら、実際に神経細胞集団のデータとして扱うためには動画データから、各細胞の位置と活動の時系列を抽出する必要がある。実際には、データにはノイズが多く含まれるため、得られるのは確定した情報ではなく、データを元に推定された情報である。大容量の記録データから、高速・高精度に細胞位置と活動時系列を推定することが重要である。

従来手法として主に用いられている手法は、まず主成分分析などを用いて類似した活動パターンを見せる領域を見つけ、いくつかの処理をすることで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は相対的な値であるため、何らかのベースラインを規定する必要がある。

時間軸に対するカルシウム濃度はスパイク（図３参照）によりデルタ関数的に増加し、その後、減衰する。u_c(0)=v_c(0)よりu_c,v_cは一方の時系列からもう一方が線形変換で導出可能である。

ベースラインは平均あるいはmedianで初期設定する。

v_c(t)=0として、式（４）よりσ²を設定する。

但し、細胞の個数Cおよび細胞の形状a_c(x,y)の求め方は規定されておらず、これらの初期条件を得る何らかの手順を別途設定する必要があるが、解が適切に収束するためには初期条件が解に十分近い必要がある。

次いで、ステップＳ１０４において、カルシウム濃度／スパイク時系列演算部１０５が、カルシウム濃度vおよびスパイク時系列uを求める。上記の式(1),(2),(3)において、a_c,a_baselineおよびパラメータσ,λ_cを固定して、u,vに関するMAP(事後分布最大化)推定を行うと、式（２）の条件下で次の二次計画問題を解くことになる。このとき、ベースラインが固定されており、かつ非負拘束条件が設定されているため、固定されたベースラインの影響を強く受け系全体としての適切な解が得られない。

Log障壁法とニュートン法を用いてこの最小化問題を解く。Ｃ＝１ならば３重対角行列になり比較的容易に解けるが、Ｃ＞１では３重ブロック対角行列を考える必要があり、精度よく解くことは困難である。

次いで、ステップＳ１０５において、正規化部１０７が、各細胞候補のスパイク時系列と細胞形状を正規化する処理を行う。

各細胞の信号はf_c(t,x,y)=a_c(x,y)v_c(t)と表されるが、

としても信号全体は変わらない。

しかし、式（７）は最小二乗の項とuに対するペナルティ項の和になっているため、式（２）より、s_cを大きくすればするほどペナルティ項が小さくなり目的関数の値は小さくなってしまう。これではペナルティ項の意味がないため、何らかの形で正規化する必要がある。そこで、s_c=maxv_cとする。しかし、拘束条件s_c=maxv_cをいれて問題を解いたわけではないため、正規化後は式（７）の解になっていない。

ステップＳ１０６で、細胞形状／ベースライン演算部が、細胞形状a_cおよびベースラインa_baselineを求める。ここではu,vを固定して、式（７）が最小となるa_c,a_baselineを求める。最尤法を行うと、ピクセルごとに独立して考えることができて、（９）式をそれぞれ解けばよい。

これは、拘束条件がないので簡単に解けるが、a_c(x,y)には負の値も含まれるため細胞の形状はノイズの影響を強く受け、安定しない。

ステップＳ１０７で、パラメータ更新部１１１が、パラメータの更新を行う。式（４）、（５）に従いパラメータを更新する。しかしながら、真の解で規定されるべきパラメータを現在の仮の値で近似しており、事前分布が変数から定義されることになるため、目的関数を正しく設定できない。

ステップＳ１０４とＳ１０７の処理を、反復・収束判定部１１７が、反復・収束判定する。ステップＳ１０８で処理が終了すると、出力部１２１が結果を出力する。

正規化処理の説明においても述べたように、そもそも問題設定が不良のため最終的な収束判定が不能であるため、各ステップを適当な回数繰り返したところで打ち切ることになる。

このように、非特許文献２の問題点としては、アルゴリズム中で細胞数が固定であること、初期条件として細胞の形状を正しい形状に十分近くしなければならないこと、示された解法では速度および精度が不十分であること、多数の細胞が存在する際に適切な収束が得られないこと、などである。

本発明は、画像に含まれる細胞数が不明な状態から初期条件を個別に設定することなく高速に十分な精度で安定して細胞数およびそれぞれの細胞の形状、活動時系列を推定する技術を提供することを目的とする。

本発明では、パラメータの非負拘束条件と適切な事前分布の設定による正規化項（ペナルティ項）をモデルに導入し得られる解の精度を向上させた上で、多数の細胞候補を画像中に重なりを持たせて密に埋め尽くす様な初期条件設定と、細胞候補の取捨選択の工夫により、細胞数とパラメータを正しく推定できる装置（プログラム）を開発した。主たるアルゴリズムとしては、形状と時間変化を交互に反復して改善する反復法を用い、途中でノイズレベル以下の小さな信号しか持たない細胞候補を取り除き、類似した活動状態を持つ細胞候補を統合することで計算量を削減しながら画像中に含まれる活動した細胞をすべて正しく検出する。各ステップで必要な非線形計画の解法には、主双対内点法および大規模一次方程式の高速な解法としての共役勾配法を用いた独自の並列化実装を用いる。

最終的に、画像に含まれる細胞数が不明な状態から初期条件を個別に設定することなく高速に十分な精度で安定して細胞数およびそれぞれの細胞の形状、活動時系列を推定することができるプログラムが得られた。

従来法２では、細胞数が固定かつ初期条件が厳しいという問題点がある。そこで、多数の細胞候補を画像中に重なりを持たせて密に埋め尽くす様に初期条件を設定し、解が疎になるような条件をモデルに付け加える。前者の工夫により、真の細胞はいずれかの細胞候補の十分近くにあることになるため取りこぼすことなく細胞を検出でき、かつ後者の工夫とクロストークをきちんと考慮したモデルを用いている効果により多数の細胞候補があっても必要以上に真の細胞が複数の細胞候補に分割されてしまうことを防ぐことができる。実際、アルゴリズムの反復を繰り返すと関連度自動決定と呼ばれる効果により大きな信号を持つ細胞候補と、ノイズレベル以下の小さな信号しか持たない細胞候補が生じるようになるため、小さな信号の細胞候補を取り除き、類似した活動状態を持つ細胞候補を統合することで最終的に画像中に含まれる活動した細胞をすべて正しく検出できる。ここで新しく提案した手法は、観測画像が各細胞の活動の和であるというモデルに対して適切な拘束条件を入れることで、反復法が正しく収束することをサポートしている。

この手法では、多数の細胞候補を用意するという性質上、多くのメモリおよび計算量を必要とすることになるが、画像を適切なサイズに分割して前処理しそれらの結果を統合する方針をとることで計算量を大幅に削減することができ、個々のステップでの適切な解法と並列化による高速化した実装により現実的な実行時間で目的を達成することができるようになっている。

最終的に、画像に含まれる細胞数が不明な状態から初期条件を個別に設定することなく高速に十分な精度で安定して細胞数およびそれぞれの細胞の形状、活動時系列を推定することができる様になった。

本発明の一観点によれば、細胞の活動を非負の値で表現される細胞形状a_c(x,y)と非負の値で表現されるスパイク時系列u_c(t)より導出されるカルシウム濃度変化v_c(t)の積により表現し、観測データf(t,x,y)を複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データとモデルと事前分布を元にしたベイズ理論に基づく推定により想定される細胞の位置と活動を求める情報処理装置であって、細胞のスパイク生成に起因する細胞内のカルシウム濃度の上昇に伴って起こる蛍光タンパク質からの発光を位置(x,y)における時刻tでの蛍光強度f(t,x,y)として読み取るデータ読み取り部と、スパイク時系列u_c(t)とカルシウム濃度変化v_c(t)との関係式と細胞形状a_c(x,y)とスパイク時系列u_c(t)に関する事前分布を設定するパラメータ設定部と、前記データ読み取り部により読みとられた位置(x,y)における時刻tでの蛍光強度f(t,x,y)に基づいて細胞形状a_c(x,y)の細胞候補を配置し、ベースラインを蛍光強度f(t,x,y)の統計量に基づいて設定する初期設定部と、前記ベイズ理論に基づく推定で導出される最小化問題の目的関数を減少させ、細胞形状a_c(x,y)、カルシウム濃度変化v_c(t)およびベースラインをより真の解に近い推定値に更新する演算部と、前記演算部における更新で目的関数の変化が十分小さくなった場合に目的関数が収束したと判定して処理を終了する収束判定部と、を有することを特徴とする情報処理装置が提供される。

前記演算部において、空間ベースラインa_baseline(x,y)に加えて時間ベースラインv_baseline(t)を導入したモデルに基づいた推定値更新を行うことを特徴とする。

また、前記パラメータ設定部において、スパイク時系列に対するカルシウム濃度変化のインパルス応答に基づいた変換行列により、スパイク時系列とカルシウム濃度変化との関係式を設定し、前記演算部において、変換行列に基づいた推定値更新を行うことを特徴とする。

また、前記パラメータ設定部において、細胞cの細胞形状a_c(x,y)とスパイク時系列u_c(t)の事前分布を細胞の平均サイズ、スパイクの発火率、それらの積である信号強度を考慮した指数分布として設定し、前記演算部において、設定された事前分布に基づいた推定値更新を行うことを特徴とする。

また、前記演算部においてスパイク時系列u_c(t)に加えて細胞形状a_c(x,y)に対する非負拘束条件を導入したモデルに基づいた推定値更新を行うことを特徴とする。

これにより、疎な解を得られるような条件で推定を行うことで関連度自動決定の効果が働き、不要な細胞候補は自動的に縮退し（信号が小さくなり）細胞数が自動的に決定される。

また、前記初期設定部において、想定される細胞のサイズより大きめの特定の領域を覆うような細胞候補を十分な重なりを持たせて画像全体を埋め尽くすように真の細胞の位置に無関係に多数配置することを特徴とする。

さらに、前記スパイク時系列／時間方向のベースライン演算部と、前記細胞形状／空間方向のベースライン演算部と、の出力の少なくともいずれか一方に対して縮退した細胞候補の除去を行う縮退細胞候補除去部を有することを特徴とする。

さらに、前記スパイク時系列／時間方向のベースライン演算部と、前記細胞形状／空間方向のベースライン演算部と、の出力の少なくともいずれか一方に対して類似した細胞候補を統合する類似細胞統合部を有することを特徴とする。

また、前記演算部は、定められたモデル、拘束条件、事前分布のもとでMAP推定の目的関数を減少させる演算を行う。尚、MAP推定はベイズ推定という枠組みの中の比較的簡易な推定方法の例示であり、広く「ベイズ理論に基づく推定」を用いることができる。具体的には、定められたモデル、拘束条件、事前分布のもとで事後分布を近似したテスト関数が真の事後分布に近づくような演算を行う。

このベイズ推定法による定量的カルシウム動態推定については、以下発表を参照することができる。
角田敬正・織田善晃・大森敏明・井上雅司・宮川博義・岡田真人・青西亨、信学技報, vol. 111, no. 483, NC2011-175, pp. 317-322, 2012年3月.
http://www.ieice.org/ken/paper/20120316a0pt/

さらに、前記演算部は、細胞形状a_c(x,y)と空間ベースラインa_baseline(x,y)を固定し、スパイク時系列の変数値が非負という拘束条件のもとで、MAP推定の目的関数を最小化するカルシウム濃度変化v_c(t)と時間ベースラインv_baseline(t)を求める(C+1)T次元の二次計画問題を解くカルシウム濃度変化／時間方向ベースライン演算部と、カルシウム濃度変化v_c(t)と時間ベースラインv_baseline(t)を固定し、細胞形状の変数値が非負という拘束条件のもとで、MAP推定の目的関数を最小化する細胞形状a_c(x,y)と空間ベースラインa_baseline(x,y)を求めるXY個の(C+1)次元の二次計画問題を解く細胞形状／空間方向ベースライン演算部とを有することを特徴とする。

前記カルシウム濃度変化／時間方向ベースライン演算部において、当該二次計画問題を主双対内点法と（前処理付き）共役勾配法を用いて反復的に解くことを特徴とする。

本発明の他の観点によれば、細胞の活動を非負の値で表現される細胞形状a_c(x,y)と非負の値で表現されるスパイク時系列u_c(t)より導出されるカルシウム濃度変化v_c(t)の積により表現し、観測データf(t,x,y)を複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データとモデルと事前分布を元にしたベイズ理論に基づく推定により想定される細胞の位置と活動を求める情報処理装置を用いた情報処理方法であって、細胞のスパイク生成に起因する細胞内のカルシウム濃度の上昇に伴って起こる蛍光タンパク質からの発光位置(x,y)における時刻tでの蛍光強度f(t,x,y)を読み取るデータ読み取りステップと、スパイク時系列u_c(t)とカルシウム濃度変化v_c(t)との関係式と細胞形状a_c(x,y)とスパイク時系列u_c(t)に関する事前分布を設定するパラメータ設定ステップと、前記データ読み取り部により読みとられた位置(x,y)における時刻tでの蛍光強度f(t,x,y)に基づいて細胞形状a_c(x,y)の細胞候補を配置し、ベースラインを蛍光強度f(t,x,y)の統計量に基づいて設定する初期設定ステップと、前記ベイズ理論に基づく推定で導出される最小化問題の目的関数を減少させ、細胞形状a_c(x,y)、カルシウム濃度変化v_c(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は相対的な値であるため、何らかのベースラインを規定する必要がある。解析モデル記憶部２１に記憶される変数を更新する本実施の形態によるモデルは、既存手法のモデルに対して以下の変更を行ったモデルである。

１）時間方向のベースラインv_baseline(t)を導入
既存手法のモデルでは、式（１）のように空間方向のみについてベースラインを考慮したが、本実施の形態では、時系列推定時に時間ベースラインv_baseline(t)を同時に推定することで、解が正しく収束するようにする。

２）
既存手法では、u_cとv_cの対応を式（2）の形に限定していたが、本実施の形態では、より一般的な形に拡張することで、現実のカルシウム応答に近づけより精度の高い解が得られるようにする。

拡張したモデルでは、カルシウム濃度v_c=[v_c(1) v_c(2) ・・・ v_c(T)]^Tはインパルス応答 [g₁,g₂,・・・]とスパイク時系列u_c=[u_c(1) u_c(2) ・・・ u_c(T)]^Tとのたたみ込みで記述される。

たとえば、既存の方法ではインパルス応答を指数関数で{1,γ,γ²,・・・,γ^T}と定義していた。このとき、R=1である。

現実的な拡張としては、カルシウム濃度の立ち上がりにデルタ関数ではなく有限の値の時定数γ₁での立ち上がりを仮定した{1,(γ₁ ²-γ₂ ²)/(γ₁-γ₂),(γ₁ ³-γ₂ ³)/(γ₁-γ₂),・・・,(γ₁ ^T-γ₂ ^T)/(γ₁-γ₂)}という形のインパルス応答を定義して以下の式(14)を用いることが考えられる。このとき、R=2 である。

このように、本実施の形態によれば、わずかな計算量の増加でモデルを拡張することができる。

３）a_cに対する非負拘束条件を導入
上記既存手法では、式(3)および(5)により結果的にペナルティ項が定数となってしまい、解が一意に定まらない。スパイク時系列u_c(t)には非負条件が課されているが、細胞の形状a_c(x,y)の拘束条件が規定されておらず負の値もとり得るため、細胞とは実際には無関係な領域の影響を受けてしまう。本実施の形態では、本来の物理的な条件を反映させ細胞形状a_c(x,y)にも非負の拘束条件を導入する。対応する細胞に対する相関が十分強い領域のみが値を持ち、それ以外がすべて0となる（解が疎になる）ため、細胞の形状が明確になる（図３参照）。

既存手法では、式(3)の事前分布のハイパーパラメータが式(5)により変数自身により決定される自己一貫性が要求されており、解が一意に定まらなかった。また、スパイク時系列のみに事前分布が設定されているため、信号を形成する細胞形状とスパイク時系列の比も一意に定まらなかった。本実施の形態では、細胞形状とスパイク時系列に関して式(16)の形式の事前分布を導入する。

であり、a₀,u₀,s₀は事前分布のハイパーパラメータである。この事前分布はモデルにおいてスパイク時系列の発火率および細胞形状の領域サイズさらにそれらの積で表される細胞の信号強度に対してハイパーパラメータで表現される事前知識を用いることに対応し、MAP推定においてはa_c,u_c,a_cu_cに関するL1ノルム正規化項を加えることに対応する。この事前分布の導入により、MAP推定における解が一意に定まり、式(12)(15)と併せてより疎な解を得ることができる。

目的関数とパラメータ設定処理においては、拘束条件（式(12),(15)）、尤度関数（式(10)）、事前分布（式(16)(17)）のMAP推定は、拘束条件（式(12),(15)）の元で目的関数式(18)を最小化する細胞数Cと変数a_c,a_baseline,v_c,v_baselineを求める問題となる。

σ²は真の解が得られたときの二乗誤差であり現時点では不明であるため、全体を定数σ²倍した上で、パラメータρ_a,ρ_u,ρを用いて次のように置き換える。

ただし、任意の正の定数sに対して、

としても目的関数の値も個々の細胞の信号も変わらずa_c,u_cの比を変化させるだけであるため、ρ_a,ρ_uの比は解の実質的な形に影響しない。このため、パラメータγ,ρ_a,ρ_u,ρを定めて式(19)を最小化する解を求めておいてから、必要であれば正の定数sを用いて(20)の変換を行うことにより最終的なa_c,v_cのスケールを決めればよい。(19)を最小化する解は、適当なa_c,a_baselineを初期条件として設定した上で、a_c,a_baselineを固定することにより(19)をv_c,v_baselineに関する二次計画問題として解く時間変数更新ステップと、v_c,v_baselineを固定することにより式(19)をa_c,a_baselineに関する二次計画問題として解く空間変数更新ステップとを交互に繰り返すことにより、目的関数を減少させ最終的に収束させることで求める。ベースラインの初期条件としては、最小二乗法により簡単に求めることのできる細胞が存在しない場合の最適値を設定し、細胞形状に関しては想定される細胞のサイズより大きめの特定の領域を覆うような細胞候補を十分な重なりを持たせて画像全体を埋め尽くすように多数配置する。不要な細胞候補は関連度自動決定により縮退するため、自動的に正しい細胞数が得られる。

パラメータは、予想される細胞の画像上での大きさA_size、発火率U_freqおよびSN比Sに基づいて設定する。具体的にはa_c,u_cの解がそれぞれ２つの値{0,σ},{0,S}しかとらず、σの値をとるピクセル数がA_size、Sの値をとる確率がU_freqであることを想定して、ρ_a≒σA_size,ρ_u≒STU_freqと設定する。

上記のアルゴリズムを用いることで、ここで設定した問題の大域解を高精度に求めることができる。

次いで、ステップＳ３において、初期設定部３が、ベースライン、細胞の個数および細胞の形状を初期設定する。例えば、細胞が存在しないと仮定した場合の最適なベースラインを初期条件として設定する。

ここで、想定される細胞のサイズより大きめの特定の領域を覆うような細胞候補を十分な重なりを持たせて画像全体を埋め尽くすように真の細胞の位置に無関係に多数配置する。疎な解を得られるような条件でMAP推定を行うことで関連度自動決定の効果が働き、不要な細胞候補は自動的に縮退し（信号が小さくなり）細胞数が自動的に決定される。

次いで、ステップＳ４において、スパイク時系列／時間方向のベースライン演算部５が、スパイク時系列および時間方向のベースラインを求める。

スパイク時系列／時間方向のベースライン演算部５は、拘束条件（式(2)）のもとで、下記の目的関数（式（１７））を最小化するv_c,v_baseline,u_cを求めるCT+T次元の二次計画問題を解く。

a_c,a_baselineを固定すると、拘束条件(12)(15)のもとで目的関数(21)を最小化するv_c,v_baselineを求めるCT+T次元の二次計画問題となる。この問題を主双対内点法と前処理付き共役勾配法を組み合わせて解く。

ただし、e_tはt番目の要素が1となる規定ベクトルとし、パラメータは次の式で表される。

・主双対内点法
二次計画問題(21)の基本的な解法として主双対内点法を導入する。主双対内点法では変数v_c,v_baseline,u_cおよびスラック変数s_cに対して初期条件を設定し、それらの値を目的の二次計画問題の解に近づくように更新していく。更新する方向は次の線形方程式の解により決定される。ここで、βは制御パラメータ、eは要素がすべて１のベクトルである。

このとき、(21)からΔv_baseline,Δs,Δuを消去すると、CT次元の線形方程式(28)が得られる。

ただし、

である。

CT次元の線形方程式は通常の解法では膨大なメモリと計算量が必要となるが、共役勾配法では対象となる行列と任意のベクトルの積が計算できればよいので、(28)のA’およびQ_ijのみ記憶し、その構造を利用して消費メモリO(C²+CTR)および計算量O(C²T+CTR)で行列とベクトルの積を計算することにより、（前処理付きの）共役勾配法を適用する。通常CT×CT次元の行列とCT次元のベクトルの積の計算量および消費メモリはO(C²T²)である。

一方、ステップＳ９において、細胞形状／空間ベースライン演算部５が、細胞形状および空間方向のベースラインを求める。

v_c,v_baselineを固定すると、(15)の拘束条件のもとで目的関数(19)を最小化するa_c,a_baselineを求めるXY個の(C+1)次元の二次計画問題となる。この問題を主双対内点法を用いて解く。ここでは、時間方向と違い問題の次元は比較的小さいので線形方程式の解法にはコレスキー分解を用いた直接法を用いる。このとき、各問題を並列処理で解くことで速度を向上させることができる。

次いで、ステップＳ４、ステップＳ１０の後に、ステップＳ５又はステップＳ１０において、縮退細胞候補除去部１１が、縮退した細胞候補の除去を行う。時系列および細胞形状を求める各ステップについてステップ後に細胞候補の除去を試みる。この際のアルゴリズムとして、縮退した細胞候補を含んだ形で計算を進めても良いが、縮退した細胞候補に関しても通常の細胞候補と同様の計算量がかかるため、適切に早めに除去していくことが望ましい。このため、maxa_c×maxv_c<10^-4σの場合に細胞候補を除去すると良い。

また必要に応じて、スパイク時系列／時間方向のベースライン演算部５と、細胞形状／空間方向のベースライン演算部７と、の出力の少なくともいずれか一方に対して、類似した細胞候補を統合する類似細胞統合部１３を任意に設けても良い（ステップＳ６，ステップＳ１１）。

類似細胞統合部１３は、複数の細胞候補の時系列の比較を行い十分に類似性が高いと判定した場合に、細胞形状を足し合わせた細胞候補を新たに導入し、元の細胞候補を除去することにより、同一細胞が複数の細胞候補に分割されることを防ぐことができる。

ただし、関連度自動決定の効果によりこのような統合処理を必要としない場合もあるため、関連度自動決定の効果に依存して、これを設けるかどうかを決めれば良い。

次いで、ステップＳ７又はステップＳ１２において、反復・収束判定部１５が、収束判定を行う。目的関数の変化が十分小さくかつ除去された細胞候補がない場合に収束したと判定して終了する（ステップＳ８、Ｓ１３）。ステップＳ７で収束していないと判定されると、ステップＳ９に進み、ステップＳ１２で収束していないと判定されると、ステップＳ４に戻る。

尚、収束の判定は、例えば、細胞がない場合の二乗誤差をσ₀ ²、目的関数を０．５σ₀ ²としたとき、各ステップでの目的関数の変化が、閾値０．００１σ₀ ²以下であれば収束したとして終了することができる。この閾値の値は、経験則で決めることができるが、この方法に限定されるものではない。

以上の処理を繰り返すことで収束すると、細胞数Cと変数a_c,a_baseline,v_c,v_baselineおよびu_cすべてが全て求まり、これらの値を出力部１７が結果を出力することができる。

上記の実施の形態において、添付図面に図示されている構成等については、これらに限定されるものではなく、本発明の効果を発揮する範囲内で適宜変更することが可能である。その他、本発明の目的の範囲を逸脱しない限りにおいて適宜変更して実施することが可能である。また、本発明の各構成要素は、任意に取捨選択することができ、取捨選択した構成を具備する発明も本発明に含まれるものである。

また、本実施の形態で説明した機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより各部の処理を行ってもよい。尚、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。

本発明は、情報処理装置として利用可能である。

Ａ…情報処理装置、１…データ読み込み部、３…初期設定部、５…スパイク時系列／時間ベースライン演算部、７…細胞形状／空間ベースライン演算部、１１…縮退細胞除去部、１５…反復・収束判定部、１７…出力部、２１…解析モデル記憶部。

Claims

細胞の活動を非負の値で表現される細胞形状a_c(x,y)と非負の値で表現されるスパイク時系列u_c(t)より導出されるカルシウム濃度変化v_c(t)の積により表現し、観測データf(t,x,y)を複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データとモデルと事前分布を元にしたベイズ理論に基づく推定により想定される細胞の位置と活動を求める情報処理装置であって、
細胞のスパイク生成に起因する細胞内のカルシウム濃度の上昇に伴って起こる蛍光タンパク質からの発光を位置(x,y)における時刻tでの蛍光強度f(t,x,y)として読み取るデータ読み取り部と、
スパイク時系列u_c(t)とカルシウム濃度変化v_c(t)との関係式と細胞形状a_c(x,y)とスパイク時系列u_c(t)に関する事前分布を設定するパラメータ設定部と、
前記データ読み取り部により読みとられた位置(x,y)における時刻tでの蛍光強度f(t,x,y)に基づいて細胞形状a_c(x,y)の細胞候補を配置し、ベースラインを蛍光強度f(t,x,y)の統計量に基づいて設定する初期設定部と、
前記ベイズ理論に基づく推定で導出される最小化問題の目的関数を減少させ、細胞形状a_c(x,y)、カルシウム濃度変化v_c(t)およびベースラインをより真の解に近い推定値に更新する演算部と、
前記演算部における更新で目的関数の変化が十分小さくなった場合に目的関数が収束したと判定して処理を終了する収束判定部とを有することを特徴とする情報処理装置。
前記演算部において、
空間ベースラインa_baseline(x,y)に加えて時間ベースラインv_baseline(t)を導入したモデルに基づいた推定値更新を行うことを特徴とする請求項１に記載の情報処理装置。
前記パラメータ設定部において、
スパイク時系列に対するカルシウム濃度変化のインパルス応答に基づいた変換行列により、スパイク時系列とカルシウム濃度変化との関係式を設定し、
前記演算部において、
変換行列に基づいた推定値更新を行うことを特徴とする請求項１又は２に記載の情報処理装置。
前記パラメータ設定部において、

前記演算部において、
設定された事前分布に基づいた推定値更新を行うことを特徴とする請求項１から３までのいずれか１項に記載の情報処理装置。
前記初期設定部において、
想定される細胞のサイズより大きめの特定の領域を覆うような細胞候補を十分な重なりを持たせて画像全体を埋め尽くすように真の細胞の位置に無関係に多数配置することを特徴とする請求項１から４までのいずれか１項に記載の情報処理装置。
細胞の活動を非負の値で表現される細胞形状a_c(x,y)と非負の値で表現されるスパイク時系列u_c(t)より導出されるカルシウム濃度変化v_c(t)の積により表現し、観測データf(t,x,y)を複数の細胞活動の線形和とベースラインを用いて記述するモデルに基づいて、観測データとモデルと事前分布を元にしたベイズ理論に基づく推定により想定される細胞の位置と活動を求める情報処理装置を用いた情報処理方法であって、
細胞のスパイク生成に起因する細胞内のカルシウム濃度の上昇に伴って起こる蛍光タンパク質からの発光位置(x,y)における時刻tでの蛍光強度f(t,x,y)を読み取るデータ読み取りステップと、
スパイク時系列u_c(t)とカルシウム濃度変化v_c(t)との関係式と細胞形状a_c(x,y)とスパイク時系列u_c(t)に関する事前分布を設定するパラメータ設定ステップと、
前記データ読み取りステップにより読みとられた位置(x,y)における時刻tでの蛍光強度f(t,x,y)に基づいて細胞形状a_c(x,y)の細胞候補を配置し、ベースラインを蛍光強度f(t,x,y)の統計量に基づいて設定する初期設定ステップと、
前記ベイズ理論に基づく推定で導出される最小化問題の目的関数を減少させ、細胞形状a_c(x,y)、カルシウム濃度変化v_c(t)およびベースラインをより真の解に近い推定値に更新する演算ステップと、
前記演算ステップにおける更新で目的関数の変化が十分小さくなった場合に目的関数が収束したと判定して処理を終了する収束判定ステップと
を有することを特徴とする情報処理方法。
コンピュータに請求項６に記載の情報処理方法を実行させるためのプログラム。