Inverse Current-Source Density Method in 3D: Reconstruction Fidelity, Boundary Effects, and Influence of Distant Sources

Abstract

Estimation of the continuous current-source density in bulk tissue from a finite set of electrode measurements is a daunting task. Here we present a methodology which allows such a reconstruction by generalizing the one-dimensional inverse CSD method. The idea is to assume a particular plausible form of CSD within a class described by a number of parameters which can be estimated from available data, for example a set of cubic splines in 3D spanned on a fixed grid of the same size as the set of measurements. To avoid specificity of particular choice of reconstruction grid we add random jitter to the points positions and show that it leads to a correct reconstruction. We propose different ways of improving the quality of reconstruction which take into account the sources located outside the recording region through appropriate boundary treatment. The efficiency of the traditional CSD and variants of inverse CSD methods is compared using several fidelity measures on different test data to investigate when one of the methods is superior to the others. The methods are illustrated with reconstructions of CSD from potentials evoked by stimulation of a bunch of whiskers recorded in a slab of the rat forebrain on a grid of 4×5×7 positions.

Appendices

Appendices

Linear iCSD Method

In this appendix we construct the F matrix for the linear distribution of CSD. Consider a grid of points (i,j,k), where i = 1..n _x , j = 1..n _y , k = 1..n _z . Let us number the points with a multi-index α ≡ (i,j,k) and write (i,j,k) ≡ (x _α ,y _α ,z _α ). Let V be the set of 8 vectors (v ₁,v ₂,v ₃), $v_i \in \lbrace 0,1 \rbrace$. The grid has n = n _x n _y n _z nodes and there are m = (n _x  − 1)(n _y  − 1)(n _z  − 1) boxes spanned by these nodes. We index the boxes so that the corners of the box number α are α + v for v ∈ V. Let B denote the set of all the allowed indices α numbering the boxes and G stand for all the grid points. Let C denote the vector of CSD values at the nodes, that is C _α  = C(x _α ,y _α ,z _α ) for α ∈ G. We assume that CSD in boxes is given by a linear approximation.

Take a point (x,y,z) in box number α and let δx = x − x _α , δy = y − y _α , δz = z − z _α . The value of CSD at this point obtained with the linear interpolation is given by:

$$C\left( {x,y,z} \right) = \sum\limits_{\upsilon \in V} {\left[ {1 - \upsilon _1 + \left( {2\upsilon _1 - 1} \right)\delta x} \right] \times \left[ {1 - \upsilon _2 + \left( {2\upsilon _2 - 1} \right)\delta y} \right] \times \left[ {1 - \upsilon _3 + \left( {2\upsilon _3 - 1} \right)\delta z} \right]C_{\alpha + \upsilon } .} $$

Therefore the distribution inside the box is a linear combination of 8 functions f _l , l = 1..8: f ₁ = 1, f ₂ = δx, f ₃ = δy, f ₄ = δz, f ₅ = δx δy, f ₆ = δx δz, f ₇ = δy δz, f ₈ = δx δy δz, with coefficients depending linearly on the values of C at the nodes of the box:

$$C\left( {x,y,z} \right) = \sum\limits_{\beta \in G} {\sum\limits_{l = 1}^8 {E_{\alpha \beta }^1 f_l C_\beta } } .$$

The coefficients $E^l_{\alpha\beta}$ are non-zero only for β − α ∈ V and follow from the above formula, e.g. $E^1_{\alpha,\beta}=1$ for β − α = (0,0,0), otherwise $E^1_{\alpha,\beta}=0$, etc. The potential generated by such a distribution of current-source density at some point $(\tilde{x},\tilde{y},\tilde{z})$ is

$$ \Phi(\tilde{x},\tilde{y},\tilde{z}) = \sum_{\alpha\in B}\sum_{\beta\in G}\sum_{l=1}^8 F^l_\alpha(\tilde{x},\tilde{y},\tilde{z}) E^l_{\alpha\beta} C_\beta \;, $$

where

$$ F^l_\alpha(\tilde{x},\tilde{y},\tilde{z})= \frac{1}{4\pi\sigma} \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \frac{f_l(x,y,z) \,\, \mathrm{d} z \, \mathrm{d} y \, \mathrm{d} x} {\sqrt{(\tilde{x}-x_\alpha-x)^2 + (\tilde{y}-y_\alpha-y)^2 + (\tilde{z}-z_\alpha-z)^2}} \;. $$

If we now take as $(\tilde{x},\tilde{y},\tilde{z})$ one of the grid points γ then

$$ \Phi_\gamma = \sum_{\beta\in G} F_{\gamma\beta} C_\beta, $$

where

$$F_{\gamma\beta} = \sum_{\alpha\in B}\sum_{l=1}^8 F^l_\alpha(x_\gamma, y_\gamma,z_\gamma) E^l_{\alpha\beta}. $$

Thus F _γβ represents the direct and indirect contributions to the total potential at point γ from the CSD associated with the grid point β.

Spline iCSD Method

The construction of the F matrix for the spline distribution is in principle very similar to the case of linear distribution, but the level of complication is much higher. Similarly as in the previous appendix we have n nodes and m boxes, but now the interpolating function in each box is the three-dimensional cubic spline. That means there are 4×4×4 = 64 base functions. Therefore there will be 64 each of E ^ℓ and F ^ℓ matrices. The main challenge here is to construct the E ^ℓ matrices, that means to describe how a CSD value associated with a given node influences the interpolating splines in all the boxes.

First we briefly remind the construction of one-dimensional spline (Press et al. 1992). Suppose we have values of a function f at points x = 1 .. n _x . For x such that j ≤ x ≤ j + 1 define P ₁(x) = j + 1 − x, P ₂(x) = x − j. The formula

$$ \label{jonel} f(x) = P_1(x) f(j) + P_2(x) f(j+1) $$

(7)

gives a linear approximation, that means an approximation with a continuous function. In case of cubic splines we want more: we want also the first and second derivatives to be continuous. This can be done if we replace the right hand side of Eq. 7 with a third-degree polynomial:

$$ \label{eq8} f(x) = P_1(x) f(j{\kern1.2pt}) + P_2(x) f(j+1) + P_3(x) f^{\prime\prime}(j{\kern1.2pt}) + P_4(x) f^{\prime\prime}(j+1) \;, $$

(8)

where $P_3(x) = \frac16 (P_1(x)^3 -P_1(x))$, $P_4(x) = \frac16 (P_2(x)^3 - P_2(x))$. This formula guarantees that both f and its second derivative are continuous. The values of the second derivative f″ of f at nodes are a priori unknown, but can be calculated from the condition that also the first derivative is continuous. This condition leads to the following system of equations (j = 2.. n _x  − 1):

$$\frac{1}{6}f\prime \prime \left( {j - 1} \right) + \frac{1}{3}f\prime \prime \left( j \right) + \frac{1}{6}f\prime \prime \left( {j + 1} \right) = f\left( {j + 1} \right) - 2f\left( j \right) + f\left( {j - 1} \right).$$

There are n unknown f ^″(j) and only n − 2 equations, therefore we need two additional conditions. There are several commonly used choices for these conditions. One of them is to add equations

$$ f^{\prime\prime}(1) = 0\;, \qquad f^{\prime\prime}(n) = 0\;, $$

which leads to so-called “natural splines”. The splines implemented in Matlab use a different set of conditions, namely

$$f\prime \prime \left( 3 \right) - f\prime \prime \left( 2 \right) = f\prime \prime \left( 2 \right) - f\prime \prime \left( 1 \right),\,f\prime \prime \left( n \right) - f\prime \prime \left( {n - 1} \right) = f\prime \prime \left( {n - 1} \right) - f\prime \prime \left( {n - 2} \right).$$

These conditions mean that the third derivative is continuous at x = 2 and x = n − 1 and they lead to what is known as “not-a-knot” splines. The important thing is that in both cases the values of f ^″ at the nodes can be obtained from f(j), j = 1.. n, by a linear operation which we call G:

$$ f^{\prime\prime}(i) = \sum_{j=1}^n G_{ij} f(j)\;. $$

The matrix G is different for “natural” and “not-a-knot” splines.

The three-dimensional spline interpolation is obtained simply by performing three one-dimensional splines. The complication is that we do not want the values of the interpolating function at some points, but the coefficients standing by the base functions. We found that it is convenient to choose base functions which are products of the polynomials P ₁, P ₂, P ₃, P ₄ of variables x, y and z, that means P ₁(x)P ₁(y)P ₁(z), P ₁(x)P ₁(y)P ₂(z), ..., P ₄(x)P ₄(y)P ₄(z). To extract the coefficients we start with the spline in z direction:

$$ \label{spl1} f(x,y,z) = P_1(z) f(x,y,j) + P_2(z) f(x,y,j+1) + P_3(z) f_{zz}(x,y,j) + P_4(z) f_{zz}(x,y,j+1)\;, $$

(9)

where f _zz stands for $\frac{\partial^2f}{\partial z^2}$ and is given by $f_{zz}(x,y,j{\kern1.2pt}) = \sum_{i=1}^{n_z} G^z_{ji} f(x,y,i)$. Therefore we reduce the problem to two-dimensional splines in the xy-plane. We continue with

$$ \label{spl2} f(x,y,j\,) = P_1(y) f(x,i,j\,) + P_2(y) f(x,i+1,j) + P_3(y) f_{yy}(x,i,j\,) + P_4(z) f_{yy} (x,i+1,j\,)\;, $$

(10)

and so on. In the end we get the coefficients standing by the base functions as combinations of f(i,j,k) (values of f at the nodes) and the matrices G ^x, G ^y, G ^z. Then we construct the matrices $E^{pqr}_{\alpha\beta}$, α ∈ B, β ∈ G, 1 ≤ p,q,r ≤ 4. The number $E^{pqr}_{\alpha\beta}$ is the coefficient standing by P _p (x)P _q (y)P _r (z) in box α resulting from a unit CSD at the node β. The construction of 64 $F^{pqr}_{\gamma\alpha}$ matrices (each of size n by m), where γ ∈ G and α ∈ B, is simple:

$$F_{\gamma \alpha }^{pqr} = \frac{1}{{4\pi \sigma }}\int_0^1 {\int_0^1 {\int_0^1 {\frac{{P_p \left( x \right)P_q \left( y \right)\operatorname{P} _r \left( z \right)}}{{\sqrt {\left( {x_\gamma - x_\alpha - x} \right)^2 + \left( {y_\gamma - y_\alpha - y} \right)^2 + \left( {z_\gamma - z_\alpha - z} \right)^2 } }}dzdydx,} } } $$

and the full matrix F is now

$$ F = \sum_{p,q,r=1}^4 F^{pqr} E^{pqr} \;, $$

or

$$ F_{\gamma\beta} = \sum_{p,q,r=1}^4 \sum_{\alpha=1}^m F^{pqr}_{\gamma\alpha} E^{pqr}_{\alpha\beta} \;. $$