From 85486cfa122c806ec08c0a33dc58e5694720972b Mon Sep 17 00:00:00 2001 From: Jannik Luboeinski Date: Tue, 15 Aug 2023 10:54:54 +0200 Subject: [PATCH 01/80] Tests for particle diffusion added (for different morphologies; also pointing to the fact that the diffusion across segments of different radius is still erroneous; see https://github.com/arbor-sim/arbor/issues/2145) --- .gitignore | 3 + python/test/fixtures.py | 12 + python/test/readme.md | 4 +- python/test/unit/test_diffusion.py | 238 ++++++++++++++++++ test/unit/diffusion/neuron_with_diffusion.mod | 23 ++ .../unit/diffusion/synapse_with_diffusion.mod | 24 ++ 6 files changed, 302 insertions(+), 2 deletions(-) create mode 100644 python/test/unit/test_diffusion.py create mode 100644 test/unit/diffusion/neuron_with_diffusion.mod create mode 100644 test/unit/diffusion/synapse_with_diffusion.mod diff --git a/.gitignore b/.gitignore index 22631b86a3..9178ba0f25 100644 --- a/.gitignore +++ b/.gitignore @@ -105,3 +105,6 @@ _skbuild # generated by test scripts results + +# environment setting script +set_arbor_env* diff --git a/python/test/fixtures.py b/python/test/fixtures.py index 365d69fa8c..2d3e82f13e 100644 --- a/python/test/fixtures.py +++ b/python/test/fixtures.py @@ -172,6 +172,18 @@ def dummy_catalogue(repo_path): cat_path = _build_cat("dummy", path) return arbor.load_catalogue(str(cat_path)) +@_singleton_fixture +@repo_path() +def diffusion_catalogue(repo_path): + """ + Fixture that returns an `arbor.catalogue` + which contains mechanisms `neuron_with_diffusion` + and `synapse_with_diffusion`. + """ + path = repo_path / "test" / "unit" / "diffusion" + cat_path = _build_cat("diffusion", path) + return arbor.load_catalogue(str(cat_path)) + @_fixture class empty_recipe(arbor.recipe): diff --git a/python/test/readme.md b/python/test/readme.md index 515366df87..613e455a82 100644 --- a/python/test/readme.md +++ b/python/test/readme.md @@ -39,6 +39,6 @@ In subfolders `unit`/`unit_distributed`: ## Naming convention -- modules: `test_xxxs` (ending with `s` since module can consist of multiple classes) -- class(es): `TestXxxs` (ending with `s` since class can consist of multiple test functions) +- modules: `test_xxxs` (if suitable, ending with `s` since module can consist of multiple classes) +- class(es): `TestXxxs` (if suitable, ending with `s` since class can consist of multiple test functions) - functions: `test_yyy` diff --git a/python/test/unit/test_diffusion.py b/python/test/unit/test_diffusion.py new file mode 100644 index 0000000000..78ff7809f8 --- /dev/null +++ b/python/test/unit/test_diffusion.py @@ -0,0 +1,238 @@ +# -*- coding: utf-8 -*- + +import unittest +import arbor as A +import numpy as np +from .. import fixtures + +""" +Tests for the concentration and amount of diffusive particles across time and morphology. +Three different morphological structures are considered: 1 segment ("soma only"), 2 segments +("soma with dendrite"), and 3 segments ("soma with two dendrites"). + +NOTE: Internally, Arbor only knows concentrations. Thus, particle amounts have to be computed + from concentrations by integrating over the volume of the morphology. The total amount + of particles should be conserved unless there is deliberate injection or removal of + particles. +""" + +# --------------------------------------------------------------------------------------- +# recipe class +class recipe(A.recipe): + def __init__(self, cat, cell, probes, inc_0, inc_1, dec_0): + A.recipe.__init__(self) + self.the_cell = cell + self.the_probes = probes + self.the_props = A.neuron_cable_properties() + self.the_props.catalogue = cat # use the provided catalogue of diffusion mechanisms + self.the_props.set_ion("s", 1, 0, 0, 0) # use diffusive particles "s" + self.inc_0 = inc_0 # increase in particle amount at 0.1 s (in 1e-18 mol) + self.inc_1 = inc_1 # increase in particle amount at 0.5 s (in 1e-18 mol) + self.dec_0 = dec_0 # decrease in particle amount at 1.5 s (in 1e-18 mol) + + def num_cells(self): + return 1 + + def cell_kind(self, gid): + return A.cell_kind.cable + + def cell_description(self, gid): + return self.the_cell + + def probes(self, gid): + return self.the_probes + + def global_properties(self, kind): + return self.the_props + + def event_generators(self, gid): + g = [A.event_generator("syn_exc_A", self.inc_0, A.explicit_schedule([0.1])), + A.event_generator("syn_exc_B", self.inc_1, A.explicit_schedule([0.5])), + A.event_generator("syn_inh", -self.dec_0, A.explicit_schedule([1.5]))] + return g + +# --------------------------------------------------------------------------------------- +# test class +class TestDiffusion(unittest.TestCase): + + # Constructor (overridden) + def __init__(self, args): + super(TestDiffusion, self).__init__(args) + + self.runtime = 3.00 # runtime of the whole simulation in ms + self.dt = 0.01 # duration of one timestep in ms + self.dev = 0.01 # accepted relative deviation for `assertAlmostEqual` + + # Method to run an Arbor simulation with diffusion across different segments + def simulate_diffusion(self, cat, _num_segs, _num_cvs_per_seg, _length, _r_1, _r_2 = 0., _r_3 = 0.): + + # --------------------------------------------------------------------------------------- + # set the main parameters and calculate geometrical measures + num_segs = _num_segs # number of segments + num_cvs_per_seg = _num_cvs_per_seg # number of CVs per segment + + length = _length # length of the whole setup (in case of 1 or 2 segments, one branch) in µm + radius_1 = _r_1 # radius of the first segment in µm + if num_segs > 1: + radius_2 = _r_2 # radius of the second segment in µm + else: + radius_2 = 0 + if num_segs > 2: + radius_3 = _r_3 # radius of the third segment in µm + else: + radius_3 = 0 + + length_per_seg = length / num_segs # axial length of a segment in µm + area_tot = 2 * np.pi * (radius_1 + radius_2 + radius_3) * length_per_seg # surface area of the whole setup in µm^2 (excluding the circle-shaped ends, since Arbor does not consider current flux there) + area_per_cv = area_tot / (num_segs*num_cvs_per_seg) # surface area of one cylindrical CV in µm^2 (excluding the circle-shaped ends, since Arbor does not consider current flux there) + volume_tot = np.pi * (radius_1**2 + radius_2**2 + radius_3**2) * length_per_seg # volume of the whole setup in µm^3 + volume_per_cv = volume_tot / (num_segs*num_cvs_per_seg) # volume of one cylindrical CV in µm^3 + + inc_0 = 600 # first increase in particle amount (in 1e-18 mol) + inc_1 = 1200 # second increase in particle amount (in 1e-18 mol) + dec_0 = 1400 # decrease in particle amount (in 1e-18 mol) + diffusivity = 1 # diffusivity (in m^2/s) + + # --------------------------------------------------------------------------------------- + # set up the morphology + tree = A.segment_tree() + if num_segs == 1: + _ = tree.append(A.mnpos, + A.mpoint(-length/2, 0, 0, radius_1), + A.mpoint(+length/2, 0, 0, radius_1), + tag=0) + + labels = A.label_dict({"soma-region" : "(tag 0)", + "soma-center" : "(on-components 0.5 (region \"soma-region\"))", + "soma-end" : "(on-components 1.0 (region \"soma-region\"))"}) + elif num_segs == 2: + s = tree.append(A.mnpos, + A.mpoint(-length/2, 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/2, 0, 0, radius_2), + tag=1) + + labels = A.label_dict({"soma-region" : "(tag 0)", + "dendriteA-region" : "(tag 1)", + "soma-center" : "(on-components 0.5 (region \"soma-region\"))", + "soma-end" : "(on-components 1.0 (region \"soma-region\"))", + "dendriteA-center" : "(on-components 0.5 (region \"dendriteA-region\"))"}) + elif num_segs == 3: + + s = tree.append(A.mnpos, + A.mpoint(-1/3*length, 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(+1/3*length, 0, 0, radius_2), + tag=1) + _ = tree.append(s, + A.mpoint(-1/3*length, 0, 0, radius_3), + A.mpoint(-2/3*length, 0, 0, radius_3), + tag=2) + + + labels = A.label_dict({"soma-region" : "(tag 0)", + "dendriteA-region" : "(tag 1)", + "dendriteB-region" : "(tag 2)", + "soma-center" : "(on-components 0.5 (region \"soma-region\"))", + "soma-end" : "(on-components 1.0 (region \"soma-region\"))", + "dendriteA-center" : "(on-components 0.5 (region \"dendriteA-region\"))", + "dendriteB-center" : "(on-components 0.5 (region \"dendriteB-region\"))"}) + else: + raise ValueError(f"Specified number of segments not supported.") + morph = A.morphology(tree) + + # --------------------------------------------------------------------------------------- + # decorate the morphology with mechanisms + dec = A.decor() + if num_segs < 3: + dec.discretization(A.cv_policy(f'(fixed-per-branch {num_segs*num_cvs_per_seg} (branch 0))')) + elif num_segs == 3: + dec.discretization(A.cv_policy(f'(replace (fixed-per-branch {num_cvs_per_seg} (branch 0)) ' + \ + f'(fixed-per-branch {num_cvs_per_seg} (branch 1)) ' + \ + f'(fixed-per-branch {num_cvs_per_seg} (branch 2)))')) + dec.set_ion("s", int_con=0.0, diff=diffusivity) + if num_segs == 1: + dec.place('"soma-end"', A.synapse("synapse_with_diffusion"), "syn_exc_A") + dec.place('"soma-end"', A.synapse("synapse_with_diffusion"), "syn_exc_B") + dec.place('"soma-end"', A.synapse("synapse_with_diffusion"), "syn_inh") + elif num_segs == 2: + dec.place('"dendriteA-center"', A.synapse("synapse_with_diffusion"), "syn_exc_A") + dec.place('"soma-end"', A.synapse("synapse_with_diffusion"), "syn_exc_B") + dec.place('"soma-end"', A.synapse("synapse_with_diffusion"), "syn_inh") + elif num_segs == 3: + dec.place('"dendriteA-center"', A.synapse("synapse_with_diffusion"), "syn_exc_A") + dec.place('"dendriteB-center"', A.synapse("synapse_with_diffusion"), "syn_exc_B") + dec.place('"soma-end"', A.synapse("synapse_with_diffusion"), "syn_inh") + dec.paint('(all)', A.density("neuron_with_diffusion")) + + # --------------------------------------------------------------------------------------- + # set probes + prb = [A.cable_probe_ion_diff_concentration('"soma-center"', "s"), + A.cable_probe_density_state('"soma-center"', "neuron_with_diffusion", "sV"), + A.cable_probe_density_state_cell("neuron_with_diffusion", "sV")] + + # --------------------------------------------------------------------------------------- + # prepare the simulation + cel = A.cable_cell(tree, dec, labels) + rec = recipe(cat, cel, prb, inc_0, inc_1, dec_0) + sim = A.simulation(rec) + + # --------------------------------------------------------------------------------------- + # set handles + hdl_s = sim.sample((0, 0), A.regular_schedule(self.dt)) # s at "soma-center" + hdl_sV = sim.sample((0, 1), A.regular_schedule(self.dt)) # sV at "soma-center" + hdl_sV_all = sim.sample((0, 2), A.regular_schedule(self.dt)) # sV (cell-wide array) + + # --------------------------------------------------------------------------------------- + # run the simulation + sim.run(dt=self.dt, tfinal=self.runtime) + + # --------------------------------------------------------------------------------------- + # retrieve data and do the testing + data_s = sim.samples(hdl_s)[0][0] + times = data_s[:,0] + data_sV = sim.samples(hdl_sV)[0][0] + tmp_data = sim.samples(hdl_sV_all)[0][0] + data_sV_total = np.zeros_like(tmp_data[:, 0]) + num_cvs = len(tmp_data[0,:])-1 + for i in range(len(tmp_data[0,:])-1): # compute the total amount of particles by summing over all CVs of the whole neuron + data_sV_total += tmp_data[:, i+1] + + s_lim_expected = (inc_0+inc_1-dec_0) / volume_tot # total particle amount of s divided by total volume + s_max_expected = (inc_0+inc_1) / volume_tot # total particle amount of s divided by total volume + + if num_segs < 3: + self.assertEqual(morph.num_branches, 1) # expected number of branches: 1 + else: + self.assertEqual(morph.num_branches, 3) # expected number of branches: 3 + self.assertEqual(num_cvs, num_segs*num_cvs_per_seg) # expected total number of CVs + self.assertAlmostEqual(data_sV[-1,1]/volume_per_cv, s_lim_expected, delta=self.dev*s_lim_expected) # lim_{t->inf}(s) [estimated] + self.assertAlmostEqual(data_s[-1,1], s_lim_expected, delta=self.dev*s_lim_expected) # lim_{t->inf}(s) [direct] + self.assertAlmostEqual(np.max(data_s[:,1]), s_max_expected, delta=self.dev*s_max_expected) # max_{t}(s) [direct] + self.assertAlmostEqual(data_sV[-1,1]*num_segs*num_cvs_per_seg, inc_0+inc_1-dec_0, delta=self.dev*(inc_0+inc_1-dec_0)) # lim_{t->inf}(s⋅V) [estimated] + self.assertAlmostEqual(data_sV_total[-1], inc_0+inc_1-dec_0, delta=self.dev*(inc_0+inc_1-dec_0)) # lim_{t->inf}(s⋅V) [direct] + self.assertAlmostEqual(np.max(data_sV_total), inc_0+inc_1, delta=self.dev*(inc_0+inc_1)) # max_{t}(s⋅V) [direct] + + # Test: simulations with equal radii + @fixtures.diffusion_catalogue() + def test_diffusion_equal_radii(self, diffusion_catalogue): + + self.simulate_diffusion(diffusion_catalogue, 1, 600, 10, 4) # 1 segment with radius 4 µm + self.simulate_diffusion(diffusion_catalogue, 2, 300, 10, 4, 4) # 2 segments with radius 4 µm + self.simulate_diffusion(diffusion_catalogue, 3, 200, 10, 4, 4, 4) # 3 segments with radius 4 µm + + ''' TODO: not succeeding as of Arbor v0.9.0: + # Test: simulations with different radii + @fixtures.diffusion_catalogue() + def test_diffusion_different_radii(self, diffusion_catalogue): + + self.simulate_diffusion(diffusion_catalogue, 2, 300, 10, 4, 6) # 2 segments with radius 4 µm and 6 µm + self.simulate_diffusion(diffusion_catalogue, 3, 200, 10, 4, 6, 6) # 3 segments with radius 4 µm and 6 µm + ''' + diff --git a/test/unit/diffusion/neuron_with_diffusion.mod b/test/unit/diffusion/neuron_with_diffusion.mod new file mode 100644 index 0000000000..41fbe9a3fd --- /dev/null +++ b/test/unit/diffusion/neuron_with_diffusion.mod @@ -0,0 +1,23 @@ +NEURON { + SUFFIX neuron_with_diffusion + USEION s READ sd +} + +PARAMETER { + area : surface area of the CV (in µm^2, internal variable) + diam : CV diameter (in µm, internal variable) +} + +ASSIGNED { + volume : volume of the CV (conversion factor between concentration and particle amount, in µm^3) + sV : particle amount in the CV (in 1e-18 mol) +} + +INITIAL { + volume = area*diam/4 : = area*r/2 = 2*pi*r*h*r/2 = pi*r^2*h + sV = sd * volume +} + +BREAKPOINT { + sV = sd * volume : read and normalize particle amount +} diff --git a/test/unit/diffusion/synapse_with_diffusion.mod b/test/unit/diffusion/synapse_with_diffusion.mod new file mode 100644 index 0000000000..946f164018 --- /dev/null +++ b/test/unit/diffusion/synapse_with_diffusion.mod @@ -0,0 +1,24 @@ +NEURON { + POINT_PROCESS synapse_with_diffusion + USEION s WRITE sd +} + +PARAMETER { + area : surface area of the CV (in µm^2, internal variable) + diam : CV diameter (in µm, internal variable) +} + +ASSIGNED { + volume : volume of the CV (conversion factor between concentration and particle amount, in µm^3) +} + +INITIAL { + volume = area*diam/4 : = area*r/2 = 2*pi*r*h*r/2 = pi*r^2*h +} + +BREAKPOINT { +} + +NET_RECEIVE(weight) { + sd = sd + weight * area / volume / 1000 +} From 6062dafb98d8c2db6337a96bf1d155740d5a8344 Mon Sep 17 00:00:00 2001 From: Jannik Luboeinski Date: Tue, 15 Aug 2023 11:26:08 +0200 Subject: [PATCH 02/80] Black style formatting applied --- python/test/fixtures.py | 1 + python/test/unit/test_diffusion.py | 303 ++++++++++++++++++----------- 2 files changed, 195 insertions(+), 109 deletions(-) diff --git a/python/test/fixtures.py b/python/test/fixtures.py index 2d3e82f13e..b978481c83 100644 --- a/python/test/fixtures.py +++ b/python/test/fixtures.py @@ -172,6 +172,7 @@ def dummy_catalogue(repo_path): cat_path = _build_cat("dummy", path) return arbor.load_catalogue(str(cat_path)) + @_singleton_fixture @repo_path() def diffusion_catalogue(repo_path): diff --git a/python/test/unit/test_diffusion.py b/python/test/unit/test_diffusion.py index 78ff7809f8..3e1ea2620a 100644 --- a/python/test/unit/test_diffusion.py +++ b/python/test/unit/test_diffusion.py @@ -24,11 +24,13 @@ def __init__(self, cat, cell, probes, inc_0, inc_1, dec_0): self.the_cell = cell self.the_probes = probes self.the_props = A.neuron_cable_properties() - self.the_props.catalogue = cat # use the provided catalogue of diffusion mechanisms - self.the_props.set_ion("s", 1, 0, 0, 0) # use diffusive particles "s" - self.inc_0 = inc_0 # increase in particle amount at 0.1 s (in 1e-18 mol) - self.inc_1 = inc_1 # increase in particle amount at 0.5 s (in 1e-18 mol) - self.dec_0 = dec_0 # decrease in particle amount at 1.5 s (in 1e-18 mol) + self.the_props.catalogue = ( + cat # use the provided catalogue of diffusion mechanisms + ) + self.the_props.set_ion("s", 1, 0, 0, 0) # use diffusive particles "s" + self.inc_0 = inc_0 # increase in particle amount at 0.1 s (in 1e-18 mol) + self.inc_1 = inc_1 # increase in particle amount at 0.5 s (in 1e-18 mol) + self.dec_0 = dec_0 # decrease in particle amount at 1.5 s (in 1e-18 mol) def num_cells(self): return 1 @@ -46,11 +48,14 @@ def global_properties(self, kind): return self.the_props def event_generators(self, gid): - g = [A.event_generator("syn_exc_A", self.inc_0, A.explicit_schedule([0.1])), - A.event_generator("syn_exc_B", self.inc_1, A.explicit_schedule([0.5])), - A.event_generator("syn_inh", -self.dec_0, A.explicit_schedule([1.5]))] + g = [ + A.event_generator("syn_exc_A", self.inc_0, A.explicit_schedule([0.1])), + A.event_generator("syn_exc_B", self.inc_1, A.explicit_schedule([0.5])), + A.event_generator("syn_inh", -self.dec_0, A.explicit_schedule([1.5])), + ] return g + # --------------------------------------------------------------------------------------- # test class class TestDiffusion(unittest.TestCase): @@ -59,90 +64,123 @@ class TestDiffusion(unittest.TestCase): def __init__(self, args): super(TestDiffusion, self).__init__(args) - self.runtime = 3.00 # runtime of the whole simulation in ms - self.dt = 0.01 # duration of one timestep in ms - self.dev = 0.01 # accepted relative deviation for `assertAlmostEqual` + self.runtime = 3.00 # runtime of the whole simulation in ms + self.dt = 0.01 # duration of one timestep in ms + self.dev = 0.01 # accepted relative deviation for `assertAlmostEqual` # Method to run an Arbor simulation with diffusion across different segments - def simulate_diffusion(self, cat, _num_segs, _num_cvs_per_seg, _length, _r_1, _r_2 = 0., _r_3 = 0.): + def simulate_diffusion( + self, cat, _num_segs, _num_cvs_per_seg, _length, _r_1, _r_2=0.0, _r_3=0.0 + ): # --------------------------------------------------------------------------------------- # set the main parameters and calculate geometrical measures - num_segs = _num_segs # number of segments - num_cvs_per_seg = _num_cvs_per_seg # number of CVs per segment + num_segs = _num_segs # number of segments + num_cvs_per_seg = _num_cvs_per_seg # number of CVs per segment - length = _length # length of the whole setup (in case of 1 or 2 segments, one branch) in µm - radius_1 = _r_1 # radius of the first segment in µm + length = _length # length of the whole setup (in case of 1 or 2 segments, one branch) in µm + radius_1 = _r_1 # radius of the first segment in µm if num_segs > 1: - radius_2 = _r_2 # radius of the second segment in µm + radius_2 = _r_2 # radius of the second segment in µm else: radius_2 = 0 if num_segs > 2: - radius_3 = _r_3 # radius of the third segment in µm + radius_3 = _r_3 # radius of the third segment in µm else: radius_3 = 0 - length_per_seg = length / num_segs # axial length of a segment in µm - area_tot = 2 * np.pi * (radius_1 + radius_2 + radius_3) * length_per_seg # surface area of the whole setup in µm^2 (excluding the circle-shaped ends, since Arbor does not consider current flux there) - area_per_cv = area_tot / (num_segs*num_cvs_per_seg) # surface area of one cylindrical CV in µm^2 (excluding the circle-shaped ends, since Arbor does not consider current flux there) - volume_tot = np.pi * (radius_1**2 + radius_2**2 + radius_3**2) * length_per_seg # volume of the whole setup in µm^3 - volume_per_cv = volume_tot / (num_segs*num_cvs_per_seg) # volume of one cylindrical CV in µm^3 - - inc_0 = 600 # first increase in particle amount (in 1e-18 mol) - inc_1 = 1200 # second increase in particle amount (in 1e-18 mol) - dec_0 = 1400 # decrease in particle amount (in 1e-18 mol) - diffusivity = 1 # diffusivity (in m^2/s) + length_per_seg = length / num_segs # axial length of a segment in µm + area_tot = ( + 2 * np.pi * (radius_1 + radius_2 + radius_3) * length_per_seg + ) # surface area of the whole setup in µm^2 (excluding the circle-shaped ends, since Arbor does not consider current flux there) + area_per_cv = area_tot / ( + num_segs * num_cvs_per_seg + ) # surface area of one cylindrical CV in µm^2 (excluding the circle-shaped ends, since Arbor does not consider current flux there) + volume_tot = ( + np.pi * (radius_1 ** 2 + radius_2 ** 2 + radius_3 ** 2) * length_per_seg + ) # volume of the whole setup in µm^3 + volume_per_cv = volume_tot / ( + num_segs * num_cvs_per_seg + ) # volume of one cylindrical CV in µm^3 + + inc_0 = 600 # first increase in particle amount (in 1e-18 mol) + inc_1 = 1200 # second increase in particle amount (in 1e-18 mol) + dec_0 = 1400 # decrease in particle amount (in 1e-18 mol) + diffusivity = 1 # diffusivity (in m^2/s) # --------------------------------------------------------------------------------------- # set up the morphology tree = A.segment_tree() if num_segs == 1: - _ = tree.append(A.mnpos, - A.mpoint(-length/2, 0, 0, radius_1), - A.mpoint(+length/2, 0, 0, radius_1), - tag=0) - - labels = A.label_dict({"soma-region" : "(tag 0)", - "soma-center" : "(on-components 0.5 (region \"soma-region\"))", - "soma-end" : "(on-components 1.0 (region \"soma-region\"))"}) + _ = tree.append( + A.mnpos, + A.mpoint(-length / 2, 0, 0, radius_1), + A.mpoint(+length / 2, 0, 0, radius_1), + tag=0, + ) + + labels = A.label_dict( + { + "soma-region": "(tag 0)", + "soma-center": '(on-components 0.5 (region "soma-region"))', + "soma-end": '(on-components 1.0 (region "soma-region"))', + } + ) elif num_segs == 2: - s = tree.append(A.mnpos, - A.mpoint(-length/2, 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/2, 0, 0, radius_2), - tag=1) - - labels = A.label_dict({"soma-region" : "(tag 0)", - "dendriteA-region" : "(tag 1)", - "soma-center" : "(on-components 0.5 (region \"soma-region\"))", - "soma-end" : "(on-components 1.0 (region \"soma-region\"))", - "dendriteA-center" : "(on-components 0.5 (region \"dendriteA-region\"))"}) + s = tree.append( + A.mnpos, + A.mpoint(-length / 2, 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 / 2, 0, 0, radius_2), + tag=1, + ) + + labels = A.label_dict( + { + "soma-region": "(tag 0)", + "dendriteA-region": "(tag 1)", + "soma-center": '(on-components 0.5 (region "soma-region"))', + "soma-end": '(on-components 1.0 (region "soma-region"))', + "dendriteA-center": '(on-components 0.5 (region "dendriteA-region"))', + } + ) elif num_segs == 3: - - s = tree.append(A.mnpos, - A.mpoint(-1/3*length, 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(+1/3*length, 0, 0, radius_2), - tag=1) - _ = tree.append(s, - A.mpoint(-1/3*length, 0, 0, radius_3), - A.mpoint(-2/3*length, 0, 0, radius_3), - tag=2) - - - labels = A.label_dict({"soma-region" : "(tag 0)", - "dendriteA-region" : "(tag 1)", - "dendriteB-region" : "(tag 2)", - "soma-center" : "(on-components 0.5 (region \"soma-region\"))", - "soma-end" : "(on-components 1.0 (region \"soma-region\"))", - "dendriteA-center" : "(on-components 0.5 (region \"dendriteA-region\"))", - "dendriteB-center" : "(on-components 0.5 (region \"dendriteB-region\"))"}) + + s = tree.append( + A.mnpos, + A.mpoint(-1 / 3 * length, 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(+1 / 3 * length, 0, 0, radius_2), + tag=1, + ) + _ = tree.append( + s, + A.mpoint(-1 / 3 * length, 0, 0, radius_3), + A.mpoint(-2 / 3 * length, 0, 0, radius_3), + tag=2, + ) + + labels = A.label_dict( + { + "soma-region": "(tag 0)", + "dendriteA-region": "(tag 1)", + "dendriteB-region": "(tag 2)", + "soma-center": '(on-components 0.5 (region "soma-region"))', + "soma-end": '(on-components 1.0 (region "soma-region"))', + "dendriteA-center": '(on-components 0.5 (region "dendriteA-region"))', + "dendriteB-center": '(on-components 0.5 (region "dendriteB-region"))', + } + ) else: raise ValueError(f"Specified number of segments not supported.") morph = A.morphology(tree) @@ -151,31 +189,45 @@ def simulate_diffusion(self, cat, _num_segs, _num_cvs_per_seg, _length, _r_1, _r # decorate the morphology with mechanisms dec = A.decor() if num_segs < 3: - dec.discretization(A.cv_policy(f'(fixed-per-branch {num_segs*num_cvs_per_seg} (branch 0))')) + dec.discretization( + A.cv_policy(f"(fixed-per-branch {num_segs*num_cvs_per_seg} (branch 0))") + ) elif num_segs == 3: - dec.discretization(A.cv_policy(f'(replace (fixed-per-branch {num_cvs_per_seg} (branch 0)) ' + \ - f'(fixed-per-branch {num_cvs_per_seg} (branch 1)) ' + \ - f'(fixed-per-branch {num_cvs_per_seg} (branch 2)))')) + dec.discretization( + A.cv_policy( + f"(replace (fixed-per-branch {num_cvs_per_seg} (branch 0)) " + + f"(fixed-per-branch {num_cvs_per_seg} (branch 1)) " + + f"(fixed-per-branch {num_cvs_per_seg} (branch 2)))" + ) + ) dec.set_ion("s", int_con=0.0, diff=diffusivity) if num_segs == 1: dec.place('"soma-end"', A.synapse("synapse_with_diffusion"), "syn_exc_A") dec.place('"soma-end"', A.synapse("synapse_with_diffusion"), "syn_exc_B") dec.place('"soma-end"', A.synapse("synapse_with_diffusion"), "syn_inh") elif num_segs == 2: - dec.place('"dendriteA-center"', A.synapse("synapse_with_diffusion"), "syn_exc_A") + dec.place( + '"dendriteA-center"', A.synapse("synapse_with_diffusion"), "syn_exc_A" + ) dec.place('"soma-end"', A.synapse("synapse_with_diffusion"), "syn_exc_B") dec.place('"soma-end"', A.synapse("synapse_with_diffusion"), "syn_inh") elif num_segs == 3: - dec.place('"dendriteA-center"', A.synapse("synapse_with_diffusion"), "syn_exc_A") - dec.place('"dendriteB-center"', A.synapse("synapse_with_diffusion"), "syn_exc_B") + dec.place( + '"dendriteA-center"', A.synapse("synapse_with_diffusion"), "syn_exc_A" + ) + dec.place( + '"dendriteB-center"', A.synapse("synapse_with_diffusion"), "syn_exc_B" + ) dec.place('"soma-end"', A.synapse("synapse_with_diffusion"), "syn_inh") - dec.paint('(all)', A.density("neuron_with_diffusion")) + dec.paint("(all)", A.density("neuron_with_diffusion")) # --------------------------------------------------------------------------------------- # set probes - prb = [A.cable_probe_ion_diff_concentration('"soma-center"', "s"), - A.cable_probe_density_state('"soma-center"', "neuron_with_diffusion", "sV"), - A.cable_probe_density_state_cell("neuron_with_diffusion", "sV")] + prb = [ + A.cable_probe_ion_diff_concentration('"soma-center"', "s"), + A.cable_probe_density_state('"soma-center"', "neuron_with_diffusion", "sV"), + A.cable_probe_density_state_cell("neuron_with_diffusion", "sV"), + ] # --------------------------------------------------------------------------------------- # prepare the simulation @@ -185,9 +237,11 @@ def simulate_diffusion(self, cat, _num_segs, _num_cvs_per_seg, _length, _r_1, _r # --------------------------------------------------------------------------------------- # set handles - hdl_s = sim.sample((0, 0), A.regular_schedule(self.dt)) # s at "soma-center" - hdl_sV = sim.sample((0, 1), A.regular_schedule(self.dt)) # sV at "soma-center" - hdl_sV_all = sim.sample((0, 2), A.regular_schedule(self.dt)) # sV (cell-wide array) + hdl_s = sim.sample((0, 0), A.regular_schedule(self.dt)) # s at "soma-center" + hdl_sV = sim.sample((0, 1), A.regular_schedule(self.dt)) # sV at "soma-center" + hdl_sV_all = sim.sample( + (0, 2), A.regular_schedule(self.dt) + ) # sV (cell-wide array) # --------------------------------------------------------------------------------------- # run the simulation @@ -196,43 +250,74 @@ def simulate_diffusion(self, cat, _num_segs, _num_cvs_per_seg, _length, _r_1, _r # --------------------------------------------------------------------------------------- # retrieve data and do the testing data_s = sim.samples(hdl_s)[0][0] - times = data_s[:,0] + times = data_s[:, 0] data_sV = sim.samples(hdl_sV)[0][0] tmp_data = sim.samples(hdl_sV_all)[0][0] data_sV_total = np.zeros_like(tmp_data[:, 0]) - num_cvs = len(tmp_data[0,:])-1 - for i in range(len(tmp_data[0,:])-1): # compute the total amount of particles by summing over all CVs of the whole neuron - data_sV_total += tmp_data[:, i+1] - - s_lim_expected = (inc_0+inc_1-dec_0) / volume_tot # total particle amount of s divided by total volume - s_max_expected = (inc_0+inc_1) / volume_tot # total particle amount of s divided by total volume + num_cvs = len(tmp_data[0, :]) - 1 + for i in range( + len(tmp_data[0, :]) - 1 + ): # compute the total amount of particles by summing over all CVs of the whole neuron + data_sV_total += tmp_data[:, i + 1] + + s_lim_expected = ( + inc_0 + inc_1 - dec_0 + ) / volume_tot # total particle amount of s divided by total volume + s_max_expected = ( + inc_0 + inc_1 + ) / volume_tot # total particle amount of s divided by total volume if num_segs < 3: - self.assertEqual(morph.num_branches, 1) # expected number of branches: 1 + self.assertEqual(morph.num_branches, 1) # expected number of branches: 1 else: - self.assertEqual(morph.num_branches, 3) # expected number of branches: 3 - self.assertEqual(num_cvs, num_segs*num_cvs_per_seg) # expected total number of CVs - self.assertAlmostEqual(data_sV[-1,1]/volume_per_cv, s_lim_expected, delta=self.dev*s_lim_expected) # lim_{t->inf}(s) [estimated] - self.assertAlmostEqual(data_s[-1,1], s_lim_expected, delta=self.dev*s_lim_expected) # lim_{t->inf}(s) [direct] - self.assertAlmostEqual(np.max(data_s[:,1]), s_max_expected, delta=self.dev*s_max_expected) # max_{t}(s) [direct] - self.assertAlmostEqual(data_sV[-1,1]*num_segs*num_cvs_per_seg, inc_0+inc_1-dec_0, delta=self.dev*(inc_0+inc_1-dec_0)) # lim_{t->inf}(s⋅V) [estimated] - self.assertAlmostEqual(data_sV_total[-1], inc_0+inc_1-dec_0, delta=self.dev*(inc_0+inc_1-dec_0)) # lim_{t->inf}(s⋅V) [direct] - self.assertAlmostEqual(np.max(data_sV_total), inc_0+inc_1, delta=self.dev*(inc_0+inc_1)) # max_{t}(s⋅V) [direct] + self.assertEqual(morph.num_branches, 3) # expected number of branches: 3 + self.assertEqual( + num_cvs, num_segs * num_cvs_per_seg + ) # expected total number of CVs + self.assertAlmostEqual( + data_sV[-1, 1] / volume_per_cv, + s_lim_expected, + delta=self.dev * s_lim_expected, + ) # lim_{t->inf}(s) [estimated] + self.assertAlmostEqual( + data_s[-1, 1], s_lim_expected, delta=self.dev * s_lim_expected + ) # lim_{t->inf}(s) [direct] + self.assertAlmostEqual( + np.max(data_s[:, 1]), s_max_expected, delta=self.dev * s_max_expected + ) # max_{t}(s) [direct] + self.assertAlmostEqual( + data_sV[-1, 1] * num_segs * num_cvs_per_seg, + inc_0 + inc_1 - dec_0, + delta=self.dev * (inc_0 + inc_1 - dec_0), + ) # lim_{t->inf}(s⋅V) [estimated] + self.assertAlmostEqual( + data_sV_total[-1], + inc_0 + inc_1 - dec_0, + delta=self.dev * (inc_0 + inc_1 - dec_0), + ) # lim_{t->inf}(s⋅V) [direct] + self.assertAlmostEqual( + np.max(data_sV_total), inc_0 + inc_1, delta=self.dev * (inc_0 + inc_1) + ) # max_{t}(s⋅V) [direct] # Test: simulations with equal radii @fixtures.diffusion_catalogue() def test_diffusion_equal_radii(self, diffusion_catalogue): - self.simulate_diffusion(diffusion_catalogue, 1, 600, 10, 4) # 1 segment with radius 4 µm - self.simulate_diffusion(diffusion_catalogue, 2, 300, 10, 4, 4) # 2 segments with radius 4 µm - self.simulate_diffusion(diffusion_catalogue, 3, 200, 10, 4, 4, 4) # 3 segments with radius 4 µm - - ''' TODO: not succeeding as of Arbor v0.9.0: + self.simulate_diffusion( + diffusion_catalogue, 1, 600, 10, 4 + ) # 1 segment with radius 4 µm + self.simulate_diffusion( + diffusion_catalogue, 2, 300, 10, 4, 4 + ) # 2 segments with radius 4 µm + self.simulate_diffusion( + diffusion_catalogue, 3, 200, 10, 4, 4, 4 + ) # 3 segments with radius 4 µm + + """ TODO: not succeeding as of Arbor v0.9.0: # Test: simulations with different radii @fixtures.diffusion_catalogue() def test_diffusion_different_radii(self, diffusion_catalogue): self.simulate_diffusion(diffusion_catalogue, 2, 300, 10, 4, 6) # 2 segments with radius 4 µm and 6 µm self.simulate_diffusion(diffusion_catalogue, 3, 200, 10, 4, 6, 6) # 3 segments with radius 4 µm and 6 µm - ''' - + """ From 890f3151528d17c2fb89a1e1c46be9acd397f13f Mon Sep 17 00:00:00 2001 From: Jannik Luboeinski Date: Wed, 16 Aug 2023 12:02:01 +0200 Subject: [PATCH 03/80] Black style formatting (v23.7.0) applied --- python/test/unit/test_diffusion.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/python/test/unit/test_diffusion.py b/python/test/unit/test_diffusion.py index 3e1ea2620a..79e220acca 100644 --- a/python/test/unit/test_diffusion.py +++ b/python/test/unit/test_diffusion.py @@ -16,6 +16,7 @@ particles. """ + # --------------------------------------------------------------------------------------- # recipe class class recipe(A.recipe): @@ -59,7 +60,6 @@ def event_generators(self, gid): # --------------------------------------------------------------------------------------- # test class class TestDiffusion(unittest.TestCase): - # Constructor (overridden) def __init__(self, args): super(TestDiffusion, self).__init__(args) @@ -72,7 +72,6 @@ def __init__(self, args): def simulate_diffusion( self, cat, _num_segs, _num_cvs_per_seg, _length, _r_1, _r_2=0.0, _r_3=0.0 ): - # --------------------------------------------------------------------------------------- # set the main parameters and calculate geometrical measures num_segs = _num_segs # number of segments @@ -97,7 +96,7 @@ def simulate_diffusion( num_segs * num_cvs_per_seg ) # surface area of one cylindrical CV in µm^2 (excluding the circle-shaped ends, since Arbor does not consider current flux there) volume_tot = ( - np.pi * (radius_1 ** 2 + radius_2 ** 2 + radius_3 ** 2) * length_per_seg + np.pi * (radius_1**2 + radius_2**2 + radius_3**2) * length_per_seg ) # volume of the whole setup in µm^3 volume_per_cv = volume_tot / ( num_segs * num_cvs_per_seg @@ -150,7 +149,6 @@ def simulate_diffusion( } ) elif num_segs == 3: - s = tree.append( A.mnpos, A.mpoint(-1 / 3 * length, 0, 0, radius_1), @@ -302,7 +300,6 @@ def simulate_diffusion( # Test: simulations with equal radii @fixtures.diffusion_catalogue() def test_diffusion_equal_radii(self, diffusion_catalogue): - self.simulate_diffusion( diffusion_catalogue, 1, 600, 10, 4 ) # 1 segment with radius 4 µm From 1bd43cca499b085dfe76a28c2d3e0113c90a9a87 Mon Sep 17 00:00:00 2001 From: Jannik Luboeinski Date: Wed, 16 Aug 2023 15:49:26 +0200 Subject: [PATCH 04/80] Cleanup --- python/test/unit/test_diffusion.py | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/python/test/unit/test_diffusion.py b/python/test/unit/test_diffusion.py index 79e220acca..02ee406764 100644 --- a/python/test/unit/test_diffusion.py +++ b/python/test/unit/test_diffusion.py @@ -10,8 +10,8 @@ Three different morphological structures are considered: 1 segment ("soma only"), 2 segments ("soma with dendrite"), and 3 segments ("soma with two dendrites"). -NOTE: Internally, Arbor only knows concentrations. Thus, particle amounts have to be computed - from concentrations by integrating over the volume of the morphology. The total amount +NOTE: Internally, Arbor only knows concentrations. Thus, particle amounts have to be computed + from concentrations by integrating over the volume of the morphology. The total amount of particles should be conserved unless there is deliberate injection or removal of particles. """ @@ -89,12 +89,6 @@ def simulate_diffusion( radius_3 = 0 length_per_seg = length / num_segs # axial length of a segment in µm - area_tot = ( - 2 * np.pi * (radius_1 + radius_2 + radius_3) * length_per_seg - ) # surface area of the whole setup in µm^2 (excluding the circle-shaped ends, since Arbor does not consider current flux there) - area_per_cv = area_tot / ( - num_segs * num_cvs_per_seg - ) # surface area of one cylindrical CV in µm^2 (excluding the circle-shaped ends, since Arbor does not consider current flux there) volume_tot = ( np.pi * (radius_1**2 + radius_2**2 + radius_3**2) * length_per_seg ) # volume of the whole setup in µm^3 @@ -180,7 +174,7 @@ def simulate_diffusion( } ) else: - raise ValueError(f"Specified number of segments not supported.") + raise ValueError(f"Specified number of segments ({num_segs}) not supported.") morph = A.morphology(tree) # --------------------------------------------------------------------------------------- @@ -248,7 +242,6 @@ def simulate_diffusion( # --------------------------------------------------------------------------------------- # retrieve data and do the testing data_s = sim.samples(hdl_s)[0][0] - times = data_s[:, 0] data_sV = sim.samples(hdl_sV)[0][0] tmp_data = sim.samples(hdl_sV_all)[0][0] data_sV_total = np.zeros_like(tmp_data[:, 0]) From 72b0d6138e18420276590807bd4c9ca348135c80 Mon Sep 17 00:00:00 2001 From: Jannik Luboeinski Date: Wed, 16 Aug 2023 15:51:54 +0200 Subject: [PATCH 05/80] Black reformatting --- python/test/unit/test_diffusion.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/python/test/unit/test_diffusion.py b/python/test/unit/test_diffusion.py index 02ee406764..25825665d6 100644 --- a/python/test/unit/test_diffusion.py +++ b/python/test/unit/test_diffusion.py @@ -174,7 +174,9 @@ def simulate_diffusion( } ) else: - raise ValueError(f"Specified number of segments ({num_segs}) not supported.") + raise ValueError( + f"Specified number of segments ({num_segs}) not supported." + ) morph = A.morphology(tree) # --------------------------------------------------------------------------------------- From 5b9ecd33a49e70c46f4e1acdebcef7743bfc1925 Mon Sep 17 00:00:00 2001 From: Jannik Luboeinski Date: Thu, 24 Aug 2023 15:16:50 +0200 Subject: [PATCH 06/80] Revision: list of dictionaries for stimulation added; small improvements --- .gitignore | 3 - python/test/readme.md | 4 +- python/test/unit/test_diffusion.py | 119 +++++++++++++----- test/unit/diffusion/neuron_with_diffusion.mod | 2 +- .../unit/diffusion/synapse_with_diffusion.mod | 4 +- 5 files changed, 91 insertions(+), 41 deletions(-) diff --git a/.gitignore b/.gitignore index 9178ba0f25..22631b86a3 100644 --- a/.gitignore +++ b/.gitignore @@ -105,6 +105,3 @@ _skbuild # generated by test scripts results - -# environment setting script -set_arbor_env* diff --git a/python/test/readme.md b/python/test/readme.md index 613e455a82..ac4cc673e0 100644 --- a/python/test/readme.md +++ b/python/test/readme.md @@ -39,6 +39,6 @@ In subfolders `unit`/`unit_distributed`: ## Naming convention -- modules: `test_xxxs` (if suitable, ending with `s` since module can consist of multiple classes) -- class(es): `TestXxxs` (if suitable, ending with `s` since class can consist of multiple test functions) +- modules: `test_xxxs` (if applicable use the plural with `s` as modules can comprise multiple classes) +- class(es): `TestXxxs` (if applicable use the plural with `s` as classes can comprise multiple test functions) - functions: `test_yyy` diff --git a/python/test/unit/test_diffusion.py b/python/test/unit/test_diffusion.py index 25825665d6..0b0be11e40 100644 --- a/python/test/unit/test_diffusion.py +++ b/python/test/unit/test_diffusion.py @@ -20,7 +20,14 @@ # --------------------------------------------------------------------------------------- # recipe class class recipe(A.recipe): - def __init__(self, cat, cell, probes, inc_0, inc_1, dec_0): + # Constructor + # - cat: catalogue of custom mechanisms + # - cell: cell description + # - probes: list of probes + # - inject_remove: list of dictionaries of the form [ {"time" :