Increase max number of ions per compartment

[wthun]

There's a relatively low limit (128) on the number of ions (and consequently RxD species) due to using ```long´´´ for bitvector comparisons. This PR extends this to 10000 (or an arbitrary number chosen at compile time) by using std::bitset.

This approach tried to minimize other changes. I guess there could be a dynamic limit by using vectors of booleans or the dynamic bitset from boost instead.

Rationale: Issue #2947 (comment)

Minimal test:

import neuron
from neuron import h
import neuron.crxd as rxd

soma = h.Section(name="soma")
soma.L = soma.diam = 10

species_list = []

print("Creating regions")
region = rxd.Region(soma, nrn_region="i")

print("Creating species")
for idx in range(0,200):
    species_name = f"species{idx}"
    spec = rxd.Species(region,
                       d=0,
                       initial=1,
                       charge=0,
                       name=species_name)
    species_list.append(spec)

print("Inserting the first channel writing a particular species (previously crashed due to length limit)")
soma.insert("caldyn_ms")


print("Initializes and runs...")
h.finitialize()
neuron.run(100)

print("Inserting second mechanism writing the same ion, this should produce a warning: \n")
soma.insert("sneaky_caldyn_ms")

caldyn_ms.mod

TITLE Calcium dynamics for L and T calcium pool

NEURON {
    SUFFIX caldyn_ms
    USEION cal READ ical, cali WRITE cali VALENCE 2
    RANGE pump, cainf, taur, drive, depth
}

UNITS {
    (molar) = (1/liter) 
    (mM) = (millimolar)
    (um) = (micron)
    (mA) = (milliamp)
    (msM) = (ms mM)
    FARADAY = (faraday) (coulomb)
}

PARAMETER {
    drive = 10000 (1)
    depth = 0.2  (um)
    cainf = 70e-6 (mM)
    taur = 43 (ms)
    kt = 1e-4 (mM/ms)
    kd = 1e-4 (mM)
    pump = 0.02
}

STATE { cali (mM) }

INITIAL { cali = cainf }

ASSIGNED {
    ical (mA/cm2)
    drive_channel (mM/ms)
    drive_pump (mM/ms)
}
    
BREAKPOINT {
    SOLVE state METHOD cnexp
}


DERIVATIVE state { 
    
    : force concentration to stay above cainf by only pumping if larger
    drive_channel = -drive*ical/(2*FARADAY*depth)
    drive_pump    = -kt*(cali-cainf)/(cali+kd)
    
    if (drive_channel <= 0.) { drive_channel = 0. }
    
    cali' = drive_channel + pump*drive_pump + (cainf-cali)/taur
}

sneaky_caldyn_ms.mod

TITLE Calcium dynamics for L and T calcium pool

NEURON {
    SUFFIX sneaky_caldyn_ms
    USEION cal READ ical, cali WRITE cali VALENCE 2
    RANGE pump, cainf, taur, drive, depth
}

UNITS {
    (molar) = (1/liter) 
    (mM) = (millimolar)
    (um) = (micron)
    (mA) = (milliamp)
    (msM) = (ms mM)
    FARADAY = (faraday) (coulomb)
}

PARAMETER {
    drive = 10000 (1)
    depth = 0.2  (um)
    cainf = 70e-6 (mM)
    taur = 43 (ms)
    kt = 1e-4 (mM/ms)
    kd = 1e-4 (mM)
    pump = 0.02
}

STATE { cali (mM) }

INITIAL { cali = cainf }

ASSIGNED {
    ical (mA/cm2)
    drive_channel (mM/ms)
    drive_pump (mM/ms)
}
    
BREAKPOINT {
    SOLVE state METHOD cnexp
}


DERIVATIVE state { 
    
    : force concentration to stay above cainf by only pumping if larger
    drive_channel = -drive*ical/(2*FARADAY*depth)
    drive_pump    = -kt*(cali-cainf)/(cali+kd)
    
    if (drive_channel <= 0.) { drive_channel = 0. }
    
    cali' = drive_channel + pump*drive_pump + (cainf-cali)/taur
}

[1uc]

(deleted.)

[wthun]

(deleted.)

I see that the comment is deleted, but I did try using ASAN. Is it this line the issue refers to

nrn/src/nrnoc/init.cpp

Line 784 in c918f6b

return s->subtype + i * 1000;

?

Not sure how to test if this becomes a problem since it would just cause a normal overwrite?

[1uc]

The comment was incorrect and therefore deleted. I know 1000 is a special number and things are separated by "enough", where enough means 1000 (somehow somewhere).

What I would do is try with more than 1000 species, just to be sure. Then it would also be good to integrate the test into the automated test suite. With some luck you can put the MOD files and python script in the directory test/rxd and they should be picked up:

nrn/test/CMakeLists.txt

Lines 321 to 324 in c918f6b

	nrn_add_test_group(
	NAME rxdmod_tests
	MODFILE_PATTERNS test/rxd/*.mod
	SCRIPT_PATTERNS test/rxd/*.py test/rxd/*/*.py test/rxd/3d/*.asc test/rxd/testdata/**/*.dat)

[1uc]

@wthun I noticed that the PR if from wthun:master into master. This means that when you merge other, unrelated changes into wthun:master they appear as part of this PR.

The original PR seems like a good idea. NRN moves a little slower than small software projects; and there's vacations going on. It'll need review from others as well, and I was hoping one could set up some additional testing (possibly by someone from our side).

To keep unrelated changes separate it would be helpful if you could create a PR from a dedicated branch into our master. Naturally, you can merge those branches into your copy of NRN which you can use if/until NRN has the features you need for your research.

Edit: the new removed changes seem interesting too. (In a separate PR.)

[wthun]

@wthun I noticed that the PR if from wthun:master into master. This means that when you merge other, unrelated changes into wthun:master they appear as part of this PR.

The original PR seems like a good idea. NRN moves a little slower than small software projects; and there's vacations going on. It'll need review from others as well, and I was hoping one could set up some additional testing (possibly by someone from our side).

To keep unrelated changes separate it would be helpful if you could create a PR from a dedicated branch into our master. Naturally, you can merge those branches into your copy of NRN which you can use if/until NRN has the features you need for your research.

Edit: the new removed changes seem interesting too. (In a separate PR.)

Thanks! Yes, I mixed up my directories, fixing this.

[nrnhines]

I believe this change is sound but needs a test for interactions with some other parts of the code when an ion with mechanism index > 1000 is actually used. I'll add this test in a branch off this PR.

[1uc]

@wthun thank you for a very nice contribution.

In #3081, @nrnhines developed a test that shows that there will be an issue with the magic 1000. We understand that there's interesting science that requires more Ions and would like to work towards enabling that usecase. However, it feels non-trivial.

I'd like to propose the following:

1. We merge a version of your contribution with a small max_ions, since it's a nice refactoring.
2. (optional) We (e.g. me) implement CMake code to allow picking the value NRN_INTERNAL_MAX_IONS.
3. We continue working towards fully supporting the usecase.

@wthun is there a value of max_ions that's way less than 1000, e.g. 128 or 256, and large enough to allow you to progress? Second question: have you been able to use the modified version of NEURON with a large number of chemical species, i.e. >1000? (If it's also failing for you, then there's no value in adding a CMake option.)

[1uc]

Very nice refactoring!

[1uc]

In order to format NEURON code run:

./external/coding-conventions/bin/format

from the root directory. This will download the required version of all formatters in a venv called .bbp-project; and then run them on everything.

[codecov]

Codecov Report

Attention: Patch coverage is 81.81818% with 2 lines in your changes missing coverage. Please review.

Project coverage is 67.39%. Comparing base (e462814) to head (dc11e93).
Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
src/nrnoc/eion.cpp	81.81%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #3055      +/-   ##
==========================================
- Coverage   67.40%   67.39%   -0.02%     
==========================================
  Files         573      573              
  Lines      104927   104927              
==========================================
- Hits        70729    70716      -13     
- Misses      34198    34211      +13     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[nrnhines]

We merge a version of your contribution with a small max_ions, since it's a nice refactoring.

I agree with @1uc but speculate that in the end the elimination of the magic number 1000 will be straightforward (but tedious) with a natural limit on number of mechanisms of maxint/2 just by changing the coding of semantics for the two ion cases.

[nrnhines]

Will deal with the magic number 1000 in a later PR along with additional tests.

[wthun]

@1uc Sounds good! I ran the formatting command and set the number of ions to 256. Our (@Hjorthmedh) use case had >60 ions/species (per compartment) for the D1 and M4 signal transduction cascades in striatal projection neurons, but this number could grow quickly when more receptors are included or different cascades are used for different cells, so maybe 256 would be a good margin?

[sonarcloud]

Quality Gate passed

Issues
1 New issue
0 Accepted issues

Measures
0 Security Hotspots
0.0% Coverage on New Code
0.0% Duplication on New Code

See analysis details on SonarCloud

[1uc]

Yes, 256 is unlikely to cause issues. If it allows you to progress, then that buys us some time to fix the other issues; and longer term will be able to allow many more ions.

[nrnhines]

@wthun

... Our (@Hjorthmedh) use case had >60 ions/species (per compartment) for the D1 and M4 signal transduction cascades in striatal projection neurons, but this number could grow quickly when more receptors are included or different cascades are used for different cells,...

I think I have more or less finished up the PR's that allow an unlimited number of types of ions that can exist within a compartment (#3081, #3089, #3097). However I am curious about how the ions are being used within your model. Are they used by mod files or exclusively used within in the RxD framework? (If more convenient, I'm happy to carry on this discussion via email ([email protected]). ) Behind my question is a speculation about whether RxD should/could provide some features that would avoid the fairly elaborate apparatus behind full ion mechanism behavior provided by the mod file USEION statement.

[wthun]

@nrnhines
Thanks! That was quick.

We've initially represented signal transduction cascades within RxD, and their downstream effects by having membrane mechanisms read the concentration of e.g. PKA. Synaptic release to trigger the cascades also use .mod files along with include_flux from RxD to couple them.

However, we've also considered representing the signal transduction cascades in mod files due to run time requirements. In our benchmarks, RxD brings a 10-20x slow down compared to having no cascade, while the same model represented as ODE:s in a mod file brings only a slow down of 1.45x. The latter didn't write the ions to the membrane, we only tested the run time required to compute the reaction by itself. The ODE representation also seem to scale linearly with the number of ODE:s, while the RxD representation seem to scale cubically with the number of species. Solving for the Jacobian essentially accounts for the full slow down in benchmarks.

We've speculated that maybe we could generate nmodl code from our RxD specifixation (e.g. using the nmodl framework) but still harmonize with the rest of the RxD features. But maybe that's fairly tangential to this question?