I was notified that there might be a bug in the LaminarCurrentSourceDensity class, where a user had the problem that the CSD did not sum to 0 at each time step, as one would expect for equal volumes that contain all the membrane currents.
The user sent me a slightly modified version of example_1, and I was able to reproduce the problem with the Mainen Sejnowski cell model, but not really understand the cause. Below is a minimal working example with a simple cell:
import numpy as np
import matplotlib.pyplot as plt
from lfpykit import CellGeometry, LaminarCurrentSourceDensity
cell = CellGeometry(x=np.array([[50, 50], [50, -50]]),
y=np.array([[0, 0], [0, 0]]),
z=np.array([[-100, -50], [-50, 50]]),
d=np.array([1, 1]))
I_m = np.array([[1.,], [-1]])
z = np.array([[-150., 0.], [0., 150.]])
r = np.array([200., 200.,])
# instantiate electrode, get linear response matrix
csd_lam = LaminarCurrentSourceDensity(cell=cell, z=z, r=r)
M = csd_lam.get_transformation_matrix()
csd = M @ I_m
print(f"CSD:\n{csd}\n")
print(f"Sum of CSD:\n{np.max(np.sum(csd,axis=0))}\n")
print(f"Sum of Imem:\n{np.max(np.sum(I_m,axis=0))}\n")
ax1 = plt.subplot(111, aspect=1, xlim=[-200, 200])
ax1.plot(cell.x.T, cell.z.T)
for i_, z_ in enumerate(z):
ax1.axhspan(z_[0], z_[1], alpha=0.3, fc=plt.cm.viridis(i_ / z.shape[0]))
This gives the output:
CSD:
[[ 3.42950578e-08]
[-1.87565899e-08]]
Sum of CSD:
1.5538467857419034e-08
Sum of Imem:
0.0
If the cell is shifted along the x-axis, this should not change anything, but the CSD instead sum to zero in this case:
cell = CellGeometry(x=np.array([[50, 50], [50, -50]]) - 50,
y=np.array([[0, 0], [0, 0]]),
z=np.array([[-100, -50], [-50, 50]]),
d=np.array([1, 1]))
gives:
CSD:
[[ 2.65258238e-08]
[-2.65258238e-08]]
Sum of CSD:
0.0
Sum of Imem:
0.0
I'm quite sure the problem is that the Pr is not commputed correctly
|
Pr, Pz, hit = _PrPz(r0=R[k, 0], z0=self.cell.z[k, 0], |
and that this is related to that differences in the radial coordinate does not necessarily reflect cartesian lengths. In the example cell model, the dendrite extends from x=50 to x=-50, which has the same radial distance from the center, but is still a distance of 100 µm. This is why the shift works, since now the radial distance also changes. Or something like that.
A potentially useful diagnostics and future unit test is
V = np.pi * r**2 * np.diff(z, axis=-1).flatten() # Volume
fraction_matrix = M.T * V
print(fraction_matrix)
which, for the cell that gives the error is:
[[1. 0. ]
[0.35355339 0.35355339]]
I think each row here should sum to one, as is reflects how currents should be distributed between the CSD regions.
It shouldn't be too hard to find a fix, but I will not have time to look at this in the next few days, so if anybody feels inclined to fix it before that, please go ahead!
I was notified that there might be a bug in the LaminarCurrentSourceDensity class, where a user had the problem that the CSD did not sum to 0 at each time step, as one would expect for equal volumes that contain all the membrane currents.
The user sent me a slightly modified version of example_1, and I was able to reproduce the problem with the Mainen Sejnowski cell model, but not really understand the cause. Below is a minimal working example with a simple cell:
This gives the output:
CSD:
[[ 3.42950578e-08]
[-1.87565899e-08]]
Sum of CSD:
1.5538467857419034e-08
Sum of Imem:
0.0
If the cell is shifted along the x-axis, this should not change anything, but the CSD instead sum to zero in this case:
gives:
CSD:
[[ 2.65258238e-08]
[-2.65258238e-08]]
Sum of CSD:
0.0
Sum of Imem:
0.0
I'm quite sure the problem is that the Pr is not commputed correctly
LFPykit/lfpykit/models.py
Line 2084 in 8980988
and that this is related to that differences in the radial coordinate does not necessarily reflect cartesian lengths. In the example cell model, the dendrite extends from x=50 to x=-50, which has the same radial distance from the center, but is still a distance of 100 µm. This is why the shift works, since now the radial distance also changes. Or something like that.
A potentially useful diagnostics and future unit test is
which, for the cell that gives the error is:
[[1. 0. ]
[0.35355339 0.35355339]]
I think each row here should sum to one, as is reflects how currents should be distributed between the CSD regions.
It shouldn't be too hard to find a fix, but I will not have time to look at this in the next few days, so if anybody feels inclined to fix it before that, please go ahead!