Inject/decay of diffusing particles does not work with custom mechanism

[jlubo]

Injecting diffusing particles with the built-in mechanism inject and letting their concentration decay with the built-in mechanism decay works fine. But if, instead of the built-in mechanism inject, an identical mechanism custom_inject from a custom catalogue is used, diffusion does not work anymore.

Please see the attached plots and code: diffusion_with_inject_and_decay_2023-06-02.zip. The example is somewhat similar to the existing diffusion example. I derived it from Shirin's code for calcium diffusion. I've tested this with Arbor v0.8.1.

[thorstenhater]

That would be extremely strange, but I'll take a look. Did you try adding custom_inject to the vanilla example?

[jlubo]

For the vanilla example the results differ as well:

results_builtin.pdf results_custom.pdf

[thorstenhater]

Ok, that is curious. Just to double check: Are both build using the same compiler flags? That is in particular: ARB_VECTORIZE?

[jlubo]

Yes, the only thing that is changed is the mechanism, i.e.

mech_inject = arbor.mechanism('inject/x=test_ion', {'alpha' : 100})

is exchanged with

mech_inject = arbor.mechanism('custom_inject/x=test_ion', {'alpha' : 100})

[Shirin1993]

Here is another version of the problem : I tried a custom injection in combination with a custom decay and also a built-in decay. It seems that :

• the built-in decay doesnt work, since it results in zero concentration.

• the custom decay also doesnt work since changing the parameter tau has no effect on diffusion.

I have attached the related codes and plots

inj_decay_in_2mod_files.zip

[thorstenhater]

I am trying to reproduce this locally (M1 MacBook Pro).

Baseline

Debug build SIMD=ON; with the defaultdiffusion.py

|----------+----------+----------+----------+----------+----------+----------+----------+----------|
| Time/ms  | 0.0      | 0.125    | 0.25     | 0.375    | 0.5      | 0.625    | 0.75     | 0.875    |
|----------+----------+----------+----------+----------+----------+----------+----------+----------|
|    0.000 |  100.000 |   60.400 |    1.000 |    1.000 |   36.368 |    1.000 |    1.000 |    1.000 |
|    0.100 |   59.804 |   36.976 |    1.049 |    0.866 |   21.526 |    0.863 |    0.608 |    0.606 |
|    0.200 |   35.783 |   22.611 |    0.902 |    0.678 |   12.746 |    0.671 |    0.372 |    0.367 |
|    0.300 |   21.420 |   13.813 |    0.706 |    0.502 |    7.549 |    0.492 |    0.228 |    0.223 |
|    0.400 |   12.828 |    8.430 |    0.524 |    0.357 |    4.473 |    0.348 |    0.141 |    0.135 |
|----------+----------+----------+----------+----------+----------+----------+----------+----------|

(I did, however, reformat the table to fit here)

Replacing the catalogue

@@ -11,7 +16,7 @@ class recipe(A.recipe):
         self.the_cell = cell
         self.the_probes = probes
         self.the_props = A.neuron_cable_properties()
-        self.the_props.catalogue = A.default_catalogue()
+        self.the_props.catalogue = A.load_catalogue('local-catalogue.so')

     def num_cells(self):
         return 1

where the local catalogue just has copies of the original mechanisms.

Debug build, SIMD=ON

| Time/ms  | 0.0      | 0.125    | 0.25     | 0.375    | 0.5      | 0.625    | 0.75     | 0.875    |
|----------+----------+----------+----------+----------+----------+----------+----------+----------|
|    0.000 |  100.000 |   60.400 |    1.000 |    1.000 |   36.368 |    1.000 |    1.000 |    1.000 |
|    0.100 |   59.804 |   36.976 |    1.049 |    0.866 |   21.526 |    0.863 |    0.608 |    0.606 |
|    0.200 |   35.783 |   22.611 |    0.902 |    0.678 |   12.746 |    0.671 |    0.372 |    0.367 |
|    0.300 |   21.420 |   13.813 |    0.706 |    0.502 |    7.549 |    0.492 |    0.228 |    0.223 |
|    0.400 |   12.828 |    8.430 |    0.524 |    0.357 |    4.473 |    0.348 |    0.141 |    0.135 |
|----------+----------+----------+----------+----------+----------+----------+----------+----------|

Release build, SIMD=ON

|----------+----------+----------+----------+----------+----------+----------+----------+----------|
| Time/ms  | 0.0      | 0.125    | 0.25     | 0.375    | 0.5      | 0.625    | 0.75     | 0.875    |
|----------+----------+----------+----------+----------+----------+----------+----------+----------|
|    0.000 |  100.000 |   60.400 |    1.000 |    1.000 |   36.368 |    1.000 |    1.000 |    1.000 |
|    0.100 |   59.804 |   36.976 |    1.049 |    0.866 |   21.526 |    0.863 |    0.608 |    0.606 |
|    0.200 |   35.783 |   22.611 |    0.902 |    0.678 |   12.746 |    0.671 |    0.372 |    0.367 |
|    0.300 |   21.420 |   13.813 |    0.706 |    0.502 |    7.549 |    0.492 |    0.228 |    0.223 |
|    0.400 |   12.828 |    8.430 |    0.524 |    0.357 |    4.473 |    0.348 |    0.141 |    0.135 |
|----------+----------+----------+----------+----------+----------+----------+----------+----------|

Release, SIMD=OFF

|----------+----------+----------+----------+----------+----------+----------+----------+----------|
| Time/ms  | 0.0      | 0.125    | 0.25     | 0.375    | 0.5      | 0.625    | 0.75     | 0.875    |
|----------+----------+----------+----------+----------+----------+----------+----------+----------|
|    0.000 |  100.000 |   60.400 |    1.000 |    1.000 |   36.368 |    1.000 |    1.000 |    1.000 |
|    0.100 |   59.804 |   36.976 |    1.049 |    0.866 |   21.526 |    0.863 |    0.608 |    0.606 |
|    0.200 |   35.783 |   22.611 |    0.902 |    0.678 |   12.746 |    0.671 |    0.372 |    0.367 |
|    0.300 |   21.420 |   13.813 |    0.706 |    0.502 |    7.549 |    0.492 |    0.228 |    0.223 |
|    0.400 |   12.828 |    8.430 |    0.524 |    0.357 |    4.473 |    0.348 |    0.141 |    0.135 |
|----------+----------+----------+----------+----------+----------+----------+----------+----------|

Using the default and the copied catalogue

@@ -11,7 +16,7 @@ class recipe(A.recipe):
         self.the_cell = cell
         self.the_probes = probes
         self.the_props = A.neuron_cable_properties()
-        self.the_props.catalogue = A.default_catalogue()
+        self.the_props.catalogue.extend(A.load_catalogue('local-catalogue.so'), 'local-')

     def num_cells(self):
         return 1
@@ -38,7 +43,7 @@ _ = tree.append(s, A.mpoint(3, 0, 0, 1), A.mpoint(33, 0, 0, 1), tag=3)

 dec = A.decor()
 dec.set_ion("na", int_con=0.0, diff=0.005)
-dec.place("(location 0 0.5)", A.synapse("inject/x=na", {"alpha": 200.0}), "Zap")
+dec.place("(location 0 0.5)", A.synapse("local-inject/x=na", {"alpha": 200.0}), "Zap")
 dec.paint("(all)", A.density("decay/x=na"))
 dec.discretization(A.cv_policy("(max-extent 5)"))

Only one experiment here: Release, SIMD=OFF

|----------+----------+----------+----------+----------+----------+----------+----------+----------|
| Time/ms  | 0.0      | 0.125    | 0.25     | 0.375    | 0.5      | 0.625    | 0.75     | 0.875    |
|----------+----------+----------+----------+----------+----------+----------+----------+----------|
|    0.000 |  100.000 |   60.400 |    1.000 |    1.000 |   36.368 |    1.000 |    1.000 |    1.000 |
|    0.100 |   59.804 |   36.976 |    1.049 |    0.866 |   21.526 |    0.863 |    0.608 |    0.606 |
|    0.200 |   35.783 |   22.611 |    0.902 |    0.678 |   12.746 |    0.671 |    0.372 |    0.367 |
|    0.300 |   21.420 |   13.813 |    0.706 |    0.502 |    7.549 |    0.492 |    0.228 |    0.223 |
|    0.400 |   12.828 |    8.430 |    0.524 |    0.357 |    4.473 |    0.348 |    0.141 |    0.135 |
|----------+----------+----------+----------+----------+----------+----------+----------+----------|

Conclusion

I cannot reproduce this here. Are you sure all the parameters are identical?

[jlubo]

Here is the output of the vanilla example (diffusion.py) in in Arbor v0.9.0.

Built-in catalogue:

Configuration: {'mpi': False, 'mpi4py': False, 'gpu': None, 'vectorize': False, 'profiling': False, 'neuroml': True, 'bundled': True, 'version': '0.9.0', 'source': '2023-08-10T11:47:38+02:00 217c77686e3674adc101db34f947d32785bc788c', 'build_config': 'RELEASE', 'arch': 'native', 'prefix': '/home/jlubo/arbor_v0.9.0', 'python_lib_path': '', 'binary_path': 'bin', 'lib_path': 'lib', 'data_path': 'share', 'CXX': '/usr/bin/c++', 'pybind-version': '2.10.1', 'timestamp': 'Aug 10 2023 16:53:29'}
Sodium concentration (NaD/mM)
|----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------|
| Time (ms)            | (cable 0 0 0.125)    | (cable 0 0.125 0.25) | (cable 0 0.25 0.375) | (cable 0 0.375 0.5)  | (cable 0 0.5 0.625)  | (cable 0 0.625 0.75) | (cable 0 0.75 0.875) | (cable 0 0.875 1)    |
|----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------|
|                0.000 |              100.000 |               60.400 |                1.000 |                1.000 |               36.368 |                1.000 |                1.000 |                1.000 |
|                0.100 |               59.804 |               36.976 |                1.049 |                0.866 |               21.526 |                0.863 |                0.608 |                0.606 |
|                0.200 |               35.783 |               22.611 |                0.902 |                0.678 |               12.746 |                0.671 |                0.372 |                0.367 |
|                0.300 |               21.420 |               13.813 |                0.706 |                0.502 |                7.549 |                0.492 |                0.228 |                0.223 |
|                0.400 |               12.828 |                8.430 |                0.524 |                0.357 |                4.473 |                0.348 |                0.141 |                0.135 |
|----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------|

Custom catalogue:

Configuration: {'mpi': False, 'mpi4py': False, 'gpu': None, 'vectorize': False, 'profiling': False, 'neuroml': True, 'bundled': True, 'version': '0.9.0', 'source': '2023-08-10T11:47:38+02:00 217c77686e3674adc101db34f947d32785bc788c', 'build_config': 'RELEASE', 'arch': 'native', 'prefix': '/home/jlubo/arbor_v0.9.0', 'python_lib_path': '', 'binary_path': 'bin', 'lib_path': 'lib', 'data_path': 'share', 'CXX': '/usr/bin/c++', 'pybind-version': '2.10.1', 'timestamp': 'Aug 10 2023 16:53:29'}
Sodium concentration (NaD/mM)
|----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------|
| Time (ms)            | (cable 0 0 0.125)    | (cable 0 0.125 0.25) | (cable 0 0.25 0.375) | (cable 0 0.375 0.5)  | (cable 0 0.5 0.625)  | (cable 0 0.625 0.75) | (cable 0 0.75 0.875) | (cable 0 0.875 1)    |
|----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------|
|                0.000 |              100.000 |               60.400 |                1.000 |                1.000 |               36.368 |                1.000 |                1.000 |                1.000 |
|                0.100 |               59.804 |               36.976 |                1.049 |                0.866 |               21.526 |                0.863 |                0.608 |                0.606 |
|                0.200 |               35.783 |               22.611 |                0.902 |                0.678 |               12.746 |                0.671 |                0.372 |                0.367 |
|                0.300 |               21.420 |               13.813 |                0.706 |                0.502 |                7.549 |                0.492 |                0.228 |                0.223 |
|                0.400 |               12.828 |                8.430 |                0.524 |                0.357 |                4.473 |                0.348 |                0.141 |                0.135 |
|----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------+----------------------|

So now this example seems to work. Maybe something in v0.9.0 was changed? Which version were you using @thorstenhater ? And nevertheless, for the original example in this issue, the problem persists...

[thorstenhater]

@Shirin1993 This is the current main branch and I used the diffusion example in the simplest variation reported to exhibit the issue. Can we conclude, that some alteration of inject vs the default/inject must be to be blame or the specifics of its use?

[jlubo]

Okay. Regarding the persisting problem, yeah, I guess some use of inject is an issue here. I've tested the original example both with and without SIMD (in Arbor v0.9.0, as I said) and I get the same result:

Built-in catalogue

Custom catalogue

[jlubo]

Might somehow the use of decor.paint be required for the diffusive ion species?

Because after adding decor.paint("(all)", ion_name="test_ion", int_con=100.0, diff=1e-5) I get this:

Built-in catalogue

Custom catalogue

[thorstenhater]

That would indicate that there is a difference between those two examples apart from the mechanisms used?! Or do you suggest that the custom catalogue adds the requirement of that call to paint? Still, the behaviour measured for spine is still different?!

[jlubo]

Yeah, the mechanisms are exactly the same.

[thorstenhater]

So, what is the delta, then?

[thorstenhater]

Now. An extended discussion with @jlubo hinted at the problem and I finally was able to fix this. Please test #2221

[jlubo]

Nice, now I get the expected results! With https://github.com/arbor-sim/arbor/pull/2221, the first example works as it is:

Built-in catalogue

Custom catalogue

For Shirin's code, the custom decay mechanism must be replaced by the new decay mechanism as in https://github.com/arbor-sim/arbor/pull/2221:

NEURON {
    SUFFIX calcium_decay
    USEION ca WRITE cad
    RANGE tau
}

PARAMETER { tau = 5 }

BREAKPOINT {
   SOLVE dca METHOD cnexp
}

DERIVATIVE dca {
   cad' = -tau*cad
}

With that it works!

And the vanilla example still works as well (strangely, already without replacing the custom decay mechanism by the new one from https://github.com/arbor-sim/arbor/pull/2221, also see here).

[thorstenhater]

Solved by #2221