Tutorial 4 - Stimulating nerves with NRV

In this tutorial, we will create a 2-fascicle nerve, populate it with axons and stimulate it with intra- and extra-fascicular electrodes.

As before, we start by importing the NRV package as well as numpy and matplotlib:

import nrv
import matplotlib.pyplot as plt
import numpy as np

Nerve creation

First, we need to create our nerve object using the NRV’s nerve-class. This object contains the geometrical properties of the nerve. NRV currently only supports cylindrical shapes for nerve, thus a diameter (nerve_d) and a length (nerve_l) must be specified at the nerve creation. The Outer_D parameter can also be specified. It refers to the saline solution bath diameter in which the nerve is plunged into.

outer_d = 5         # in mm
nerve_d = 500       # in um
nerve_l = 5000      # in um
nerve = nrv.nerve(length=nerve_l, diameter=nerve_d, Outer_D=outer_d)

Then, we will add two fascicles to the nerve. Fascicles in NRV are cylindrical shapes defined by their diameter and their (y,z) coordinates in space. The (0,0) coordinate is aligned with the center of the nerve. Fascicle are defined with the NRV’s fascicle-class. The ID parameters of the fascicle-object tags each fascicle of the model which will facilitate the post-simulation analysis. Fascicle are incorporated one by one to the nerve-object using the add_fascicle-method. We can now plot a 2-D section of the nerve with the plotmethod of the nerve-object to visualize it.

fasc1_d = 200       # in um
fasc1_y = -100      # in um
fasc1_z = 0         # in um

fasc2_d = 100       # in um
fasc2_y = 100       # in um
fasc2_z = 0         # in um

#create the fascicle objects
fascicle_1 = nrv.fascicle(diameter=fasc1_d,ID=1)
fascicle_2 = nrv.fascicle(diameter=fasc2_d, ID=2)

#Add the fascicles to the nerve
nerve.add_fascicle(fascicle=fascicle_1, y=fasc1_y, z=fasc1_z)
nerve.add_fascicle(fascicle=fascicle_2, y=fasc2_y, z=fasc2_z)

#plot
fig, ax = plt.subplots(1, 1, figsize=(6,6))
nerve.plot(ax)
ax.set_xlabel("z-axis (µm)")
ax.set_ylabel("y-axis (µm)")

Populate fascicles with axons

Now that our nerve is created, it is time to populate them with axons. The first step is to create the population. In this example, we will populate each fascicle with myelinated and unmyelinated axons, with a total of 100 axons in each fascicle. To create a realistic axon diameter distribution, we use the NRV’s create_axon_population-method. The function take for arguments:

• the number of axon in the population n_ax

• the proportion of unmyelinated fibers in the population percent_unmyel

• the myelinated axon distribution M_stat

• the unmyelinated axon distribution U_stat

Available myelinated and unmyelinated axon distributions are described in xxx.

The create_axon_populationmethod returns four numpy arrays:

• axons_diameterswhich contains the diameter of each axon of the population.

• axon_type containing a ‘1’ value for indexes where the axon is myelinated, else ‘0’.

• M_diam_list the diameter of myelinated axons only

• U_diam_list the diameter of unmyelinated axons only

n_ax = 100      #size of the axon population
axons_diameters, axons_type, M_diam_list, U_diam_list = nrv.create_axon_population(n_ax, percent_unmyel=0.7, M_stat="Ochoa_M", U_stat="Ochoa_U",)

Once the population is generated, we can fill the fascicle with it using the fill_with_population-method of the fascicle-object. If the (y,z) coordinates of each axon is not explicitly specified in the fill_with_population-method, the NRV’s build-in axon packing algorithm will be automatically called to place each axon within the fascicle. The delta parameter of method indicates the minimum distance between two axons, and between one axon and the border of the fascicle (in \(\mu m\)).

NOTE

The axon packing algorithm can take several minutes to run for large axon population.

fascicle_1.fill_with_population(axons_diameters, axons_type, delta=5)

Let’s repeat this operation for the 2nd fascicle and plot the nerve again:

axons_diameters, axons_type, M_diam_list, U_diam_list = nrv.create_axon_population(n_ax, percent_unmyel=0.7, M_stat="Ochoa_M", U_stat="Ochoa_U",)
fascicle_2.fill_with_population(axons_diameters, axons_type, delta=5)

#Plot the nerve again.
fig, ax = plt.subplots(1, 1, figsize=(6,6))
nerve.plot(ax)
ax.set_xlabel("z-axis (µm)")
ax.set_ylabel("y-axis (µm)")

We see that now our nerve is populated with fibers. Since the fascicle_1and fascicle_2-objects are attached to the nerve-object, any modification to one of them will be propagated to the nerve. The packed population in fascicle_2 was too large to fit in the fascicle. Any axon outside the fascicle are automatically discarded. In fascicle fascicle_1 the population is however too small to fill the entire fascicle. We can fix this by simple calling the fit_population_to_sizemethod of fascicle_1, where deltaspecifies the minimum distance between the fascicle border and the axons.

NOTE

The axon population can be automatically fitted to the fascicle size by setting the fit_to_sizeboolean parameter to Truein the fill_with_population method.

fascicle_1.fit_population_to_size(delta = 2)
fig, ax = plt.subplots(1, 1, figsize=(6,6))
nerve.plot(ax)
ax.set_xlabel("z-axis (µm)")
ax.set_ylabel("y-axis (µm)")

While we are here, we can also define stimulation parameters of the axons. For example, we can specify the computational model of the myelinated and unmyelinated fibers. You can refer to the previous tutorials for a thorough overview of the fiber’s simulation parameters available.

m_model = 'MRG'
um_model = 'Rattay_Aberham'
u_param = {"model": um_model}
m_param = {"model": m_model}

#For fascicle1
fascicle_1.set_axons_parameters(unmyelinated_only=True,**u_param)
fascicle_1.set_axons_parameters(myelinated_only=True,**m_param)

#For fascicle2
fascicle_2.set_axons_parameters(unmyelinated_only=True,**u_param)
fascicle_2.set_axons_parameters(myelinated_only=True,**m_param)

Extracellular stimulation context

Now we will define everything related to the extracellular stimulation. First, we need to create a FEM_stimulation-object. In this object, we can specify the conductivity of each material of the FEM stimulation. Available material conductivities are specified in xxx.

extra_stim = nrv.FEM_stimulation(endo_mat="endoneurium_ranck",      #endoneurium conductivity
                                 peri_mat="perineurium",            #perineurium conductivity
                                 epi_mat="epineurium",              #epineurium conductivity
                                 ext_mat="saline")                  #saline solution conductivity

Adding intracellular electrodes

First, we will run some simulation with 3 intrafascicular LIFE-like electrodes, using the LIFE_electrode NRV’s object. In NRV, LIFEs are defined by a diameter (life_d), an active-site length (life_length) and a (x,y,z) spatial coordinates. A label and an ID can also be specified to facilitate post-simulation analysis. In this example we aligned the LIFEs x-position to the middle of the nerve, and set their (y,z) coordinates such that:

• LIFE_0 is located inside the nerve but outside the fascicles

• LIFE_1 is located inside fascicle_1

• LIFE_2 is located inside fascicle_2

The electrodes are attached to the extra_stim FEM_stimulation-object with the add_electrode-method. The method also requires to link the electrode to a NRV stimulus-object. For that, we created a dummy stimulus dummy_stimthat we will change later.

life_d = 25                                 #LIFE diamter in um
life_length = 1000                          #LIFE active-site length in um
life_x_offset = (nerve_l-life_length)/2     #x position of the LIFE (centered)

life_y_c_0 = 0                              #LIFE_0 y-coordinate (in um)
life_z_c_0 = 150                            #LIFE_0 z-coordinate (in um)
life_y_c_1 = fasc1_y                        #LIFE_1 y-coordinate (in um)
life_z_c_1 = fasc1_z                        #LIFE_1 z-coordinate (in um)
life_y_c_2 = fasc2_y                        #LIFE_2 y-coordinate (in um)
life_z_c_2 = fasc2_z                        #LIFE_1 z-coordinate (in um)

elec_0 = nrv.LIFE_electrode("LIFE_0", life_d, life_length, life_x_offset, life_y_c_0, life_z_c_0, ID = 0) # LIFE in neither of the two fascicles
elec_1 = nrv.LIFE_electrode("LIFE_1", life_d, life_length, life_x_offset, life_y_c_1, life_z_c_1, ID = 1) # LIFE in the fascicle 1
elec_2 = nrv.LIFE_electrode("LIFE_2", life_d, life_length, life_x_offset, life_y_c_2, life_z_c_2, ID = 2) # LIFE in the fascicle 2

#Dummy stimulus
dummy_stim = nrv.stimulus()
dummy_stim.pulse(0, 0.1, 1)

#Attach electrodes to the extra_stim object
extra_stim.add_electrode(elec_0, dummy_stim)
extra_stim.add_electrode(elec_1, dummy_stim)
extra_stim.add_electrode(elec_2, dummy_stim)

Last, we attach extra_stim-object to the nerve with the attach_extracellular_stimulation-method:

nerve.attach_extracellular_stimulation(extra_stim)

Let’s see how our nerve with electrodes now looks like:

fig, ax = plt.subplots(1, 1, figsize=(6,6))
nerve.plot(ax)
ax.set_xlabel("z-axis (µm)")
ax.set_ylabel("y-axis (µm)")

The three LIFEs now are showing up, and we can make sure that their positions within the nerve are corrects. We also note that axon overlapping with the electrodes are removed.

Simulating the nerve

Now it’s time to run some simulations!

First, we set up a few flags:

• nerve.save_results = False disables the automatic saving of the simulation results in a folder

• nerve.return_parameters_only = False makes sure that all simulation results are returned to the nerve_resultsdictionnary.

• nerve.verbose = True so it looks cool

NOTE

Saving simulation results in a folder and returning simulation parameters only can avoid excessive RAM memory usage for large nerve simulation. By default, nerve.save_results and nerve.return_parameters_only are set to False i.e. results are not saved in a folder and all simulation results are available in nerve_results.

Simulation duration is set with the t_sim parameter (in ms). We can also specify a postproc_function which will be applied to each axon’s simulation results. This is particularly useful to remove unused data and save up some memory. In this example we will use the is_recruited function. More details here xxx.

NOTE

This cell takes several minutes to run.

nerve.save_results = False
nerve.return_parameters_only = False
nerve.verbose = True
nerve_results = nerve(t_sim=1,postproc_script = "is_recruited")         #Run the simulation

We can plot the nerve again and highlight axons that are recruited:

fig, ax = plt.subplots(1, 1, figsize=(6,6))
nerve_results.plot_recruited_fibers(ax)
ax.set_xlabel("z-axis (µm)")
ax.set_ylabel("y-axis (µm)")

No fiber activated are activated, of course the electrodes are stimulating with the dummy_stimstimulus! Let’s change the stimulus of LIFE_2 (in fascicle_2) with a 100µs-long 60µA monophasic cathodic pulse:

t_start = 0.1       #start of the pulse, in ms
t_pulse = 0.1       #duration of the pulse, in ms
amp_pulse = 60      #amplitude of the pulse, in uA

pulse_stim = nrv.stimulus()
pulse_stim.pulse(t_start, -amp_pulse, t_pulse)      #cathodic pulse

fig, ax = plt.subplots()                            #plot it
pulse_stim.plot(ax) #
ax.set_ylabel("Amplitude (µA)")
ax.set_xlabel("Time (ms)")

We can change the stimulus of LIFE_2 by calling change_stimulus_from_electrode of the nerve-object with the LIFE_2 ID and the new stimulus. We then re-run the simulation and plot the activated fibers.

nerve.change_stimulus_from_electrode(ID_elec=2,stimulus=pulse_stim)
nerve_results = nerve(t_sim=3,postproc_script = "is_recruited")

fig, ax = plt.subplots(1, 1, figsize=(6,6))
nerve_results.plot_recruited_fibers(ax)
ax.set_xlabel("z-axis (µm)")
ax.set_ylabel("y-axis (µm)")

Now we see some activation some fibers being recruited! All myelinated fibers in the fascicle_2 are recruited, as few as a few unmyelinated ones. Some myelinated fibers are also recruited in fascicle_1 but no unmyelinated ones. We can get the ratio of activated fiber in each fascicle using NRV’s built-in methods.

NOTE

Note that FEM is not recomputed between this simulation run and the previous. Indeed, as long as we don’t change any geometrical properties of the model, we only need to run the FEM solver once. This is automatically handled by the framework.

fasc_results = nerve_results.get_fascicle_results(ID = 1)                                #get results in fascicle 1
unmyel = fasc_results.get_recruited_axons('unmyelinated', normalize = True)              #get ratio of unmyelinated axon activated in fascicle 1
myel = fasc_results.get_recruited_axons('myelinated', normalize = True)                  #get ratio of myelinated axon activated in fascicle 1

print(f"Proportion of unmyelinated recruited in fascicle_1: {unmyel*100}%")
print(f"Proportion of myelinated recruited in fascicle_1: {myel*100}%")

fasc_results = nerve_results.get_fascicle_results(ID = 2)                                #get results in fascicle 2
unmyel = fasc_results.get_recruited_axons('unmyelinated', normalize = True)              #get ratio of unmyelinated axon activated in fascicle 2
myel = fasc_results.get_recruited_axons('myelinated', normalize = True)                  #get ratio of myelinated axon activated in fascicle 2

print(f"Proportion of unmyelinated recruited in fascicle_2: {unmyel*100}%")
print(f"Proportion of myelinated recruited in fascicle_2: {myel*100}%")

Let’s remove the stimulation in LIFE_2 and apply it via LIFE_0 instead:

nerve.change_stimulus_from_electrode(ID_elec=0,stimulus=pulse_stim)
nerve.change_stimulus_from_electrode(ID_elec=2,stimulus=dummy_stim)
nerve_results = nerve(t_sim=3,postproc_script = "is_recruited")

Let’s see how many fibers are activated now:

fasc_results = nerve_results.get_fascicle_results(ID = 1)                                #get results in fascicle 1
unmyel = fasc_results.get_recruited_axons('unmyelinated', normalize = True)              #get ratio of unmyelinated axon activated in fascicle 1
myel = fasc_results.get_recruited_axons('myelinated', normalize = True)                  #get ratio of myelinated axon activated in fascicle 1

print(f"Proportion of unmyelinated recruited in fascicle_1: {unmyel*100}%")
print(f"Proportion of myelinated recruited in fascicle_1: {myel*100}%")

fasc_results = nerve_results.get_fascicle_results(ID = 2)                                #get results in fascicle 2
unmyel = fasc_results.get_recruited_axons('unmyelinated', normalize = True)              #get ratio of unmyelinated axon activated in fascicle 2
myel = fasc_results.get_recruited_axons('myelinated', normalize = True)                  #get ratio of myelinated axon activated in fascicle 2

print(f"Proportion of unmyelinated recruited in fascicle_2: {unmyel*100}%")
print(f"Proportion of myelinated recruited in fascicle_2: {myel*100}%")

fig, ax = plt.subplots(figsize=(8, 8))
nerve_results.plot_recruited_fibers(ax)
ax.set_xlabel("z-axis (µm)")
ax.set_ylabel("y-axis (µm)")

We see that the recruitment profile in the fascicles is very different whether we stimulate with one or another electrode. We can analyze it by plotting recruitment curves.

Recruitment curves with LIFEs

To build the recruitment curve of our three electrodes, we are going to create a quick python function get_recruitment_electrodethat take as argument and electrode ID and a numpy array containing the pulse amplitude for the curve. The function return the ratio of myelinated and unmyelinated fibers recruited in each fascicle in python list.

def get_recruitment_electrode(elec_ID:int, amp_vec:np.array, nerve) -> list:

    nerve.verbose = False

    #create empty list to store results
    unmyel_fasc1,myel_fasc1,unmyel_fasc2,myel_fasc2 = ([] for i in range(4))

    #Deactivate unused electrodes
    elec_IDs = [0,1,2]
    unused_elec = [x for x in elec_IDs if elec_ID != x]
    for elec in unused_elec:
        nerve.change_stimulus_from_electrode(ID_elec=elec,stimulus=dummy_stim)

    #Loop throught amp_vec
    print(f"Stimulating nerve with LIFE_{elec_ID}")
    for idx,amp in enumerate(amp_vec):
        amp = np.round(amp,1)                                                       #get the amplitude
        print(f"Pulse amplitude set to {-amp}µA ({idx+1}/{len(amp_vec)})")
        pulse_stim = nrv.stimulus()                                                 #create a new empty stimulus
        pulse_stim.pulse(t_start, -amp, t_pulse)                                    #create a pulse with the new amplitude
        nerve.change_stimulus_from_electrode(ID_elec=elec_ID,stimulus=pulse_stim)    #attach stimulus to selected electrode
        nerve_results = nerve(t_sim=3,postproc_script = "is_recruited")             #run the simulation

        #add results to lists
        fasc_results = nerve_results.get_fascicle_results(ID = 1)
        unmyel_fasc1.append(fasc_results.get_recruited_axons('unmyelinated', normalize = True))
        myel_fasc1.append(fasc_results.get_recruited_axons('myelinated', normalize = True))
        fasc_results = nerve_results.get_fascicle_results(ID = 2)
        unmyel_fasc2.append(fasc_results.get_recruited_axons('unmyelinated', normalize = True))
        myel_fasc2.append(fasc_results.get_recruited_axons('myelinated', normalize = True))
    return(unmyel_fasc1,myel_fasc1,unmyel_fasc2,myel_fasc2)

We use this function to get the recruitment curve of each electrode with the cathodic pulse amplitude varying from 0µA to 150µA, in 20pts.

NOTE

Running this cell takes about 30min on a laptop. This code can be considerably speed-up by using NRV’s built-in parallelization capabilities on HPC.

amp_min = 0             #start at 0µA
amp_max = 100           #ends at 100µA
n_amp = 20              #20pts
amp_vec = np.linspace(amp_min,amp_max,n_amp)
nrv.parameters.set_nmod_ncore(4)            #number of core allocated to fascicle simulations
unmyel_fasc1_LIFE0,myel_fasc1_LIFE0,unmyel_fasc2_LIFE0, myel_fasc2_LIFE0 = get_recruitment_electrode(0,amp_vec,nerve)
unmyel_fasc1_LIFE1,myel_fasc1_LIFE1,unmyel_fasc2_LIFE1, myel_fasc2_LIFE1 = get_recruitment_electrode(1,amp_vec,nerve)
unmyel_fasc1_LIFE2,myel_fasc1_LIFE2,unmyel_fasc2_LIFE2, myel_fasc2_LIFE2 = get_recruitment_electrode(2,amp_vec,nerve)

del nerve, extra_stim #to avoid meshing error, known bug

Now let’s look at the results for myelinated fibers:

c_LIFE_0 = "darkcyan"
c_LIFE_1 = "orangered"
c_LIFE_2 = "seagreen"

fig, (ax1, ax2) = plt.subplots(1, 2)

ax1.plot(amp_vec,myel_fasc1_LIFE0, '-o', lw=2, color= c_LIFE_0, label = 'LIFE_0')
ax1.plot(amp_vec,myel_fasc1_LIFE1, '-o', lw=2, color= c_LIFE_1, label = 'LIFE_1')
ax1.plot(amp_vec,myel_fasc1_LIFE2, '-o', lw=2, color= c_LIFE_2, label = 'LIFE_2')
ax1.set_title("Fascicle 1 - Myelinated")

ax2.plot(amp_vec,myel_fasc2_LIFE0, '-o', lw=2, color= c_LIFE_0, label = 'LIFE_0')
ax2.plot(amp_vec,myel_fasc2_LIFE1, '-o', lw=2, color= c_LIFE_1, label = 'LIFE_1')
ax2.plot(amp_vec,myel_fasc2_LIFE2, '-o', lw=2, color= c_LIFE_2, label = 'LIFE_2')
ax2.set_title("Fascicle 2 - Myelinated")

for ax in ax1, ax2:
    ax.set_xlabel('Amplitude (µA)')
    ax.set_ylabel('Recruitment')
    ax.legend()

fig.tight_layout()

Myelinated fibers are progressively recruited when increasing the pulse amplitude. LIFE_1 recruits the entire fascicle_1 without recruiting any axon in fascicle_2. Oppositely, LIFE_2 recruits the entire fascicle_2 without recruiting any axon in fascicle_1. In other words, intrafascicular selective activation is possible with LIFE_1and LIFE_2. LIFE_0 however, located is neither of the two fascicles, can’t selectively activate one or the other fascicle.

NOTE

Proper curve analysis would require more simulation points. The presented result is for demonstration purposes only.

Let’s plot the unmyelinated fibers’ recruitment curves:

fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(amp_vec,unmyel_fasc1_LIFE0, '-o', lw=2, color= c_LIFE_0, label = 'LIFE_0')
ax1.plot(amp_vec,unmyel_fasc1_LIFE1, '-o', lw=2, color= c_LIFE_1, label = 'LIFE_1')
ax1.plot(amp_vec,unmyel_fasc1_LIFE2, '-o', lw=2, color= c_LIFE_2, label = 'LIFE_2')
ax1.set_title("Fascicle 1 - Unmyelinated")

ax2.plot(amp_vec,unmyel_fasc2_LIFE0, '-o', lw=2, color= c_LIFE_0, label = 'LIFE_0')
ax2.plot(amp_vec,unmyel_fasc2_LIFE1, '-o', lw=2, color= c_LIFE_1, label = 'LIFE_1')
ax2.plot(amp_vec,unmyel_fasc2_LIFE2, '-o', lw=2, color= c_LIFE_2, label = 'LIFE_2')
ax2.set_title("Fascicle 2 - Unmyelinated")

for ax in ax1, ax2:
    ax.set_xlabel('Amplitude (µA)')
    ax.set_ylabel('Recruitment')
    ax.legend()

fig.tight_layout()

Activation of unmyelinated fibers requires much higher pulse amplitude. Electrodes located in the fascicle recruits at most about 10% of the unmyelinated fibers in fascicle_1 and about 70% in fascicle_2. Electrode outside the fascicle or located in the other one fail at recruiting myelinated fibers.

Recruitment curves with a monopolar cuff-like electrode

Let’s create a second nerve with a cuff electrode now:

#creating the fascicles are populating them
fascicle_1_c = nrv.fascicle(diameter=fasc1_d,ID=1)
fascicle_2_c = nrv.fascicle(diameter=fasc2_d, ID=2)
fascicle_1_c.fill_with_population(axons_diameters, axons_type, delta=5, fit_to_size = True)
fascicle_2_c.fill_with_population(axons_diameters, axons_type, delta=5, fit_to_size = True)

#set simulation parameters
fascicle_1_c.set_axons_parameters(unmyelinated_only=True,**u_param)
fascicle_1_c.set_axons_parameters(myelinated_only=True,**m_param)
fascicle_2_c.set_axons_parameters(unmyelinated_only=True,**u_param)
fascicle_2_c.set_axons_parameters(myelinated_only=True,**m_param)

#create the nerve and add fascicles
nerve_cuff = nrv.nerve(length=nerve_l, diameter=nerve_d, Outer_D=outer_d)
nerve_cuff.add_fascicle(fascicle=fascicle_1_c, y=fasc1_y, z=fasc1_z)
nerve_cuff.add_fascicle(fascicle=fascicle_2_c, y=fasc2_y, z=fasc2_z)

#set the simulation flags
nerve_cuff.save_results = False
nerve_cuff.return_parameters_only = False
nerve_cuff.verbose = True

We now create a FEM stimulation context, create a cuff electrode using the CUFF_electrode-class, combine everything and add it to the nerve_cuff-object:

extra_stim_cuff = nrv.FEM_stimulation(endo_mat="endoneurium_ranck",      #endoneurium conductivity
                                 peri_mat="perineurium",            #perineurium conductivity
                                 epi_mat="epineurium",              #epineurium conductivity
                                 ext_mat="saline")                  #saline solution conductivity

contact_length=1000         # length (width) of the cuff contact, in um
contact_thickness=100       # thickness of the contact, in um
insulator_length=1500       # length (width) of the cuff insulator, on top of the contact
insulator_thickness=500     # thickness of the in insulator
x_center = nerve_l/2        # x-position of the cuff

cuff_1 = nrv.CUFF_electrode('CUFF', contact_length=contact_length,
    contact_thickness=contact_thickness, insulator_length=insulator_length,
    insulator_thickness=insulator_thickness, x_center=x_center)

extra_stim_cuff.add_electrode(cuff_1, dummy_stim)
nerve_cuff.attach_extracellular_stimulation(extra_stim_cuff)

fig, ax = plt.subplots(figsize=(8, 8))
nerve_cuff.plot(ax)

We can now simulate a recruitment curve with a cuff just like we did with the LIFE electrodes:

#create empty list to store results
unmyel_fasc1_cuff,myel_fasc1_cuff,unmyel_fasc2_cuff,myel_fasc2_cuff = ([] for i in range(4))

#Loop throught amp_vec
print("Stimulating nerve with CUFF")
for idx,amp in enumerate(amp_vec):
    amp = np.round(amp,1)                                                       #get the amplitude
    print(f"Pulse amplitude set to {-amp}µA ({idx+1}/{len(amp_vec)})")
    pulse_stim = nrv.stimulus()                                                 #create a new empty stimulus
    pulse_stim.pulse(t_start, -amp, t_pulse)                                    #create a pulse with the new amplitude
    nerve_cuff.change_stimulus_from_electrode(ID_elec=0,stimulus=pulse_stim)          #attach stimulus to selected electrode
    nerve_results = nerve_cuff(t_sim=3,postproc_script = "is_recruited")             #run the simulation

    #add results to lists
    fasc_results = nerve_results.get_fascicle_results(ID = 1)
    unmyel_fasc1_cuff.append(fasc_results.get_recruited_axons('unmyelinated', normalize = True))
    myel_fasc1_cuff.append(fasc_results.get_recruited_axons('myelinated', normalize = True))
    fasc_results = nerve_results.get_fascicle_results(ID = 2)
    unmyel_fasc2_cuff.append(fasc_results.get_recruited_axons('unmyelinated', normalize = True))
    myel_fasc2_cuff.append(fasc_results.get_recruited_axons('myelinated', normalize = True))

And plot the results:

c_fascicle_0 = "royalblue"
c_fascicle_1 = "orange"

fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(amp_vec,myel_fasc1_cuff, '-o', lw=2, color= c_fascicle_0, label = 'fascicle_0')
ax1.plot(amp_vec,myel_fasc2_cuff, '-o', lw=2, color= c_fascicle_1, label = 'fascicle_1')
ax1.set_title("Fascicle 1 - Myelinated")

ax2.plot(amp_vec,unmyel_fasc1_cuff, '-o', lw=2, color= c_fascicle_0, label = 'fascicle_0')
ax2.plot(amp_vec,unmyel_fasc2_cuff, '-o', lw=2, color= c_fascicle_1, label = 'fascicle_1')
ax2.set_title("Fascicle 1 - Unmyelinated")

for ax in ax1, ax2:
    ax.set_xlabel('Amplitude (µA)')
    ax.set_ylabel('Recruitment')
    ax.legend()

fig.tight_layout()
#plt.show()