Flaws in the diffusion process

[jlubo]

There seem to be several flaws in the diffusion process (cf. the requirements here):

1. the number of particles does not stay constant in the absence of particle injection and removal
2. the equilibrium concentration depends on the spatial discretization
3. negative concentrations can occur, which is not possible in physical reality

These points can be demonstrated with the attached example (tested with master version d1139b7).

Point 1 is shown in the figures below by the dashed line termed "sum", which is the sum over the concentrations of all compartments weighted by their volume. The sum should exhibit a discontinuous increase upon each particle injection (or decrease upon each removal) and otherwise stay constant. Its final value should be 1500, because of twice an injection of 1000 particles and once a removal of 500 particles. The outcome does not match this assumption.

Point 2 can be seen from comparing the two figures below. For the first figure, the two dendrites each consisted of 3 compartments (CVs). For the second figure, no such fine discretization was used and each dendrite as well as the soma consisted of one compartment only. The final values of the concentrations do not match, which is opposed to the expectation.

Point 3 should become obvious from the top right plot in both figures.

[thorstenhater]

Hi @jlubo,

thanks for reporting. I sadly cannot run the example on current master for two reasons:

• area is being set from external, which is refused since area is built-in (since a few commits back)
• removing the settings result in further errors in the plotting code.

Could try this with current master, potentially the area calculation is off. That might be also causing issue 2.

For the negative values, as far as I understand the issue, this can only happen if a mechanism directly accesses
Xd and doesn't check for underflow? In some interpretations of Xd this might be considered part of the design
(ie a fluctuation around the mean). Is that correct, or am I missing something here? What semantics would you have
expected?

As for the conservation issues: I cannot extract the injected values from the plots. In the second, the deviations seem
close to what you describe (1000, 1000, -500), but in the first they are seemingly closer to 3000, 3000, -1000, which
is also hinting at an issue with area normalisation. It's hard to see, but from the plot, I'd say the deviation of sum
(dashed ./. red lines) exist from about 1.5ms onwards. Thus, I suspect that the calculations for the injected amounts
are slightly off. Whether this is in our code or yours is hard to tell. Also, I'd like to run a convergence study to rule out
numerical issues due to underresolution.

[jlubo]

I've compiled an updated version of the example that is run with the newer master version 235b291. In the discretized (3 CVs per dendrite) case above, the normalization was not correct - my bad! The dendrite volume and area in that case had to be divided by the number of CVs.

I've now run the simulation again, for three different discretizations (the soma always only has 1 CV):

• 1 CV per dendrite
• 3 CVs per dendrite
• 1000 CVs per dendrite

The injected/removed amount of particles now - with the correct normalization - matches the expectation in all cases (as shown in the plots below). I've now added scales to the plots to indicate these changes in the amount of particles. Please note that the sum (black dashed line) is not supposed to match the red line, it's just by incident in the same plot. I would expect the sum to look like this:

Here are the new results:

1 CV per dendrite (fine_discr=1)

3 CVs per dendrite (fine_discr=3)

1000 CVs per dendrite (fine_discr=1000)

The issues mentioned above still seem unresolved.

Regarding the constant number of particles (point 1), the sum in the plot still does not match the expectation. I was wondering if this would be different for finer discretization (for which I haven't computed the sum yet), but from looking at the third figure it doesn't seem to be the case.

Regarding point 2, the discretization still has an (undesired) impact on the equilibrium concentration.

Regarding the negative concentrations (point 3), I guess we might leave this here because it may be left to the user to ensure that not too high an amount of particles is removed. It's also good to note that negative (relative) concentrations may be desirable in some cases, as mentioned by @thorstenhater.

I'm wondering what we can conclude for points 1 & 2 from the results now. High resolution (1000 CVs per dendrite) apparently doesn't solve the issues. I guess something might be wrong with the diffusion across CVs. Especially, I'm wondering why the sum does not match the expectation, and why in the 1000 CVs case the protein concentration converges more slowly.

[jlubo]

I've just updated the code in the previous post following our offline discussion today (the issue remains, though):

• declared ASSIGNED variables as RANGE (to declare that there is one value per CV),
• simplified the point mechanism.

[jlubo]

Here is a little update on the "diffusion confusion", following joint efforts by @thorstenhater, @Shirin1993, and me to nail down the issue.

To emphasize the main issue again: the total amount of particles is in some cases not conserved.

Meanwhile, we have found that at the connection between segments, erroneous adding or removal of particles can occur. This does not seem to happen between segments of the same radius (case a), and neither for diffusion across the CVs of a single segment, but between segments of different radius (case b).

The results below demonstrate the apparent behavior (here is the example code to generate them, tested with the brand new release v0.9.0). For both the case with two segments (i.e., one branch) and three segments (i.e., three branches), there is no issue if the radii are equal (a). But if the radii are different, issues occur (b) In the plots, $\mathrm{s}$ indicates the concentration of diffusive particles in the center of the first segment/"soma", $\mathrm{sd(n)}$ indicates the concentration at the center of different segments, and the expected increases/decreases in the amount of particles are shown by $\mathrm{inc}_0$, $\mathrm{inc}_1$, and $\mathrm{dec}_0$.

a) Same radius

Diffusion across morphology of 1 segments
radius_1 = 4.0 (µm)
length_tot = 10.0 (µm)
area_tot = 251.32741228718345 (µm^2)
area_per_cv = 0.41887902047863906 (µm^2)
volume_tot = 502.6548245743669 (µm^3)
volume_per_cv = 0.8377580409572781 (µm^3)
diffusivity = 1 (m^2/s)
num_cvs_per_seg = 600; len(fs) = 1
Number of branches: 1
Number of CVs: 600 (expected: 600)
lim(s) [estimated] = 0.7957747137822269 (mmol/L) (expected: 0.7957747154594768)
lim(s) [direct] = 0.7957747137822034 (mmol/L) (expected: 0.7957747154594768)
max(s) [direct] = 3.5809946559988144 (mmol/L) (expected: 3.580986219567645)
lim(s⋅V) [estimated] = 399.99999915692223 (1e-18 mol) (expected: 400)
lim(s⋅V) [direct] = 399.9999991569817 (1e-18 mol) (expected: 400)
max(s⋅V) [direct] = 1799.9999999181446 (1e-18 mol) (expected: 1800)

Diffusion across morphology of 2 segments
radius_1 = 4.0 (µm)
radius_2 = 4.0 (µm)
length_tot = 10.0 (µm)
area_tot = 251.32741228718345 (µm^2)
area_per_cv = 0.41887902047863906 (µm^2)
volume_tot = 502.6548245743669 (µm^3)
volume_per_cv = 0.8377580409572781 (µm^3)
diffusivity = 1 (m^2/s)
num_cvs_per_seg = 300; len(fs) = 1
Number of branches: 1
Number of CVs: 600 (expected: 600)
lim(s) [estimated] = 0.7957747139090435 (mmol/L) (expected: 0.7957747154594768)
lim(s) [direct] = 0.7957747139090259 (mmol/L) (expected: 0.7957747154594768)
max(s) [direct] = 3.5809862192155584 (mmol/L) (expected: 3.580986219567645)
lim(s⋅V) [estimated] = 399.99999922066723 (1e-18 mol) (expected: 400)
lim(s⋅V) [direct] = 399.9999992207561 (1e-18 mol) (expected: 400)
max(s⋅V) [direct] = 1799.999999923777 (1e-18 mol) (expected: 1800)

Diffusion across morphology of 3 segments
radius_1 = 4.0 (µm)
radius_2 = 4.0 (µm)
radius_3 = 4.0 (µm)
length_tot = 10.0 (µm)
area_tot = 251.32741228718345 (µm^2)
area_per_cv = 0.41887902047863906 (µm^2)
volume_tot = 502.6548245743669 (µm^3)
volume_per_cv = 0.8377580409572781 (µm^3)
diffusivity = 1 (m^2/s)
num_cvs_per_seg = 200; len(fs) = 1
Number of branches: 3
Number of CVs: 600 (expected: 600)
lim(s) [estimated] = 0.794450629478249 (mmol/L) (expected: 0.7957747154594768)
lim(s) [direct] = 0.7944506294782641 (mmol/L) (expected: 0.7957747154594768)
max(s) [direct] = 3.575027839510368 (mmol/L) (expected: 3.580986219567645)
lim(s⋅V) [estimated] = 399.33444179338454 (1e-18 mol) (expected: 400)
lim(s⋅V) [direct] = 399.3344417934238 (1e-18 mol) (expected: 400)
max(s⋅V) [direct] = 1797.0512842877124 (1e-18 mol) (expected: 1800)

b) Different radius

Diffusion across morphology of 2 segments
radius_1 = 4.0 (µm)
radius_2 = 6.0 (µm)
length_tot = 10.0 (µm)
area_tot = 314.1592653589793 (µm^2)
area_per_cv = 0.5235987755982989 (µm^2)
volume_tot = 816.8140899333463 (µm^3)
volume_per_cv = 1.3613568165555772 (µm^3)
diffusivity = 1 (m^2/s)
num_cvs_per_seg = 300; len(fs) = 1
Number of branches: 1
Number of CVs: 600 (expected: 600)
lim(s) [estimated] = 0.21764778586169836 (mmol/L) (expected: 0.4897075172058318)
lim(s) [direct] = 0.3536776520252521 (mmol/L) (expected: 0.4897075172058318)
max(s) [direct] = 1.5915494314654912 (mmol/L) (expected: 2.203683827426243)
lim(s⋅V) [estimated] = 177.77777813463098 (1e-18 mol) (expected: 400)
lim(s⋅V) [direct] = 288.8888894697973 (1e-18 mol) (expected: 400)
max(s⋅V) [direct] = 1315.8436086902598 (1e-18 mol) (expected: 1800)

Diffusion across morphology of 3 segments
radius_1 = 4.0 (µm)
radius_2 = 6.0 (µm)
radius_3 = 6.0 (µm)
length_tot = 10.0 (µm)
area_tot = 335.1032163829113 (µm^2)
area_per_cv = 0.5585053606381855 (µm^2)
volume_tot = 921.533845053006 (µm^3)
volume_per_cv = 1.53588974175501 (µm^3)
diffusivity = 1 (m^2/s)
num_cvs_per_seg = 200; len(fs) = 1
Number of branches: 3
Number of CVs: 600 (expected: 600)
lim(s) [estimated] = -0.6500050617797218 (mmol/L) (expected: 0.4340589357051691)
lim(s) [direct] = -1.1916759465961795 (mmol/L) (expected: 0.4340589357051691)
max(s) [direct] = 1.5889012631169197 (mmol/L) (expected: 1.953265210673261)
lim(s⋅V) [estimated] = -599.0016638857837 (1e-18 mol) (expected: 400)
lim(s⋅V) [direct] = -1098.1697171259289 (1e-18 mol) (expected: 400)
max(s⋅V) [direct] = 1497.9650884778123 (1e-18 mol) (expected: 1800)

I hope that this documents the whole issue for the moment. As a next step, I want to create a PR with a unit test for the diffusion process. And hopefully soon, we can get the whole thing fixed.

Further notes

• Thorsten has tested the case where the value of the particle concentration is stipulated in a small region (as opposed to injection through a point mechanism), and found that this did not make a difference.
• Setting the coordinates for the segments does sometimes not yield the expected result (see commented block in the code attached to this post).

[thorstenhater]

Hi @jlubo

a few questions

Setting the coordinates for the segments does sometimes not yield the expected result (see commented block in the code attached to this post).

Which one?

Also, some kind of legend in the plots/explanations would be appreciated.

• s
• sd
• ...

From your plots in the 'different radius' section I'd suspect X < 0 for t > 1.5ms. Can you check?

After thinking about what I found before the holidays, my current hypothesis is that we need to scale
concentrations fluxes by the relative volumes in the transport equation. That is, model particle flow across
the CV boundaries. While tempting, we also cannot simply say that Xd models particle amount as the
diffusive flux is driven by the difference in concentration.

[jlubo]

Hi @thorstenhater!

1.

Regarding the code where the segments are set, I meant this:

s = tree.append(A.mnpos,
		          A.mpoint(-length/3, 0, 0, radius_1),
		          A.mpoint(0, 0, 0, radius_1),
		          tag=0)
_ = tree.append(s,
		          A.mpoint(0, 0, 0, radius_2),
		          A.mpoint(length/3, 0, 0, radius_2),
		          tag=1)
_ = tree.append(s,
		          A.mpoint(-length/3, 0, 0, radius_3),
		          A.mpoint(-2*length/3, 0, 0, radius_3),
		          tag=2)
''' strangely, this does not work:
s = tree.append(A.mnpos,
		          A.mpoint(-length/6, 0, 0, radius_1),
		          A.mpoint(+length/6, 0, 0, radius_1),
		          tag=0)
_ = tree.append(s,
		          A.mpoint(-length/6, 0, 0, radius_2),
		          A.mpoint(-length/3, 0, 0, radius_2),
		          tag=1)
_ = tree.append(s,
		          A.mpoint(+length/6, 0, 0, radius_3),
		          A.mpoint(+length/3, 0, 0, radius_3),
		          tag=2)
'''

Maybe that would be another issue, though.

2.

For the legend, I've added a description in the previous post above the plots.

3.

It's right, the concentration becomes negative in the case with 3 segments.

4.

After thinking about what I found before the holidays, my current hypothesis is that we need to scale
concentrations fluxes by the relative volumes in the transport equation. That is, model particle flow across
the CV boundaries. While tempting, we also cannot simply say that Xd models particle amount as the
diffusive flux is driven by the difference in concentration.

I agree, that sounds like a possible solution. Especially, smoothing the transition between different radii (see red line in the sketch below) could be useful, as we've discussed offline.

[thorstenhater]

Re 1: What does 'doesn't work mean'? Gives an error?
Semantically, these coordinates you give have no impact on the location where the segment attaches, ie

p = tree.append(p, (0,0,0,1), (0,0,1,1), 1)
p = tree.append(p, (0,0,1,1), (0,0,2,1), 1)

is the same as

p = tree.append(p, (0,0,0,1), (0,0,1,1), 1)
p = tree.append(p, (0,0,1000,1), (0,0,1001,1), 1)

[jlubo]

is the same as

That's what I thought, too. But the second setting (the commented one) does not yield the correct results even with equal radii of the segments.

[thorstenhater]

About the coordinates issue:

# length = l/3
s = tree.append(A.mnpos,
		          A.mpoint(-length/3, 0, 0, radius_1),
		          A.mpoint(0, 0, 0, radius_1),
		          tag=0)
# length = l/3
_ = tree.append(s,
		          A.mpoint(0, 0, 0, radius_2),
		          A.mpoint(length/3, 0, 0, radius_2),
		          tag=1)
# length = l/3
_ = tree.append(s,
		          A.mpoint(-length/3, 0, 0, radius_3),
		          A.mpoint(-2*length/3, 0, 0, radius_3),
		          tag=2)

Whereas

# length = l/3
s = tree.append(A.mnpos,
		          A.mpoint(-length/6, 0, 0, radius_1),
		          A.mpoint(+length/6, 0, 0, radius_1),
		          tag=0)
# length = l/6
_ = tree.append(s,
		          A.mpoint(-length/6, 0, 0, radius_2),
		          A.mpoint(-length/3, 0, 0, radius_2),
		          tag=1)
# length = l/6
_ = tree.append(s,
		          A.mpoint(+length/6, 0, 0, radius_3),
		          A.mpoint(+length/3, 0, 0, radius_3),
		          tag=2)

So, I wouldn't be too surprised about different results.

[jlubo]

Absolutely right, thanks for redoing the simple math 🙈 Using

	s = tree.append(A.mnpos,
		            A.mpoint(-length/6, 0, 0, radius_1),
		            A.mpoint(+length/6, 0, 0, radius_1),
		            tag=0)
	_ = tree.append(s,
		            A.mpoint(-length/6, 0, 0, radius_2),
		            A.mpoint(-length/2, 0, 0, radius_2),
		            tag=1)
	_ = tree.append(s,
		            A.mpoint(+length/6, 0, 0, radius_3),
		            A.mpoint(+length/2, 0, 0, radius_3),
		            tag=2)

yields just the expected result.

[jlubo]

I have just created PR #2209, which basically adds the tests that are shown in the most recent plots here in this issue.