JP2014095678A - Information processing apparatus, and information processing program - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a technology to estimate, from a state where the number of cells included in an image is unknown, the number of cells, and the shapes and activity time series of the respective cells quickly in a sufficiently accurate and stable manner, without the need to separately set initial conditions.SOLUTION: An information processing apparatus represents the activity of cells by the product of the shapes of the cells a(x,y) represented by a non-negative value and a change in calcium concentration v(t) derived from the spike time series u(t) represented by a non-negative value; and determines, on the basis of a model in which observation data f(t,x,y) is described by using the linear sum of the plurality of cell activities and a baseline, the positions and activity of the cells assumed from the estimation based on the Bayesian theory obtained from the observation data, the model, and a prior distribution.

Description

  The present invention relates to information processing technology, and more particularly to information processing technology for estimating action potentials of nerve cells by image processing or the like.

  Neurons transmit information to other neurons when action potentials (spikes) are generated. In order to understand how neuronal populations in the brain perform information processing, it is necessary to simultaneously record the action potentials of multiple neurons. When an action potential is generated, the intracellular calcium concentration rises. Therefore, by recording a brain with a fluorescent protein that shines depending on the calcium concentration as a moving image with a two-photon microscope, the action potential is applied to multiple cells. You can observe what is happening.

  However, in order to actually handle the data as a neuron population, it is necessary to extract the time series of the position and activity of each cell from the moving image data. Actually, since the data contains a lot of noise, what is obtained is not fixed information but information estimated based on the data. It is important to estimate cell positions and activity time series from high-capacity recorded data at high speed and high accuracy.

  The method that is mainly used as a conventional method is to first find a region that shows a similar activity pattern using principal component analysis, etc., and then perform some processing to determine the individual cells called ROI (Region of Interest). The procedure is such that the shape is fixed, the average of the fluorescent protein reaction in the region is regarded as the activity of the cell, and the spike timing is accurately obtained if necessary (see Non-Patent Document 1 below). This non-patent document relates to an ROI determination algorithm based on principal component analysis / independent component analysis.

  Actually, the cell population has a three-dimensional structure, and there are many areas where multiple cells overlap in the image data projected in two dimensions. As a paper on the model that forms the basis of the invention, in recent years, the idea of using the sum of the activities of multiple cells as an observation image model has been proposed. The demonstration shows that the activity can be estimated (see Non-Patent Document 2 below).

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).

  The problem with the method as described in Non-Patent Document 1 is that the cell population actually has a three-dimensional structure, and there are many regions where a plurality of cells overlap in image data projected in two dimensions. Is not considering.

  Further, the method of Non-Patent Document 2 has many problems in processing actual large-capacity data.

Below, the content of the existing process currently disclosed by the nonpatent literature 2 and its problem are demonstrated, referring drawings.
(Existing processing)
FIG. 4 is a functional block diagram showing a configuration example of an information processing apparatus that performs the existing processing shown in Non-Patent Document 2, and FIG. 5 is a flowchart showing the flow of the existing processing shown in Non-Patent Document 2. is there. First, when information processing is started in step S101, the data reading unit 101 reads data in step S102. As data to be read, fluorescence intensity f (t, x, y) at time t at position (x, y) is given as data. However, it is assumed that t, x, and y are natural numbers less than or equal to T, X, and Y, respectively, and have a scale determined from the sampling rate and spatial resolution. Since the fluorescence intensity f is a relative value, it is necessary to define some baseline.

The calcium concentration with respect to the time axis increases in a delta function by a spike (see FIG. 3), and then decays. From u _c (0) = v _c (0), u _c and v _c can be derived from one time series by linear transformation.

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

Baseline is initially set as average or median.

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

v _c (t) = 0 and σ ² is set from the equation (4).

However, how to find the number of cells C and the shape of cell a _c (x, y) is not specified, and some procedure to obtain these initial conditions must be set separately, but the solution converges appropriately In order to do this, the initial conditions need to be close enough to the solution.

Next, in step S104, the calcium concentration / spike time series calculation unit 105 obtains the calcium concentration v and the spike time series u. In the above formulas (1), (2), and (3), when a _c and a _baseline and parameters σ and λ _c are fixed and MAP (posterior distribution maximization) estimation for u and v is performed, formula ( The following quadratic programming problem is solved under the condition of 2). At this time, since the baseline is fixed and the non-negative constraint condition is set, the influence of the fixed baseline is strongly received, and an appropriate solution as the entire system cannot be obtained.

  This minimization problem is solved using Log barrier method and Newton method. If C = 1, it becomes a tridiagonal matrix and can be solved relatively easily, but if C> 1, it is necessary to consider a triblock diagonal matrix and it is difficult to solve with high accuracy.

  Next, in step S105, the normalization unit 107 performs a process of normalizing the spike time series and cell shape of each cell candidate.

としても信号全体は変わらない。 The signal of each cell is expressed as f _c (t, x, y) = a _c (x, y) v _c (t)

But the whole signal doesn't change.

However, Equation (7) because it has a sum of the penalty term for least squares term and u, the equation (2), the value of the objective function penalty term decreases larger the s _c is smaller End up. This has no meaning for the penalty term, so it must be normalized in some way. Therefore, s _c = maxv _c . However, since the problem is not solved by entering the constraint condition s _c = maxv _c , the solution of the equation (7) is not obtained after normalization.

In step S106, the cell shape / baseline calculation unit obtains the cell shape _ac and the baseline a _baseline . Here, u and v are fixed, and a _c and a _baseline that minimize Equation (7) are obtained. When the maximum likelihood method is performed, each pixel can be considered independently, and Equation (9) may be solved.

This is easy to solve because there is no constraint, but since a _c (x, y) includes negative values, the shape of the cell is strongly affected by noise and is not stable.

  In step S107, the parameter update unit 111 updates the parameters. The parameters are updated according to equations (4) and (5). However, since the parameter to be specified by the true solution is approximated with the current temporary value and the prior distribution is defined from the variable, the objective function cannot be set correctly.

  The iterative / convergence determining unit 117 performs the iterative / convergent determination for the processes in steps S104 and S107. When the process ends in step S108, the output unit 121 outputs the result.

  As described in the description of the normalization process, since the problem setting is poor in the first place and the final convergence determination is impossible, each step is terminated after being repeated an appropriate number of times.

  Thus, the problems of Non-Patent Document 2 are that the number of cells is fixed in the algorithm, that the cell shape must be sufficiently close to the correct shape as an initial condition, the speed and For example, the accuracy is insufficient and proper convergence cannot be obtained when there are a large number of cells.

  The present invention is a technique for estimating the number of cells, the shape of each cell, and the time series of activity stably at a high speed with sufficient accuracy without individually setting initial conditions from a state in which the number of cells contained in the image is unknown. The purpose is to provide.

  In the present invention, the accuracy of a solution obtained by introducing a normalization term (penalty term) by setting a non-negative constraint condition of a parameter and an appropriate prior distribution into a model is improved, and a large number of cell candidates are overlapped in an image. We have developed a device (program) that can correctly estimate the number of cells and parameters by setting initial conditions that are filled in tightly and by devising selection of cell candidates. The main algorithm is to use an iterative method that iteratively and repeatedly improves shape and time change, removes cell candidates that have only small signals below the noise level, and integrates cell candidates that have similar activity states. Detects all active cells in the image correctly while reducing the amount of computation. To solve the nonlinear programming required at each step, a unique parallel implementation using the main dual interior point method and the conjugate gradient method as a fast solution method for large-scale linear equations is used.

  Finally, it is possible to stably estimate the number of cells, the shape of each cell, and the time series of each cell quickly and with sufficient accuracy without individually setting the initial conditions from the state in which the number of cells contained in the image is unknown. A program that can be obtained.

  Conventional method 2 has a problem that the number of cells is fixed and initial conditions are severe. Therefore, an initial condition is set so that a large number of cell candidates are filled with an overlap in the image, and a condition that makes the solution sparse is added to the model. Due to the effect of the former device, a true cell is sufficiently close to one of the cell candidates so that it can be detected without being missed. Even if there are a large number of cell candidates, it is possible to prevent a true cell from being divided into a plurality of cell candidates more than necessary. In fact, repeated iterations of the algorithm will result in cell candidates with a large signal due to an effect called automatic relevance determination and cell candidates with only a small signal below the noise level. By integrating cell candidates having similar activity states, it is possible to correctly detect all activated cells finally included in the image. The newly proposed method supports the convergence of the iterative method correctly by putting an appropriate constraint on the model that the observed image is the sum of the activities of each cell.

  This method requires a lot of memory and computational complexity due to the nature of preparing a large number of cell candidates, but it is a policy to divide the image into appropriate sizes, pre-process them, and integrate the results. By doing so, the amount of calculation can be greatly reduced, and the objective can be achieved in a realistic execution time by an appropriate solution in each step and a high-speed implementation by parallelization.

  Finally, it is possible to stably estimate the number of cells, the shape of each cell, and the time series of each cell quickly and with sufficient accuracy without individually setting the initial conditions from the state in which the number of cells contained in the image is unknown. I can do it now.

According to one aspect of the present invention, calcium derived from a cell shape a _c (x, y) expressed by a non-negative value of a cell activity and a spike time series u _c (t) expressed by a non-negative value Based on a model that expresses the observed data f (t, x, y) using a linear sum of multiple cell activities and a baseline, expressed as a product of the concentration change v _c (t). An information processing device that obtains the location and activity of cells that are assumed by estimation based on Bayesian theory based on the distribution, from fluorescent proteins caused by the increase in intracellular calcium concentration caused by cell spike generation Data reading unit that reads luminescence as fluorescence intensity f (t, x, y) at time t at position (x, y), and relationship between spike time series u _c (t) and calcium concentration change v _c (t) wherein the cell shape a _c (x, y) and parameters to configure a pre-distribution for spike time series u _c (t) A setting unit, the data reading unit by read position (x, y) fluorescence intensity f at time t in the (t, x, y) the cell candidate cell shape a _c (x, y) on the basis of the The initial setting unit that sets and sets the baseline based on the statistic of the fluorescence intensity f (t, x, y), and reduces the objective function of the minimization problem derived by the estimation based on the Bayesian theory, The calculation unit that updates the shape a _c (x, y), the calcium concentration change v _c (t) and the baseline to an estimated value closer to the true solution, and the update in the calculation unit sufficiently reduces the change in the objective function. And a convergence determining unit that determines that the objective function has converged and terminates the processing when the objective function is converged.

The calculation unit updates the estimated value based on a model in which a time baseline v _baseline (t) is introduced in addition to a spatial baseline a _baseline (x, y).

  Further, in the parameter setting unit, a relational expression between the spike time series and the calcium concentration change is set by a conversion matrix based on the impulse response of the calcium concentration change with respect to the spike time series, and based on the conversion matrix in the calculation unit The estimated value is updated.

Further, in the parameter setting unit, the prior distribution of the cell shape a _c (x, y) and the spike time series u _c (t) of the cell _{c represents} the average cell size, the firing rate of the spike, and the signal intensity that is the product of them. Is set as an exponential distribution in consideration of the above, and the calculation unit updates the estimated value based on the set prior distribution.

Further, the calculation unit updates the estimated value based on a model in which a non-negative constraint condition for the cell shape a _c (x, y) is introduced in addition to the spike time series u _c (t).

  In this way, the effect of automatic relevance determination works by performing estimation under conditions that can obtain a sparse solution, and unnecessary cell candidates are automatically degenerated (signal becomes smaller) and the number of cells is automatically determined. Is done.

  In the initial setting unit, a large number of cell candidates covering a specific region larger than the assumed cell size are sufficiently overlapped regardless of the true cell position so as to fill the entire image. It is characterized by arranging.

  Further, a degenerate cell candidate that removes a degenerated cell candidate with respect to at least one of the output of the spike time series / time direction baseline calculation unit and the cell shape / space direction baseline calculation unit It has a removal part.

  Further, a similar cell integration unit that integrates similar cell candidates with respect to at least one of the outputs of the spike time series / time direction baseline calculation unit and the cell shape / space direction baseline calculation unit It is characterized by having.

  In addition, the calculation unit performs a calculation for reducing the objective function of the MAP estimation based on a predetermined model, constraint conditions, and prior distribution. MAP estimation is an example of a relatively simple estimation method in the framework of Bayesian estimation, and “estimation based on Bayesian theory” can be widely used. Specifically, an operation is performed such that a test function approximating the posterior distribution under a predetermined model, constraint condition, and prior distribution approaches the true posterior distribution.

The following announcement can be referred to for quantitative calcium dynamics estimation by this Bayesian estimation method.
Takamasa Tsunoda, Yoshiaki Oda, Toshiaki Omori, Masashi Inoue, Hiroyoshi Miyagawa, Masato Okada, Jun Aonishi, IEICE Tech., Vol. 111, no. .
http://www.ieice.org/ken/paper/20120316a0pt/

Furthermore, the calculation unit fixes the cell shape a _c (x, y) and the spatial baseline a _baseline (x, y), and performs the MAP estimation under the constraint that the spike time series variable value is non-negative. Calculating calcium concentration change v _c (t) and time baseline v _baseline (t) to minimize objective function (C + 1) Calculating calcium concentration change / time direction baseline calculation unit to solve T-dimensional quadratic programming problem Cell shape a _c that minimizes the objective function of MAP estimation under the constraint that the calcium concentration change v _c (t) and the time baseline v _baseline (t) are fixed and the cell shape variable value is non-negative It has a cell shape / space direction baseline calculation unit that solves XY (C + 1) -dimensional quadratic programming problems to obtain (x, y) and a space baseline a _baseline (x, y) To do.

  In the calcium concentration change / time direction baseline calculation section, the quadratic programming problem is iteratively solved using a main dual interior point method and a conjugate gradient method (with preprocessing).

According to another aspect of the present invention, cell activity is derived from a cell shape a _c (x, y) expressed by a non-negative value and a spike time series u _c (t) expressed by a non-negative value. Based on a model that expresses the observation data f (t, x, y) using a linear sum of multiple cell activities and a baseline, expressed as a product of the calcium concentration change v _c (t). An information processing method using an information processing device that calculates the location and activity of a cell that is assumed by estimation based on Bayesian theory based on prior distribution, in which the intracellular calcium concentration is increased due to cell spike generation. A data reading step for reading the fluorescence intensity f (t, x, y) at time t at the emission position (x, y) from the fluorescent protein that accompanies, spike time series u _c (t) and calcium concentration change v _c (t) and the relationship cell shape a _c (x, y) and pre relates spike time series u _c (t) A parameter setting step of setting the fabric, the data reading unit by read position (x, y) fluorescence intensity f at time t in the (t, x, y) cell shape based on a _c (x, y ) Cell candidates, and an initial setting step for setting a baseline based on a statistic of fluorescence intensity f (t, x, y), and an objective function for a minimization problem derived by estimation based on the Bayesian theory And _calculating the cell shape a _c (x, y), the calcium concentration change v _c (t) and the baseline to an estimated value closer to the true solution, and updating the objective function by updating in the calculation step. There is provided an information processing method comprising a convergence determination step for determining that the objective function has converged when the change becomes sufficiently small and terminating the processing.

  The present invention may be a program for causing a computer to execute the information processing method described above, or a computer-readable storage medium for recording the program.

  According to the present invention, the number of cells, the shape of each cell, and the time series of activity are stably estimated with sufficient accuracy at high speed without individually setting initial conditions from a state in which the number of cells contained in the image is unknown. be able to.

It is a functional block diagram which shows the example of 1 structure of the information processing apparatus by one embodiment of this invention. It is a flowchart figure which shows the flow of the information processing by the information processing apparatus by one embodiment of this invention. It is a flowchart figure which shows the flow of the information processing by the information processing apparatus by one embodiment of this invention, and is a figure following FIG. 2A. It is a schematic diagram showing an outline of information processing by the information processing apparatus according to an embodiment of the present invention. FIG. 11 is a functional block diagram illustrating a configuration example of an information processing device described in Non-Patent Document 2. FIG. 10 is a flowchart showing a flow of information processing by the information processing apparatus described in Non-Patent Document 2.

Hereinafter, an information processing technique according to an embodiment of the present invention will be described with reference to the drawings.
The present invention improves the accuracy of a solution obtained by introducing a normalization term (penalty term) into a model by setting a non-negative constraint condition of parameters and an appropriate prior distribution, and then overlapping a large number of cell candidates in an image. The present invention relates to a technology that can correctly estimate the number of cells and parameters by setting initial conditions so as to fill them tightly and devising selection of cell candidates. The main algorithm is to use an iterative method that iteratively and repeatedly improves shape and time change, removes cell candidates that have only small signals below the noise level, and integrates cell candidates that have similar activity states. Detects all active cells in the image correctly while reducing the amount of computation. To solve the nonlinear programming required at each step, a unique parallel implementation using the main dual interior point method and the conjugate gradient method as a fast solution method for large-scale linear equations is used.

  FIG. 3 is a schematic diagram showing an outline of information processing according to the present embodiment. A moving image in which a moving image consisting of a large number of images is observed, and light emission from a fluorescent protein caused by an increase in intracellular calcium concentration when an action potential is generated at a certain time t and at that time t is observed with a two-photon microscope Based on the above, the fluorescence intensity f (t, x, y) can be obtained.

  By reading the observed value of this fluorescence intensity f (t, x, y), setting the initial conditions, and calculating the sum of the product of the cell shape and the corresponding calcium concentration of the spike train, the number of cells and each Cell shape and activity time series.

In the present embodiment, in order to solve the problems of the existing method, the following device has been devised.
1) Add a condition to the model of the existing method to clarify the problem setting so that the solution can be obtained uniquely.
2) At the same time, improve the model so that it has the solution properties that meet the objectives.
3) Assuming a state in which the number of cells and the cell position are undecided, an appropriate cell number and cell position can be correctly estimated from a large number of cell candidates by automatically determining the degree of association.
4) Show a solution based on a model that enables high-precision estimation even for large-scale data.

  1 is a functional block diagram showing an example of the configuration of the information processing apparatus (image processing apparatus) according to the present embodiment, and FIGS. 2A and 2B are functional block diagrams showing the flow of information processing. These are the typical figures which show the outline | summary of information processing.

  First, when processing is started (step S1), in step S2, the data reading unit 1 starts from a fluorescent protein caused by an increase in intracellular calcium concentration when generating an action potential at a certain time t and that time t. Loads moving image data observed with a two-photon microscope and sets parameters. As data to be read, fluorescence intensity f (t, x, y) at time t at position (x, y) is given as data. However, it is assumed that t, x, and y are natural numbers less than or equal to T, X, and Y, respectively, and have a scale determined from the sampling rate and spatial resolution. If the sampling rate is S [Hz] and the recording time is R [sec], T = SR. X and Y are the number of pixels of the image. Since the fluorescence intensity f is a relative value, it is necessary to define some baseline. The model according to the present embodiment that updates the variables stored in the analysis model storage unit 21 is a model obtained by making the following changes to the model of the existing method.

1) Introduction of _{baseline in} time direction v _baseline (t) In the model of the existing method, the baseline is considered only in the spatial direction as shown in Equation (1), but in this embodiment, the time base is used for time series estimation. By simultaneously estimating the line v _baseline (t), the solution converges correctly.

2)
In the existing method, the correspondence between u _c and v _c is limited to the form of Equation (2). However, in this embodiment, the accuracy is closer to the actual calcium response by extending to a more general form. To obtain a high solution.

In the expanded model, the calcium concentration v _c = [v _c (1) v _c (2) ・ ・ ・ v _c (T)] ^T is the impulse response [g ₁ , g ₂ , ...] and the spike time series u _c = [u _c (1) u _c (2) ... u _c (T)] Described by convolution with ^T.

For example, in the existing method, the impulse response is defined as an exponential function {1, γ, γ ² ,..., Γ ^T }. At this time, R = 1.

As a practical extension, we assumed that the calcium concentration rises with a finite time constant γ ₁ instead of a delta function {1, (γ ₁ ² -γ ₂ ² ) / (γ ₁ -γ ₂ ) , (γ ₁ ³ -γ ₂ ³ ) / (γ ₁ -γ ₂ ), ..., (γ ₁ ^T -γ ₂ ^T ) / (γ ₁ -γ ₂ )} It is conceivable to use the following formula (14). At this time, R = 2.

  Thus, according to the present embodiment, the model can be expanded with a slight increase in calculation amount.

3) Introducing a non-negative constraint condition for a _{c In the} above existing method, the penalty term becomes a constant as a result of equations (3) and (5), and the solution is not uniquely determined. The spike time series u _c (t) has a non-negative condition, but the cell shape a _c (x, y) is not defined and can be negative. Will be affected by unrelated areas. In the present embodiment, non-negative constraint conditions are also introduced to the cell shape a _c (x, y) to reflect the original physical conditions. Only a region having a sufficiently strong correlation with the corresponding cell has a value, and all others are 0 (the solution becomes sparse), so that the shape of the cell becomes clear (see FIG. 3).

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

The existing method requires self-consistency in which the hyperparameter of the prior distribution of Equation (3) is determined by the variable itself according to Equation (5), and the solution cannot be determined uniquely. In addition, since the prior distribution is set only for the spike time series, the ratio of the cell shape forming the signal to the spike time series could not be determined uniquely. In the present embodiment, a prior distribution in the form of Equation (16) is introduced with respect to the cell shape and spike time series.

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

A ₀ , u ₀ , s ₀ are hyperparameters of the prior distribution. This prior distribution corresponds to the use of prior knowledge expressed by hyperparameters in the model for the firing rate of spike time series and the region size of the cell shape and the signal strength of the cell represented by the product of them. Corresponds to adding an L1 norm normalization term for a _c , u _c , a _c u _c . By introducing this prior distribution, the solution in the MAP estimation is uniquely determined, and a sparse solution can be obtained together with the equations (12) and (15).

In objective function and parameter setting processing, MAP estimation of constraint conditions (Equations (12) and (15)), likelihood functions (Equation (10)), and prior distributions (Equations (16) and (17)) This is a problem of obtaining the number of cells C and variables a _c , a _baseline , v _c , and v _baseline that minimize the objective function equation (18) under (Equations (12) and (15)).

Since σ ² is a square error when a true solution is obtained and is unknown at this time, the whole is multiplied by a constant σ ² and replaced by the following using parameters ρ _a , ρ _u , ρ.

としても目的関数の値も個々の細胞の信号も変わらず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に関する二次計画問題として解く空間変数更新ステップとを交互に繰り返すことにより、目的関数を減少させ最終的に収束させることで求める。ベースラインの初期条件としては、最小二乗法により簡単に求めることのできる細胞が存在しない場合の最適値を設定し、細胞形状に関しては想定される細胞のサイズより大きめの特定の領域を覆うような細胞候補を十分な重なりを持たせて画像全体を埋め尽くすように多数配置する。不要な細胞候補は関連度自動決定により縮退するため、自動的に正しい細胞数が得られる。 However, for any positive constant s,

However, the ratio of ρ _a and ρ _u does not affect the substantial form of the solution because the value of the objective function and the signal of the individual cells are not changed and only the ratio of a _c and u _c is changed. For this reason, after determining the parameters γ, ρ _a , ρ _u , ρ and obtaining a solution that minimizes the equation (19), the conversion of (20) is performed using a positive constant s if necessary. Thus, the final scale of a _c and v _c can be determined. The solution to minimize (19) is the quadratic programming problem for v _c , v _baseline by setting a _c , a _baseline as an initial condition and fixing a _c , a _baseline. By alternately repeating the time variable update step solved as follows and the spatial variable update step solved as a quadratic programming problem with respect to a _c , a _baseline by fixing v _c , v _baseline , the objective function is Decrease and finally converge. As an initial condition of the baseline, an optimal value is set when there is no cell that can be easily obtained by the least square method, and the cell shape covers a specific region larger than the assumed cell size. A large number of cell candidates are arranged so as to fill the entire image with sufficient overlap. Since unnecessary cell candidates are degenerated by automatically determining the degree of association, the correct number of cells is automatically obtained.

The parameters are set based on the expected size A _size of the cell, firing rate U _freq, and SN ratio S. Specifically, the solutions of a _c and u _c have only two values {0, σ} and {0, S}, respectively, and the number of pixels taking the value of σ is A _size and the probability of taking the value of S is U Assuming _freq , ρ _a ≈σ A _size and ρ _u ≈STU _freq are set.

  By using the above algorithm, the global solution of the problem set here can be obtained with high accuracy.

  Next, in step S3, the initial setting unit 3 initializes the baseline, the number of cells, and the shape of the cells. For example, an optimal baseline when it is assumed that no cells exist is set as an initial condition.

  Here, a large number of cell candidates covering a specific area larger than the assumed cell size are arranged regardless of the true cell position so as to fill the entire image with sufficient overlap. By performing MAP estimation under conditions that can obtain a sparse solution, the effect of automatic relevance determination works, unnecessary cell candidates are automatically degenerated (signal becomes smaller), and the number of cells is automatically determined. .

  Next, in step S4, the spike time series / time direction baseline calculation unit 5 obtains the spike time series and time direction baseline.

The baseline calculation unit 5 in the spike time series / time direction calculates v _c , v _baseline , u _c that minimizes the following objective function (equation (17)) under the constraint condition (equation (2)). Solve the CT + T-dimensional quadratic programming problem.

If a _c and a _baseline are fixed, it becomes a CT + T-dimensional quadratic programming problem for obtaining v _c and v _baseline that minimizes the objective function (21) under the constraint conditions (12) and (15). This problem is solved by combining the main dual interior point method and the preconditioned conjugate gradient method.

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

Here, _et is a defined vector whose t-th element is 1, and the parameters are expressed by the following equations.

・ Main dual interior point method The main dual interior point method is introduced as a basic solution of the quadratic programming problem (21). In the main dual interior point method, initial conditions are set for variables v _c , v _baseline , u _c and slack variable s _c , and these values are updated so as to approach the solution of the target quadratic programming problem. The updating direction is determined by the solution of the following linear equation. Here, β is a control parameter, and e is a vector whose elements are all ones.

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

At this time, if Δv _baseline , Δs, Δu is deleted from (21), a CT-dimensional linear equation (28) is obtained.

ただし、

である。

However,

It is.

The CT linear equation requires a large amount of memory and computational complexity, but the conjugate gradient method only needs to be able to calculate the product of the target matrix and any vector, so A 'and Q in (28) Conjugate (with preprocessing) by memorizing only _ij and using the structure to calculate the product of matrix and vector with consumption memory O (C ² + CTR) and complexity O (C ² T + CTR) Apply the gradient method. Usually, the amount of computation and memory consumption of the product of a CT × CT dimensional matrix and a CT dimensional vector is O (C ² T ² ).

  On the other hand, in step S9, the cell shape / space baseline calculation unit 5 obtains a cell shape and a baseline in the spatial direction.

When v _c and v _baseline are fixed, XY (C + 1) -dimensional quadratic programming problems for obtaining a _c and a _baseline that minimize the objective function (19) under the constraint of (15) Become. This problem is solved using the principal dual interior point method. Here, since the dimension of the problem is relatively small unlike the time direction, the direct method using Cholesky decomposition is used to solve the linear equation. At this time, the speed can be improved by solving each problem by parallel processing.

Next, after step S4 and step S10, the degenerated cell candidate removal unit 11 removes the degenerated cell candidate in step S5 or step S10. For each step of obtaining the time series and the cell shape, an attempt is made to remove the cell candidate after the step. As an algorithm at this time, the calculation may be performed in a form including degenerated cell candidates, but degenerate cell candidates also require the same amount of calculation as normal cell candidates, so they are removed appropriately and early It is desirable. For this reason, cell candidates are preferably removed when maxa _c × maxv _c <10 ⁻⁴ σ.

  If necessary, similar cell candidates are integrated into at least one of the output of the spike time series / time direction baseline calculation unit 5 and the cell shape / space direction baseline calculation unit 7. The similar cell integration unit 13 may be optionally provided (step S6, step S11).

  When the similar cell integration unit 13 compares the time series of a plurality of cell candidates and determines that the similarity is sufficiently high, the similar cell integration unit 13 newly introduces cell candidates obtained by adding the cell shapes and removes the original cell candidates By doing so, it is possible to prevent the same cell from being divided into a plurality of cell candidates.

  However, since there is a case where such integration processing is not required due to the effect of automatic determination of the degree of association, it may be determined whether to provide this depending on the effect of automatic determination of the degree of association.

  Next, in step S7 or step S12, the iterative / convergence determination unit 15 performs convergence determination. If the change in the objective function is sufficiently small and there is no removed cell candidate, it is determined that the target function has converged, and the process ends (steps S8 and S13). If it is determined in step S7 that it has not converged, the process proceeds to step S9. If it is determined in step S12 that it has not converged, the process returns to step S4.

In the convergence determination, for example, when the square error when there is no cell is σ ₀ ² and the objective function is 0.5σ ₀ ² , the change of the objective function in each step is a threshold value of 0.001σ ₀ ² or less. If so, it can be concluded that it has converged. The threshold value can be determined by an empirical rule, but is not limited to this method.

When convergence is achieved by repeating the above processing, the number of cells C and all of the variables a _c , a _baseline , v _c , v _baseline, and u _c are obtained, and the output unit 17 can output the results of these values.

  In the above-described embodiment, the configuration and the like illustrated in the accompanying drawings are not limited to these, and can be appropriately changed within a range in which the effect of the present invention is exhibited. In addition, various modifications can be made without departing from the scope of the object of the present invention. Each component of the present invention can be arbitrarily selected, and an invention having a selected configuration is also included in the present invention.

  In addition, a program for realizing the functions described in the present embodiment is recorded on a computer-readable recording medium, and the program recorded on the recording medium is read into a computer system and executed to execute processing of each unit. May be performed. The “computer system” here includes an OS and hardware such as peripheral devices.

  The present invention can be used as an information processing apparatus.

A ... information processing device, 1 ... data reading unit, 3 ... initial setting unit, 5 ... spike time series / time baseline calculation unit, 7 ... cell shape / space baseline calculation unit, 11 ... degenerate cell removal unit, 15 ... Repetition / convergence determination unit, 17 ... output unit, 21 ... analysis model storage unit.

Claims (11)

The product of the cell shape a _c (x, y) expressed as a non-negative value of cell activity and the calcium concentration change v _c (t) derived from the spike time series u _c (t) expressed as a non-negative value Based on a model that describes observation data f (t, x, y) using a linear sum of multiple cell activities and a baseline, based on Bayesian theory based on observation data, model, and prior distribution An information processing apparatus for obtaining a cell position and activity assumed by estimation,
Data reading unit that reads light emission from the fluorescent protein accompanying the increase in intracellular calcium concentration caused by cell spike generation as fluorescence intensity f (t, x, y) at time t at position (x, y) When,
A parameter setting unit for setting a prior expression for the relationship between the spike time series u _c (t) and the calcium concentration change v _c (t), the cell shape a _c (x, y), and the spike time series u _c (t); ,
A cell candidate having a cell shape a _c (x, y) is arranged based on the fluorescence intensity f (t, x, y) at the time t at the position (x, y) read by the data reading unit, and the base An initial setting unit for setting a line based on a statistic of fluorescence intensity f (t, x, y);
Reduce the objective function of the minimization problem derived by the estimation based on the Bayesian theory, and estimate the cell shape a _c (x, y), calcium concentration change v _c (t), and baseline closer to the true solution An arithmetic unit to be updated to
An information processing apparatus comprising: a convergence determining unit that determines that the objective function has converged when the change in the objective function is sufficiently small by the update in the arithmetic unit and ends the process. In the calculation unit,
The information processing apparatus according to claim 1, wherein the estimated value is updated based on a model in which a time baseline v _baseline (t) is introduced in addition to the spatial baseline a _baseline (x, y). In the parameter setting unit,
Using the transformation matrix based on the impulse response of the calcium concentration change for the spike time series, set the relational expression between the spike time series and the calcium concentration change,
In the calculation unit,
The information processing apparatus according to claim 1, wherein the estimated value is updated based on a transformation matrix. In the parameter setting unit,

In the calculation unit,
The information processing apparatus according to any one of claims 1 to 3, wherein the estimated value is updated based on the set prior distribution. In the initial setting unit,
A number of cell candidates covering a specific area larger than the assumed cell size are arranged regardless of the true cell position so as to fill the entire image with sufficient overlap. Item 5. The information processing device according to any one of Items 1 to 4.   Further, a degenerate cell candidate that removes a degenerated cell candidate with respect to at least one of the output of the spike time series / time direction baseline calculation unit and the cell shape / space direction baseline calculation unit 6. The information processing apparatus according to claim 1, further comprising a removal unit. Further, a similar cell integration unit that integrates similar cell candidates with respect to at least one of the outputs of the spike time series / time direction baseline calculation unit and the cell shape / space direction baseline calculation unit The information processing apparatus according to claim 1, wherein the information processing apparatus includes: The computing unit is
Calcium that minimizes the objective function of MAP estimation under the constraint that the cell shape a _c (x, y) and the spatial baseline a _baseline (x, y) are fixed and the variable value of the spike time series is nonnegative A calcium concentration change / time direction baseline calculation unit for solving a CT + T-dimensional quadratic programming problem to obtain a concentration change v _c (t) and a time baseline v _baseline (t);
Cell shape a _c (which minimizes the objective function of MAP estimation under the constraint that the calcium concentration change v _c (t) and time baseline v _baseline (t) are fixed and the cell shape variable value is non-negative. a cell shape / space direction baseline calculation unit that solves XY (C + 1) -dimensional quadratic programming problems for _obtaining x, y) and a space baseline a _baseline (x, y) The information processing apparatus according to any one of claims 1 to 7. In the calcium concentration change / time direction baseline calculation unit,
The information processing apparatus according to claim 8, wherein the quadratic programming problem is iteratively solved using a main dual interior point method and a conjugate gradient method. The product of the cell shape a _c (x, y) expressed as a non-negative value of cell activity and the calcium concentration change v _c (t) derived from the spike time series u _c (t) expressed as a non-negative value Based on a model that describes observation data f (t, x, y) using a linear sum of multiple cell activities and a baseline, based on Bayesian theory based on observation data, model, and prior distribution An information processing method using an information processing device for obtaining a cell position and activity assumed by estimation,
A data reading step for reading the fluorescence intensity f (t, x, y) at time t at the light emission position (x, y) from the fluorescent protein that occurs as the intracellular calcium concentration increases due to the spike formation of the cell; ,
A parameter setting step for setting a prior expression for the relationship between the spike time series u _c (t) and the calcium concentration change v _c (t) and the cell shape a _c (x, y) and the spike time series u _c (t); ,
A cell candidate having a cell shape a _c (x, y) is arranged based on the fluorescence intensity f (t, x, y) at the time t at the position (x, y) read by the data reading unit, and the base An initial setting step for setting a line based on a statistic of fluorescence intensity f (t, x, y);
Reduce the objective function of the minimization problem derived by the estimation based on the Bayesian theory, and estimate the cell shape a _c (x, y), calcium concentration change v _c (t), and baseline closer to the true solution A calculation step to be updated to
A convergence determination step of determining that the objective function has converged when the change in the objective function is sufficiently small by the update in the calculation step, and ending the processing.   A program for causing a computer to execute the information processing method according to claim 10.

Images