From ac03756223fbd95824397162b45d82121d12d4d1 Mon Sep 17 00:00:00 2001
From: DrTaDa <tanguy.damart@protonmail.com>
Date: Mon, 13 Jul 2026 08:55:20 +0200
Subject: [PATCH 1/3] Initial commit ICMS

---
 .DS_Store                              | Bin 6148 -> 8196 bytes
 mozaik/.DS_Store                       | Bin 6148 -> 10244 bytes
 mozaik/experiments/microstimulation.py |  58 +++++
 mozaik/sheets/__init__.py              |   7 +-
 mozaik/sheets/direct_stimulator.py     | 327 +++++++++++++++++++++++--
 mozaik/sheets/vision.py                |  28 ++-
 6 files changed, 397 insertions(+), 23 deletions(-)
 create mode 100644 mozaik/experiments/microstimulation.py

diff --git a/.DS_Store b/.DS_Store
index a944017f331dada584b5ca736be4a7c48ecf9f48..91dd1d7a6aeea8b7737935f94764c5b176841c5b 100644
GIT binary patch
literal 8196
zcmeHMYitx%6u#%Sw09_D2TG-Y3}aUXq-Y<4P+r?^t5O~%?Dn-h>g>+gPMFTro!Q+&
zB`HxcfoO~merXg<{6+8;6BB}uL=C}b3^5|cANs=|#6%PEkAFOO?kus3{UtFF=O%N{
zoOAEFGw00r&D}Fg2!YPL(LzWAA%x+<rBZ|4KQunhkE<FDEtMz%{Ry$i{!Gyu&)9C}
zs@aht%0QHXC<9Rjq6|bCxE?Zqdp2)uhJ9a*MtzilC<Fh?4Dk1Zh6k7NP)-RMygIN8
zN&u9+1PB_QV;vCqXrYXUa!Sa+71tEi1B9*+mKY%1aa`lmQN}|#C8TgCAlwOrSB9`c
zf#B-srv=jq#D$FdC<9RjuFe3bPZb#>E;&f_vikkq%(!iNesl9TKr1R|RaLWUR>Srt
zN8Mu5*ZsWLspoccE#<m)rgW}n`Yp3Ktjtfij;~venF|~&!^A0j1}w*Li=A%XF}%Px
zDF&Ftq+z9Da<Z*$ZIjlvru9&hHo2~~xk<aHZOx%W6048hy?)p5etXPukBTrLd=(J6
zEUaVYG1*=rA-+jeWJ)nzEX5?vu9+j>q{>6--Fwn|2d%8@B&~zyuu|(`Xj@R-nKKI*
z(m>X=^J9)Ttjx{ob~bP8zS&{hd@|aUyl3t*(PkCQE*wEqt(VQYw!PE!EZ=gSp+VoA
z@Jnk?8yZY|Zq82x>!<aJ)RlSj_=S7Nxx3eI<~)?Kc(aICI#t`C-M-?kwzV5O?(5mH
z>%1~|UY%SoFBm9To|Uo9zM^ILqeHrvH66opM$(vS$G0A`%$zhAnxvGYkD2G|7A{(R
z%Mw*BpIusIFrTp=&RhNgCB~SzMx%P0q`brS7#Qbv*WfD}!HC+ZE@x~vw<mHY_QEW$
zR98ug&Qc!cN-WVVw@Qj9rS>BpgzIYgUZxmQpRH$1TbSDA4NTc84HQJ=!dNfInKYPn
z$Hp)iPuLP_H&f2Cef+}PbkFY{H=WWgCDbj9rM*!-XZ98mjU|umj0Mja?cAl?uzWWc
zWTUrHRfqUe9fn~V81&+Wi|W*r+Rv1)(MPeIf`$0UXCY|7q*w`ovG|?hH1>g^1!qB$
z=a<Gq$#=V@WI1UkounUGH%T5N$H@tDl6*+ckx$8Y<R@~8{6_wSDyW9pPy-4qgr%?o
zTA&>^z(&{v_rn9w16yGbbjUy!T$q3ZFbR*qVK@rM;7NE2o`o0SMR*BbhS%VAcmv*o
zx8XE=1RujE@Fjc;7vLg%50~I~xD0<#ptETWokLZ+jJD8L+D13fjdU}O(<DvNKDv*N
z&~cg<G0+nG_DzqC9zlKD8wJ@Hei5W+W?-KtZ|T{(t!Mj<#J&r)oNx2%Z(h)F>m4hb
zSFhXDc?FdM#3kJRoeu#8rV%HgK!B2R_zGgmS|4)6oS2E9N;T&)JI`k~#$-*Ig)9>4
zvgN2&A%le4v`W(yio_LaORFX;6-Y}d@LjVCRS1>#^|Gu~B0<I2W))QuDv2&!hSQB*
zd*J^-ekQ+=%K!-cI%t3;&;+YtJ#-+@lhB7?AA&T3-GmVs1q*WELjj5i{loAmJPt?T
z33wWwf#=|PIF7(S0k6QTa1u@-_}@hMpN2E=9=s1{;T(JhpTieOgRhVVKOzZ!Ed%zB
zGGKGSPXl%#nR9H{89{rA4AWh%bFD($UjwiIyJ!CXf1O7W-CmS|C<Ff{15nwO?CQi7
zW_pY9T3f^8K|FZl^QMFhT(Aquag;$hj`GSMhIp>wRPGZG<&=<U!S<hj2nc@JUxQ!Y
N(f%Lp|2Odd{}+W|hyefq

delta 175
zcmZp1XfcprU|?W$DortDU=RQ@Ie-{MGjdEU6q~50$jCG?VE1GL8J5ZX0uqedCu;~4
zs)<)u8yTAGC|DZQ>L^rO8X4#)m>8SY)^c))D(hPZ#b@W_=H+)zJ|G~^*frT){KCeD
z2-d~y92|noK!brmfE!4-g6!B>_?>w&zl@~<NSKj<feA{pfM_5G<Qk9z88*lB%wYxq
DNZuee

diff --git a/mozaik/.DS_Store b/mozaik/.DS_Store
index e7926137e04bac406b1509347078eacd360faaa2..ef4bd066a94baadd2c231a9954250a209a2b4b52 100644
GIT binary patch
literal 10244
zcmeHMYitx%6uxI#+L@NIQ%a?P3}aWN$fA8ff%4jRi`w!iVYjcPP-l0>cEWU~&dlx>
zD@lop2}EOj@Jpj;;xB@)n3xcJBx(pgV~7zk{?H%(ASRlKfBfUQb7u)%EI$ymBr}_P
z&zy7bxwH4&^W8J&4j}|OvRV@%3L%8YLa8(ttADs?J1tjSTp+1H3aC$rPBO$JHgSl1
zMek4$S|GGQXo1iIp#?$<%mxc!%%&G#$!J)I76>g6THule=>DK!p)?lIDMtOR11moT
zKw64o|H6CB1AHANps|2XG3u-0m@IpMKox=%1DHCtE3`Y(SU{&3W$FZ&I)UKH5M(Ie
zKOOZHzdL~#qhTFdAhf`h7NEI%1!nCwvG8|U{+`ae<7v}QWBYsQano?U#>Q`fD=D2*
zR?e66mHeK>sGU!Es+V;;RA)EklD2K83-@Ze&(QP3()_q>d8%RQj;~;7I(8ZAH!RK0
zci34=bA8<eYv2W57?$cLCtF(9MwOP<<^xe>a$R#{R9V&1df<S-*GBGLziW7(IcC|1
z*>oV-5{Otdy(2TtEthaC)3!#io$6NNyeG&w)($bq7*Y#QG2H*#0;|@pi7L&l>;9z$
zW|*Dn=HhI9OVWk)JzrQ~p`voGc#A9!rgjgd_6!&q+e#RZ=)+Q#i_6@F>`q6|;qvxp
zY%@D%xx-R*Ml~~8Q}y(A)1;kICS_fHmyR+cr*~otimKd<W1Hqq+ciAHwgv}0eZnh@
zJ*{aV<=T!H_s37E6Uj^c=FtcDj8k>DPj_6Tkys;Jd*MV?opR^$ds^0RY`?F2$F8$d
z^}HIfR$S1ZGh8EW>b-eG^F{|%H=|pcVU48F?JUoD$j}|38lpndQpfbOH47IlzHNyt
z7xylVGLTIh4`&T;zZBtl_KbS@4ncaCAJTB1+ii`GsQE2&z1+a_yQw_x=vXtcxI$hj
zNGhLn(N=7TMzL9tTp_s+vB6YVi)(mE6M9WGt(#2LCT`%RtwMi}&728q#TYLPWbCmq
zoQ%t4ad`_bo#yw_2X9kduV-Ah3V$gsckz769aSB@Cy&T1Xl&<se~nR2Rkmn`XFI+%
zdg^6)kPg+ZX}X4!Uc7KojhvMGc<CF|kqj$m;9h`h!L|fmh_DqQBoe!u-9|pt6#vdQ
z^87-3K<OL6Qqn-$NC)Xd98QwQ$x(8Q948-<GvqV!J^7hjAitBppbW~P0xCg*g|HNs
zLld;Y2G|Ij;C^@jx?w8}fC_2IfDIF{A12{ZI0%Q~2s{N(!*lQ=yaX@9EATqJ0dK<F
z@D7}WkKq&e6uyG*;2fNXAK(J~0T<y<4!8=glAFuP+%m3-Yvx+G4cta<GZ*6$T$1bM
z_HrZKIG1H};0im@-syAW4k14?4ugG`J2e9i6|t*(>$dLgHxmxeRZ-;4uf26a-R*a+
zXk5K+Q^zI$WbRz>;LDNWgJBwkd@%U_QtZeKqN37^XwfHTLaJ0wQO%#Fz3U^QBF#ZO
zF?Lx4vSo-R#zt2vio_v=8QavXh*AjxRPc3MS0W1m)3#m|rBZ|{YuhX%O9Cd|iNjFf
z@z);$KayX_ujC>CLZAlfU<pKFHLQnrbom7IqT3Hb3f*3Z5f}vn9Pl6qc?7{hcnqF^
zL+~U#1JA<q@B$n~2poe~;Wao8ClCg2ArMZ&DR>_~fYWdWK8G*hO9aH%2#B8$4!;$-
z_>LkMze*r23b$qWabOgKtUeop>~(t0m^P;{w1%NI46XlRXe~xx7*GqZnqfe_>O19Z
z52(`|yC{rKkH6Q#l`tiRp8s!|`TPI13?1r(76>hHlUM*sI}@E9_-700Yw`c4XYC4>
z2eHtL_M2kVSHa3Zjz{Vr$0I|J;~TKhzQNgye;*5EQ;hoR*vD5#%1-ld_{V@S|G&}a
MaAE#G%k%%g0m=&Nv;Y7A

delta 114
zcmZn(XfcprU|?W$DortDU=RQ@Ie-{MGqg`E6q~50$jCG?z<RQQ3d>}E0ntf<YLoLN
z_$MEgciPzShHWuB2ZtatNHI_kH;`}z;`WV&-<c=#s~9peOt#^Xo}8e<&I+=FVRJms
G9A*Fj%NFYZ

diff --git a/mozaik/experiments/microstimulation.py b/mozaik/experiments/microstimulation.py
new file mode 100644
index 000000000..293e527ec
--- /dev/null
+++ b/mozaik/experiments/microstimulation.py
@@ -0,0 +1,58 @@
+from parameters import ParameterSet
+from collections import OrderedDict
+from copy import deepcopy
+
+from mozaik.experiments import Experiment
+from mozaik.stimuli import InternalStimulus
+from mozaik.tools.distribution_parametrization import MozaikExtendedParameterSet
+from mozaik.sheets.direct_stimulator import IntraCorticalMicroStimulation
+
+
+class IntraCorticalMicroStimulationStimulus(Experiment):
+    """
+
+    Parameters
+    ----------
+    model : Model
+          The model on which to execute the experiment.
+
+    Other parameters
+    ----------------
+
+    duration : float
+                Duration of the microstimulation.
+
+    sheet_list : int
+                The list of sheets in which to do stimulation.
+
+    stimulator_parameters: ParameterSet
+                See class IntraCorticalMicroStimulation in direct_stimulator.py
+    """
+
+    required_parameters = ParameterSet({
+        'duration': float,
+        'sheet_list': list,
+        "stimulator_parameters": ParameterSet,
+    })
+
+    def __init__(self, model, parameters):
+        Experiment.__init__(self, model, parameters)
+
+        stimulators = OrderedDict()
+        for sheet in parameters.sheet_list:
+            _params = MozaikExtendedParameterSet(deepcopy(parameters))
+            stimulators[sheet] = [
+                IntraCorticalMicroStimulation(model.sheets[sheet], _params.stimulator_parameters)
+            ]
+        self.direct_stimulation = [stimulators]
+
+        _params = MozaikExtendedParameterSet(deepcopy(parameters))
+        self.stimuli.append(
+            InternalStimulus(
+                frame_duration=_params.duration, 
+                duration=_params.duration,
+                trial=0,
+                direct_stimulation_name='IntraCorticalMicroStimulation',
+                direct_stimulation_parameters=_params.stimulator_parameters
+            )
+        )
diff --git a/mozaik/sheets/__init__.py b/mozaik/sheets/__init__.py
index 1e3fb92a5..7cf67aa49 100644
--- a/mozaik/sheets/__init__.py
+++ b/mozaik/sheets/__init__.py
@@ -294,6 +294,11 @@ def get_data(self, stimulus_duration=None):
 
         """
 
+        # I am not sure I understand the implication of the new version of the present function
+        # For the ICMS paper I need a to record arbitrary signals. For example, here is what I usually record:
+        # block = self.pop.get_data(['spikes', 'v', 'g_exc', 'g_inh', 'gsyn_inh_a', 'gsyn_inh_b', 'gsyn_exc', 'gsyn_inh', 'k_trace', 'w'], clear=True)
+        # How should I modify the following code to achieve that ?
+
         if self.parameters.cell.native_nest:
             vm_name = ['V_m']
         else:
@@ -307,7 +312,7 @@ def get_data(self, stimulus_duration=None):
             gsyn_names = ['g_ex', 'g_in']
         else: 
             gsyn_names = ['gsyn_exc', 'gsyn_inh']
-            
+
         block = None
         steps = self.model.parameters.steps_get_data 
         if steps:
diff --git a/mozaik/sheets/direct_stimulator.py b/mozaik/sheets/direct_stimulator.py
index a5b9e7d44..e2b0e33f0 100644
--- a/mozaik/sheets/direct_stimulator.py
+++ b/mozaik/sheets/direct_stimulator.py
@@ -7,35 +7,39 @@
 each direct stimulator specifies how the given population is stimulated. In general each population can have several
 stimultors.
 """
-from mozaik.core import ParametrizedObject
-from parameters import ParameterSet
+from builtins import zip
+import io
+import math
+import matplotlib
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+from numba import jit
 import numpy
 import numpy as np
 import numpy.random
-import mozaik
-from mozaik.tools.stgen import StGen
-from mozaik import load_component
-from pyNN.parameters import Sequence
-from mozaik import load_component
-import math
-from mozaik.tools.circ_stat import circular_dist ,circ_mean
+from neo.core.analogsignal import AnalogSignal as NeoAnalogSignal
 import pylab
-from scipy.integrate import odeint
+from parameters import ParameterSet
 import pickle
+from probeinterface import read_prb
+from pyNN.parameters import Sequence
+from pyNN.connectors import FromListConnector
+import quantities as qt
 import scipy.interpolate
-from mpl_toolkits.mplot3d import Axes3D
+from scipy.integrate import odeint
+from scipy.spatial.distance import cdist
+from scipy import signal
+from mpi4py import MPI
+
+import mozaik
+from mozaik.tools.stgen import StGen
+from mozaik.core import ParametrizedObject
+from mozaik import load_component
+from mozaik.tools.circ_stat import circular_dist
 from mozaik.controller import Global
-import matplotlib
-from mozaik.analysis.analysis import SingleValue, AnalogSignalList
-from neo.core.analogsignal import AnalogSignal as NeoAnalogSignal
-import quantities as qt
 from mozaik.tools.units import *
-import io
-from numba import jit
-
-from builtins import zip
-
-from mpi4py import MPI
+from mozaik.analysis.analysis import SingleValue, AnalogSignalList, PerNeuronValue
+from mozaik.analysis.data_structures import SingleValueList
 
 mpi_comm = MPI.COMM_WORLD
 
@@ -1041,3 +1045,284 @@ def single_pixel(sheet, coor_x, coor_y, update_interval, parameters):
     )
     signals[np.where(x == coor_x)[0][0], np.where(y == coor_y)[1][0], :] = 1
     return signals
+
+class IntraCorticalMicroStimulation(DirectStimulator):
+    """ This class instantiates a multi-electrode array using the
+    probeinterface library. For each electrode, a subset of neurons is
+    selected at random using a distance dependent distribution centered on
+    the electrode.
+    Note that a neuron can be part of the selected population of multiple electrodes.
+
+    For now this is not mpi optimize NOR TO BE USED WITH MPI.
+    """
+
+    required_parameters = ParameterSet({
+        "current_amplitude": float,       # uA
+        "pulse_per_train": int,
+        "intertrain_interval": float,     # ms
+        "pulse_frequency": float,         # Hz
+        "k_increment_ICMS_factor": float,
+        "connectivity_profiles_path": str,
+        "probe_path": str,
+        "probe_position_offset": list,
+        "probe_rotation_angle": float,
+        "probe_rotation_axis": list,
+        "probe_2d_planes": str,
+        "probe_active_electrodes": list,
+        "recruitment_offset": float,
+        "recruitment_gain_exc": float,
+        "recruitment_gain_inh": float,
+        "stimulator_seed": int
+    })
+
+    def __init__(self, sheet, parameters):
+        """Init"""
+        DirectStimulator.__init__(self, sheet, parameters)
+
+        # ICMS specific random generator. This is used to make sure that the same cells
+        # are selected to be stimulated between runs that use the same electrode(s) and amplitude.
+        self.rng = numpy.random.default_rng(self.parameters.stimulator_seed)
+
+        self.probe = None
+        self.stimulated_cells = {}
+        self.random_propensities = None
+
+        self.instantiate_probe()
+        self.select_stimulated_cells()
+
+    def instantiate_probe(self):
+        """Instantiate a probe based on its description. Please refer to
+        https://probeinterface.readthedocs.io/en/main/format_spec.html for the
+        description file format."""
+
+        self.probe = read_prb(self.parameters.probe_path).probes[0]
+        self.probe = self.probe.to_3d(axes=self.parameters.probe_2d_planes)
+
+        if self.parameters.probe_rotation_angle:
+            self.probe.rotate(
+                theta=self.parameters.probe_rotation_angle,
+                center=None,
+                axis=self.parameters.probe_rotation_axis
+            )
+
+        if any(self.parameters.probe_position_offset):
+            self.probe.move(translation_vector=self.parameters.probe_position_offset)
+
+    def generate_activation_profile(self):
+        """Generate an activation connection profile, that is, a probability for a cell to be
+        activited as a function of distance. It is computed as the convolution of an axonal activation
+        profile (the probability that an axon gets activated at a given distance from the electrode) and 
+        a pre-computed connectivity profile (the probability that this axon is connected to a cell at a
+        given distance from the electrode).
+
+        Note that this profile is independent of the current amplitude."""
+
+        # Get the pre-computer connectivity profile
+        with open(self.parameters.connectivity_profiles_path, "rb") as fp:
+            connectivity_profile = numpy.array(pickle.load(fp)[self.sheet.name])
+
+        def axonal_activation_profile(distance, a=2.67135221, b=2.17589264, c=0.01009279):
+            """This function gives the density of activated axon as a function of the distance
+            to the stimulated electrode. The parameter a, b and c have been fitted on
+            Kumaravelu et al. (2022) Fig5D"""
+            return a / (np.abs(distance)**b + c)
+
+        # Define the unified grid we will work on
+        max_dist = 4.0 # in mm
+        num_points = 4001
+        distance_axis = np.linspace(-max_dist, max_dist, num_points)
+
+        # Create a symmetric version of the experimental data
+        connect_distances = np.array(connectivity_profile[0]) / 1000. # From mm to um
+        connect_probabilities = np.array(connectivity_profile[2]) #[2]
+        sym_connect_distances = np.concatenate([-connect_distances[:0:-1], connect_distances])
+        sym_connect_probabilities = np.concatenate([connect_probabilities[:0:-1], connect_probabilities])
+
+        # Sort the data for np.interp
+        sort_indices = np.argsort(sym_connect_distances)
+        sorted_distances = sym_connect_distances[sort_indices]
+        sorted_probabilities = sym_connect_probabilities[sort_indices]
+
+        # Interpolate the connection profile onto our uniform grid
+        connection_profile_interpolated = np.interp(
+            distance_axis,
+            sorted_distances,
+            sorted_probabilities
+        )
+
+        # Generate the axonal activation profile for the same grid
+        activation_profile = axonal_activation_profile(distance_axis)
+    
+        # Normalize the kernel such as to redistribute the probabilities without
+        # changing the overall sum of probabilites
+        activation_profile = activation_profile / np.sum(activation_profile)
+
+        # Perform the convolution and normalize it to one
+        convolved_result = signal.convolve(
+            activation_profile,
+            connection_profile_interpolated,
+            mode='same'
+        )
+
+        return [distance_axis[2000:] * 1000.0, convolved_result[2000:]] # Back to um
+        # This is to check the result without the convolution:
+        #return [distance_axis[2000:] * 1000.0, connection_profile_interpolated[2000:] / np.max(connection_profile_interpolated)]
+
+    def select_stimulated_cells(self):
+        """For each electrode of the multi-electrode array, select a subset of neurons
+        that get activate by the electrode. The selection is based on the distance
+        between the electrode and the neurons.
+
+        The selection process also ensures that for a given seed, the set of cells
+        activated at a higher amplitude is a superset of the cells activated at lower
+        amplitudes.
+        """
+
+        # Generate the activation profile
+        activation_profile = self.generate_activation_profile()
+
+        # Get the positions of the neurons
+        cells_positions = self.sheet.pop.positions.T
+        cells_positions = [list(self.sheet.vf_2_cs(cp[0], cp[1])) + [cp[2]] for cp in cells_positions]
+        n_cells = len(cells_positions)
+
+        # Get the positions of the electrodes
+        electrodes_positions = self.probe.contact_positions
+        n_electrodes = len(electrodes_positions)
+
+        # Compute the distance between each neuron and each electrode
+        cells_contact_distances = cdist(numpy.array(cells_positions), numpy.array(electrodes_positions))
+
+        # Assign a fixed random activation propensity score to each cell-electrode pair
+        # This propensity is a phenomenological representation of how close the 
+        # axon of a neuron passes by each electrode.
+        self.random_propensities = self.rng.random(size=(n_cells, n_electrodes))
+
+        # Calculate distance-dependent activation probabilities
+        _idx_distribution = numpy.digitize(cells_contact_distances, activation_profile[0]) - 1
+        _idx_distribution = numpy.clip(_idx_distribution, 0, len(activation_profile[1]) - 1)
+
+        # Calculate the amplitude-specific scaling factor, which increases as sqrt(current_amplitude)
+        if "Inh" in self.sheet.name:
+            recruitment_gain = self.parameters.recruitment_gain_inh
+        else:
+            recruitment_gain = self.parameters.recruitment_gain_exc
+        recruitment_rescale = numpy.sqrt(numpy.clip(
+            recruitment_gain * (self.parameters.current_amplitude - self.parameters.recruitment_offset),
+            0,
+            numpy.inf
+        ))
+
+        # This calculates the probability P(activate | distance, current_amplitude)
+        scaled_probabilities = activation_profile[1][_idx_distribution] * recruitment_rescale
+
+        # Compare the fixed random score with the amplitude-dependent probability
+        # Activate if random_propensity < P(activate | distance, current_amplitude)
+        _mask = self.random_propensities < scaled_probabilities
+
+        # Only keep the active electrodes
+        for j in range(n_electrodes):
+            if j not in self.parameters.probe_active_electrodes:
+                _mask[:, j] = False
+
+        # Build the dictionary of stimulated cells
+        self.stimulated_cells = {}
+        for i, cell_id in enumerate(self.sheet.pop.all_cells):
+            # If any *active* electrode activates a cell, add it to the list
+            _idx_electrodes = numpy.nonzero(_mask[i, :])[0]
+            if len(_idx_electrodes):
+                self.stimulated_cells[cell_id] = _idx_electrodes
+
+        logger.info(
+             f"ICMS: in population {self.sheet.name}, activating {len(self.stimulated_cells)} cells " +
+             f"({100 * len(self.stimulated_cells) / n_cells:.2f}% of the population) " +
+             f"at amplitude {self.parameters.current_amplitude}."
+        )
+
+    def prepare_stimulation(self, duration, offset):
+        """
+        Prepares the stimulation for the next period of model simulation lasting `duration` seconds.
+
+        Parameters
+        ----------
+        duration : double (seconds)
+                 The period for which to prepare the stimulation
+
+        offset : double (seconds)
+               The current simulator time.
+        """
+
+        k_increment = self.parameters.k_increment_ICMS_factor * self.parameters.current_amplitude
+        if "Exc" in self.sheet.name:
+            k_increment = 0
+
+        ICMS_ISI = 1000. / self.parameters.pulse_frequency
+        ISIs = []
+        propensities = []
+        for i, cell_id in enumerate(self.sheet.pop.all_cells):
+            if cell_id in self.stimulated_cells:
+                ISIs.append(ICMS_ISI)
+                activating_electrodes = self.stimulated_cells[cell_id]
+                min_random_draw = np.min(self.random_propensities[i, activating_electrodes])
+                # Cells that have lower propensities are more sensitive to K accumulation
+                propensities.append(1.0 - min_random_draw)
+            else:
+                # A stimulation ISI of 0 means no ICMS-induced spikes
+                ISIs.append(0)
+                propensities.append(0)
+
+        self.sheet.pop.set(propensity_ICMS=propensities)
+        self.sheet.pop.set(pulse_ISI_ICMS=ISIs)
+        self.sheet.pop.set(intertrain_interval_ICMS=self.parameters.intertrain_interval)
+        self.sheet.pop.set(pulse_per_train_ICMS=self.parameters.pulse_per_train)
+        self.sheet.pop.set(k_increment_ICMS=k_increment)
+
+        logger.info(
+            f"ICMS: starting stimulation at frequency {self.parameters.pulse_frequency}Hz, with an inter-train interval of {self.parameters.intertrain_interval} and {self.parameters.pulse_per_train} pulses per train."
+        )
+
+    def inactivate(self, offset):
+        """
+        Ensures any influences of the stimulation are inactivated for subsequent simulation of the model.
+
+        Parameters
+        ----------
+        offset : double (seconds)
+               The current simulator time.
+
+        Note that a subsequent call to prepare_stimulation should 'activate' the stimulator again.
+        """
+        self.sheet.pop.set(pulse_ISI_ICMS=0)
+        self.sheet.pop.set(intertrain_interval_ICMS=0)
+        self.sheet.pop.set(pulse_per_train_ICMS=0)
+        self.sheet.pop.set(k_increment_ICMS=0)
+        logger.info(f"ICMS: ending stimulation.")
+
+    def save_to_datastore(self, data_store, stimulus):
+        """Stores the electrode positions and the list of cells activated by each electrode"""
+        
+        active_electrode = '_'.join(str(e) for e in self.parameters.probe_active_electrodes)
+
+        # This "metadata" is a trick used in case there are multiple ICMS experiments in the same run.
+        # It does not work if there are several experiments with the exact same parameters.
+        metadata = f"__{self.parameters.current_amplitude}__{self.parameters.pulse_frequency}__{active_electrode}"
+
+        electrode_per_cell = [list(el)[0] for el in self.stimulated_cells.values()]
+        data_store.full_datastore.add_analysis_result(
+            PerNeuronValue(
+                values=electrode_per_cell,
+                idds=[int(idd) for idd in self.stimulated_cells.keys()],
+                value_name='electrode_active_per_cell' + metadata,
+                sheet_name=self.sheet.name,
+                value_units="None"
+            )
+        )
+
+        data_store.full_datastore.add_analysis_result(
+            SingleValueList(
+                values=self.probe.contact_positions,
+                values_unit='um',
+                value_name='probe_electrode_positions' + metadata,
+                sheet_name=self.sheet.name
+            )
+        )
diff --git a/mozaik/sheets/vision.py b/mozaik/sheets/vision.py
index a68d4c2b3..b1e2c016f 100644
--- a/mozaik/sheets/vision.py
+++ b/mozaik/sheets/vision.py
@@ -303,8 +303,34 @@ def __init__(self, model, parameters):
 
         else:
             if self.parameters.cell.native_nest:
+
+                # Hack to get the NESTML model to work:
+                from pyNN.standardmodels.base import build_translations
+                celltype = self.sim.native_cell_type(self.parameters.cell.model)
+                celltype.translations = build_translations(
+                    ('tau_syn_E',  'tau_syn_exc'),
+                    ('tau_syn_I',  'tau_syn_inh_a'),
+                )
+                celltype.units['w'] = "pA"
+                celltype.units['v'] = "mV"
+                celltype.units['k_trace'] = "mM"
+                celltype.units['gsyn_exc'] = "nS"
+                celltype.units['gsyn_inh'] = "nS"
+                celltype.units['gsyn_inh_a'] = "nS"
+                celltype.units['gsyn_inh_b'] = "nS"
+
+                # Ideally, it should go something like:
+                #from pyNN.standardmodels.base import build_translations
+                #celltype = self.sim.native_cell_type(self.parameters.cell.model)
+
+                #if len(self.parameters.cell.nestml_translation):
+                #    celltype.translations = build_translations(*list(self.parameters.cell.nestml_translation.items()))
+                #if len(self.parameters.cell.nestml_units):
+                #    for k, v in parameters.nestml_units:
+                #        celltype.units[k] = v
+
                 self.pop = self.sim.Population(int(parameters.sx * parameters.sy/1000000 * parameters.density),
-                                                   self.sim.native_cell_type(self.parameters.cell.model),
+                                                   celltype,
                                                    self.parameters.cell.params,
                                                    structure=rs,
                                                    initial_values=self.parameters.cell.initial_values,

From 02b2f7e0cfa3de07d3f3dd7108072191b64d0d07 Mon Sep 17 00:00:00 2001
From: DrTaDa <tanguy.damart@protonmail.com>
Date: Mon, 13 Jul 2026 09:02:48 +0200
Subject: [PATCH 2/3] Remove and ignore .DS_Store files

---
 .DS_Store        | Bin 8196 -> 0 bytes
 .gitignore       |   1 +
 mozaik/.DS_Store | Bin 10244 -> 0 bytes
 3 files changed, 1 insertion(+)
 delete mode 100644 .DS_Store
 delete mode 100644 mozaik/.DS_Store

diff --git a/.DS_Store b/.DS_Store
deleted file mode 100644
index 91dd1d7a6aeea8b7737935f94764c5b176841c5b..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001

literal 8196
zcmeHMYitx%6u#%Sw09_D2TG-Y3}aUXq-Y<4P+r?^t5O~%?Dn-h>g>+gPMFTro!Q+&
zB`HxcfoO~merXg<{6+8;6BB}uL=C}b3^5|cANs=|#6%PEkAFOO?kus3{UtFF=O%N{
zoOAEFGw00r&D}Fg2!YPL(LzWAA%x+<rBZ|4KQunhkE<FDEtMz%{Ry$i{!Gyu&)9C}
zs@aht%0QHXC<9Rjq6|bCxE?Zqdp2)uhJ9a*MtzilC<Fh?4Dk1Zh6k7NP)-RMygIN8
zN&u9+1PB_QV;vCqXrYXUa!Sa+71tEi1B9*+mKY%1aa`lmQN}|#C8TgCAlwOrSB9`c
zf#B-srv=jq#D$FdC<9RjuFe3bPZb#>E;&f_vikkq%(!iNesl9TKr1R|RaLWUR>Srt
zN8Mu5*ZsWLspoccE#<m)rgW}n`Yp3Ktjtfij;~venF|~&!^A0j1}w*Li=A%XF}%Px
zDF&Ftq+z9Da<Z*$ZIjlvru9&hHo2~~xk<aHZOx%W6048hy?)p5etXPukBTrLd=(J6
zEUaVYG1*=rA-+jeWJ)nzEX5?vu9+j>q{>6--Fwn|2d%8@B&~zyuu|(`Xj@R-nKKI*
z(m>X=^J9)Ttjx{ob~bP8zS&{hd@|aUyl3t*(PkCQE*wEqt(VQYw!PE!EZ=gSp+VoA
z@Jnk?8yZY|Zq82x>!<aJ)RlSj_=S7Nxx3eI<~)?Kc(aICI#t`C-M-?kwzV5O?(5mH
z>%1~|UY%SoFBm9To|Uo9zM^ILqeHrvH66opM$(vS$G0A`%$zhAnxvGYkD2G|7A{(R
z%Mw*BpIusIFrTp=&RhNgCB~SzMx%P0q`brS7#Qbv*WfD}!HC+ZE@x~vw<mHY_QEW$
zR98ug&Qc!cN-WVVw@Qj9rS>BpgzIYgUZxmQpRH$1TbSDA4NTc84HQJ=!dNfInKYPn
z$Hp)iPuLP_H&f2Cef+}PbkFY{H=WWgCDbj9rM*!-XZ98mjU|umj0Mja?cAl?uzWWc
zWTUrHRfqUe9fn~V81&+Wi|W*r+Rv1)(MPeIf`$0UXCY|7q*w`ovG|?hH1>g^1!qB$
z=a<Gq$#=V@WI1UkounUGH%T5N$H@tDl6*+ckx$8Y<R@~8{6_wSDyW9pPy-4qgr%?o
zTA&>^z(&{v_rn9w16yGbbjUy!T$q3ZFbR*qVK@rM;7NE2o`o0SMR*BbhS%VAcmv*o
zx8XE=1RujE@Fjc;7vLg%50~I~xD0<#ptETWokLZ+jJD8L+D13fjdU}O(<DvNKDv*N
z&~cg<G0+nG_DzqC9zlKD8wJ@Hei5W+W?-KtZ|T{(t!Mj<#J&r)oNx2%Z(h)F>m4hb
zSFhXDc?FdM#3kJRoeu#8rV%HgK!B2R_zGgmS|4)6oS2E9N;T&)JI`k~#$-*Ig)9>4
zvgN2&A%le4v`W(yio_LaORFX;6-Y}d@LjVCRS1>#^|Gu~B0<I2W))QuDv2&!hSQB*
zd*J^-ekQ+=%K!-cI%t3;&;+YtJ#-+@lhB7?AA&T3-GmVs1q*WELjj5i{loAmJPt?T
z33wWwf#=|PIF7(S0k6QTa1u@-_}@hMpN2E=9=s1{;T(JhpTieOgRhVVKOzZ!Ed%zB
zGGKGSPXl%#nR9H{89{rA4AWh%bFD($UjwiIyJ!CXf1O7W-CmS|C<Ff{15nwO?CQi7
zW_pY9T3f^8K|FZl^QMFhT(Aquag;$hj`GSMhIp>wRPGZG<&=<U!S<hj2nc@JUxQ!Y
N(f%Lp|2Odd{}+W|hyefq

diff --git a/.gitignore b/.gitignore
index 0651a45b2..664cce454 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,4 @@ build
 cover
 .coverage
 AGENTS.md
+.DS_Store
diff --git a/mozaik/.DS_Store b/mozaik/.DS_Store
deleted file mode 100644
index ef4bd066a94baadd2c231a9954250a209a2b4b52..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001

literal 10244
zcmeHMYitx%6uxI#+L@NIQ%a?P3}aWN$fA8ff%4jRi`w!iVYjcPP-l0>cEWU~&dlx>
zD@lop2}EOj@Jpj;;xB@)n3xcJBx(pgV~7zk{?H%(ASRlKfBfUQb7u)%EI$ymBr}_P
z&zy7bxwH4&^W8J&4j}|OvRV@%3L%8YLa8(ttADs?J1tjSTp+1H3aC$rPBO$JHgSl1
zMek4$S|GGQXo1iIp#?$<%mxc!%%&G#$!J)I76>g6THule=>DK!p)?lIDMtOR11moT
zKw64o|H6CB1AHANps|2XG3u-0m@IpMKox=%1DHCtE3`Y(SU{&3W$FZ&I)UKH5M(Ie
zKOOZHzdL~#qhTFdAhf`h7NEI%1!nCwvG8|U{+`ae<7v}QWBYsQano?U#>Q`fD=D2*
zR?e66mHeK>sGU!Es+V;;RA)EklD2K83-@Ze&(QP3()_q>d8%RQj;~;7I(8ZAH!RK0
zci34=bA8<eYv2W57?$cLCtF(9MwOP<<^xe>a$R#{R9V&1df<S-*GBGLziW7(IcC|1
z*>oV-5{Otdy(2TtEthaC)3!#io$6NNyeG&w)($bq7*Y#QG2H*#0;|@pi7L&l>;9z$
zW|*Dn=HhI9OVWk)JzrQ~p`voGc#A9!rgjgd_6!&q+e#RZ=)+Q#i_6@F>`q6|;qvxp
zY%@D%xx-R*Ml~~8Q}y(A)1;kICS_fHmyR+cr*~otimKd<W1Hqq+ciAHwgv}0eZnh@
zJ*{aV<=T!H_s37E6Uj^c=FtcDj8k>DPj_6Tkys;Jd*MV?opR^$ds^0RY`?F2$F8$d
z^}HIfR$S1ZGh8EW>b-eG^F{|%H=|pcVU48F?JUoD$j}|38lpndQpfbOH47IlzHNyt
z7xylVGLTIh4`&T;zZBtl_KbS@4ncaCAJTB1+ii`GsQE2&z1+a_yQw_x=vXtcxI$hj
zNGhLn(N=7TMzL9tTp_s+vB6YVi)(mE6M9WGt(#2LCT`%RtwMi}&728q#TYLPWbCmq
zoQ%t4ad`_bo#yw_2X9kduV-Ah3V$gsckz769aSB@Cy&T1Xl&<se~nR2Rkmn`XFI+%
zdg^6)kPg+ZX}X4!Uc7KojhvMGc<CF|kqj$m;9h`h!L|fmh_DqQBoe!u-9|pt6#vdQ
z^87-3K<OL6Qqn-$NC)Xd98QwQ$x(8Q948-<GvqV!J^7hjAitBppbW~P0xCg*g|HNs
zLld;Y2G|Ij;C^@jx?w8}fC_2IfDIF{A12{ZI0%Q~2s{N(!*lQ=yaX@9EATqJ0dK<F
z@D7}WkKq&e6uyG*;2fNXAK(J~0T<y<4!8=glAFuP+%m3-Yvx+G4cta<GZ*6$T$1bM
z_HrZKIG1H};0im@-syAW4k14?4ugG`J2e9i6|t*(>$dLgHxmxeRZ-;4uf26a-R*a+
zXk5K+Q^zI$WbRz>;LDNWgJBwkd@%U_QtZeKqN37^XwfHTLaJ0wQO%#Fz3U^QBF#ZO
zF?Lx4vSo-R#zt2vio_v=8QavXh*AjxRPc3MS0W1m)3#m|rBZ|{YuhX%O9Cd|iNjFf
z@z);$KayX_ujC>CLZAlfU<pKFHLQnrbom7IqT3Hb3f*3Z5f}vn9Pl6qc?7{hcnqF^
zL+~U#1JA<q@B$n~2poe~;Wao8ClCg2ArMZ&DR>_~fYWdWK8G*hO9aH%2#B8$4!;$-
z_>LkMze*r23b$qWabOgKtUeop>~(t0m^P;{w1%NI46XlRXe~xx7*GqZnqfe_>O19Z
z52(`|yC{rKkH6Q#l`tiRp8s!|`TPI13?1r(76>hHlUM*sI}@E9_-700Yw`c4XYC4>
z2eHtL_M2kVSHa3Zjz{Vr$0I|J;~TKhzQNgye;*5EQ;hoR*vD5#%1-ld_{V@S|G&}a
MaAE#G%k%%g0m=&Nv;Y7A


From a69216304e7886aba25a6bd5eaa949ea68adf907 Mon Sep 17 00:00:00 2001
From: DrTaDa <tanguy.damart@protonmail.com>
Date: Tue, 14 Jul 2026 12:24:49 +0200
Subject: [PATCH 3/3] Added probeinterface to the requirements

---
 requirements.txt | 1 +
 1 file changed, 1 insertion(+)

diff --git a/requirements.txt b/requirements.txt
index 248fc5aa5..f8d169855 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -27,6 +27,7 @@ param==1.5.1
 parameters==0.2.1
 Pillow==9.1.0
 pluggy==1.0.0
+probeinterface
 psutil==5.9.0
 py==1.11.0
 PyNN==0.10.0