IKV1 Final Project Report

Assignment¶

The goal of the final project is to learn to construct a more complex network inspired by state-of-the-art research, analyze the results and present them in a written form. Our motivation will be the publication

Kumar A., Cardanobile S., Rotter S. and Aertsen A., The role of inhibition in generating and controlling Parkinson's disease oscillations in the basal ganglia, 2011, Front. Syst. Neurosci. https://doi.org/10.3389/fnsys.2011.00086

In this paper, Kumar et al. explore the origin of beta oscillations involved in Parkinson's disease and how to suppress them with deep brain stimulation (DBS). We will build a smaller network with adjusted parameters that exhibits similar behaviour.

Network construction

Create 2 populations of neurons:

300 neurons (analogue to STN, excitatory)
500 neurons (analogue to GPe, inhibitory)

Neurons in both populations have the same parameters:

Membrane capacity: 30pF

Leak reversal potential: -70 mV

Reset potential: -70 mV

Spike threshold: -60 mV

Refractory period: 2 ms

Membrane time constant: 20ms

Time constant of an excitatory synapse: 1ms Time constant of an inhibitory synapse: 10ms

Connect neurons in both populations (both within and across populations) with the connection rule "fixed indegree". Set the parameters of the connections as follows:

indegree synaptic weight (pA) synaptic delay (ms)

GPe -> STN 30 -15 5

STN -> GPe 50 30 5

GPe -> GPe 30 -5 2

STN -> STN 50 50 2

	indegree	synaptic weight (pA)	synaptic delay (ms)
GPe -> STN	30	-15	5
STN -> GPe	50	30	5
GPe -> GPe	30	-5	2
STN -> STN	50	50	2

Next, create one Poisson generator with rate 800Hz and create connect it to the GPe with weight 40pA and to the STN with weight 65pA.

Create a spike recorder a connect each population to it. If you wish, you can also create and connect a voltmeter to observe the membrane potential.

Now you are ready to simulate the network and observe the activity.

Increased inhibitory input to the GPe

Create a Poisson generator with mean firing rate 800 Hz and connect it to all neurons in the GPe with weight -5 pA. Pick a start and stop time to still see some activity without this input.

Deep brain stimulation

Create a direct current generator with current amplitude -180 pA and connect it to all neurons in the STN. Pick start and stop time so that you see activity in all three conditions (no additional input, additional inhibitory Poisson input to GPe, additional Poisson input to GPe and DBS to STN).

Analysis of results
1. Plot a raster plot of spikes to see how activity changes when you switch on and off the additional external inputs.
2. Plot the distribution of single cell firing rates separately for each population in all three conditions.
3. Plot the power spectrum of the population firing rate in all three conditions. You can choose to present the power spectrum in the form of a spectrogram.
4. Plot the auto-correlation function of the population firing rate for each population and each condition.
5. Test different amplitudes of the DBS. How sensitive is the effect to changes in the DBS amplitude? Create figures to demonstrate your results.
Report
1. Summarize the results from Kumar et al 2011.
2. Present your results, include figures with captions for each point from 4.

Send a pdf to karolina.korvasova@mff.cuni.cz until January 31st, 2024.

Results from Kumar et al 2011¶

The article explores neural mechanisms underlying movement problems in Parkinson's disease, focusing on the role of oscillations in the basal ganglia. Their results and key ideas can be summed up as:

Inhibiting GPe controls basal ganglia oscillations: The research finds that the strength of inhibitory inputs from the striatum to the GPe (globus pallidus external) is a crucial in controlling oscillations within the basal ganglia. These oscillations are associated with Parkinson's movement disorders.

Increased striatal activity triggers oscillation dynamics: The study reveals that the increase in striatal activity is sufficient to trigger oscillations in the basal ganglia. This finding provides a unified explanation for various phenomena: the absence of oscillation in the healthy basal ganglia, oscillations in the dopamine-depleted state (PD), and the suppression of oscillations under DBS.

Better understanding of deep-brain-stimulation (DBS): Revolving around the treatment of PD's symptoms.

Motor impairment understanding: The research also addresses how temporary increases in the activity of striatal neurons can explain motor impairment in patients with PD and the reduced response inhibition observed in patients with DBS implants.

Simulation¶

# onsets in stimulation
start = 0
sick = 750
dbs = 1250
dbs_end = 1750
sick_end = 2250
end = 2500

import nest
import numpy

nest.ResetKernel()

# Set simulation kernel
nest.SetKernelStatus({
  "local_num_threads": 1,
  "resolution": 0.1,
  "rng_seed": 1
})

# Create nodes
stn = nest.Create("iaf_psc_alpha", 300, params={
  "C_m": 30,
  "E_L": -70,
  "V_reset": -70,
  "V_th": -60,
  "t_ref": 2,
  "tau_m": 20,
  "tau_syn_ex": 1,
  "tau_syn_in": 10,
})
gpe = nest.Create("iaf_psc_alpha", 500, params={
  "C_m": 30,
  "E_L": -70,
  "V_reset": -70,
  "V_th": -60,
  "t_ref": 2,
  "tau_m": 20,
  "tau_syn_ex": 1,
  "tau_syn_in": 10,
})
sr1 = nest.Create("spike_recorder", 1)
sr2 = nest.Create("spike_recorder", 1)
pg1 = nest.Create("poisson_generator", 1, params={
  "rate": 800,
  "start": start,
  "stop": end,
})
vm1 = nest.Create("voltmeter", 1)
vm2 = nest.Create("voltmeter", 1)
pg2 = nest.Create("poisson_generator", 1, params={
  "rate": 800,
  "start": sick,
  "stop": sick_end,
})
dc1 = nest.Create("dc_generator", 1, params={
  "amplitude": -180,
  "start": dbs,
  "stop": dbs_end,
})

# Connect nodes
nest.Connect(stn, gpe, conn_spec={
  "rule": "fixed_indegree",
  "indegree": 50,
}, syn_spec={ 
  "weight": 30,
  "delay": 5,
})
nest.Connect(stn, stn, conn_spec={
  "rule": "fixed_indegree",
  "indegree": 50,
}, syn_spec={ 
  "weight": 50,
  "delay": 2,
})
nest.Connect(stn, sr1)
nest.Connect(gpe, sr2)
nest.Connect(gpe, stn, conn_spec={
  "rule": "fixed_indegree",
  "indegree": 30,
}, syn_spec={ 
  "weight": -15,
  "delay": 5,
})
nest.Connect(gpe, gpe, conn_spec={
  "rule": "fixed_indegree",
  "indegree": 30,
}, syn_spec={ 
  "weight": -5,
  "delay": 2,
})
nest.Connect(pg1, gpe, syn_spec={ 
  "weight": 40,
})
nest.Connect(pg1, stn, syn_spec={ 
  "weight": 65,
})
nest.Connect(vm1, gpe)
nest.Connect(vm2, stn)
nest.Connect(pg2, gpe, syn_spec={ 
  "weight": -5,
})
nest.Connect(dc1, stn)

# Run simulation
nest.Simulate(2500)

response = {
  "events": [sr1.events, sr2.events, vm1.events, vm2.events, ]
}

              -- N E S T --
  Copyright (C) 2004 The NEST Initiative

 Version: 3.6.0
 Built: Nov 28 2023 13:46:52

 This program is provided AS IS and comes with
 NO WARRANTY. See the file LICENSE for details.

 Problems or suggestions?
   Visit https://www.nest-simulator.org

 Type 'nest.help()' to find out more about NEST.


Nov 30 16:01:36 SimulationManager::set_status [Info]: 
    Temporal resolution changed from 0.1 to 0.1 ms.

Nov 30 16:01:36 NodeManager::prepare_nodes [Info]: 
    Preparing 807 nodes for simulation.

Nov 30 16:01:36 SimulationManager::start_updating_ [Info]: 
    Number of local nodes: 807
    Simulation time (ms): 2500
    Number of OpenMP threads: 1
    Number of MPI processes: 1

Nov 30 16:01:39 SimulationManager::run [Info]: 
    Simulation finished.

import pickle

with open('response.pkl', 'rb') as file:
    response = pickle.load(file)

Analysis of replicated results¶

%matplotlib inline

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

sns.set_theme()

1. Plot a raster plot of spikes to see how activity changes when you switch on and off the additional external inputs¶

The plot shows spike times of every recorded neuron in each condition. From time 0-750 and 2250-2500 the simulation is considered to capture a healthy interaction (labeled healthy). Interval 750-1250 and 1750-2250 is the one showing spiking activity properties of the Parkinson's disease (labeled sick). And the deep brain stimulation is within the remaining interval: 1750-2250 (labeled DBS).

ids_STN, ids_GPe = response["events"][0]["senders"], response["events"][1]["senders"]
spikes_t_STN, spikes_t_GPe = response["events"][0]["times"], response["events"][1]["times"]

spike_data = []
for neuron_id, spike_t in zip(ids_STN, spikes_t_STN):
    spike_data.append({'Neuron ID': neuron_id, 'Spike Time': spike_t, 'Population': 'STN'})

for neuron_id, spike_t in zip(ids_GPe, spikes_t_GPe):
    spike_data.append({'Neuron ID': neuron_id, 'Spike Time': spike_t, 'Population': 'GPe'})
    
spike_df = pd.DataFrame(spike_data)

plt.figure(figsize=(12, 6))
sns.scatterplot(data=spike_df, x='Spike Time', y='Neuron ID', hue='Population', palette=['mediumspringgreen', 'hotpink'], legend=True, s=5)

def set_time_xtics():
    plt.grid(axis='x')
    plt.xticks([start, sick, dbs, dbs_end, sick_end])

    for pos in [start, sick, dbs, dbs_end, sick_end]:
        plt.axvline(x=pos, color='white', linewidth=2, zorder=0)

set_time_xtics()
plt.title('Raster Plot of Spikes for STN and GPe Populations', fontsize=16)
plt.legend(title='Population', loc='upper right')
plt.xlabel('Time (ms)')
plt.ylabel('Neuron ID')

plt.show()

2. Plot the distribution of single cell firing rates separately for each population in all three conditions¶

To see the distributions of firing rates among the conditions, I've considered two approaches.

Let's plot it simply with averaged firing rates for each neuron in each condition - the bins are the intervals where conditions take place.

healthy_df = spike_df[(spike_df['Spike Time'] <= sick) | (spike_df['Spike Time'] > sick_end)].copy()
sick_df = spike_df[((spike_df['Spike Time'] > sick) & (spike_df['Spike Time'] <= dbs)) |\
                   ((spike_df['Spike Time'] > dbs_end) & (spike_df['Spike Time'] <= sick_end))].copy()
dbs_df = spike_df[(spike_df['Spike Time'] > dbs) & (spike_df['Spike Time'] <= dbs_end)].copy()

healthy_df.loc[:, 'Firing Rate'] = healthy_df.groupby('Neuron ID')['Spike Time'].transform('count') # / 1
sick_df.loc[:, 'Firing Rate'] = sick_df.groupby('Neuron ID')['Spike Time'].transform('count') # / 1
dbs_df.loc[:, 'Firing Rate'] = dbs_df.groupby('Neuron ID')['Spike Time'].transform('count') * 2 # / 0.5

healthy_df.loc[:, 'Status'] = "Healthy"
sick_df.loc[:, 'Status'] = "Sick"
dbs_df.loc[:, 'Status'] = "DBS"

combined_df = pd.concat([healthy_df, sick_df, dbs_df])

plt.figure(figsize=(12, 6))
sns.violinplot(data=combined_df, x='Firing Rate', y='Status', hue='Population', palette=['mediumspringgreen', 'hotpink'], \
               inner='box', dodge=False)

plt.title('Distributions of Firing Rates for STN and GPe Populations in each Condition', fontsize=16)
plt.xlabel('Firing Rate (Hz)')
plt.ylabel('Condition')
plt.show()

Second approach involves extracting and plotting smooth firing rate. I've used Alpha function to extract it.

import numpy as np

tau = 25.0  # Time constant for alpha function in ms
alpha = 1/tau # Rate constant for alpha function in ms^-1
dt = 1.0  # Time step for continuous representation in ms
total_duration = 2500 # Total duration of simulation in ms

def smooth_spike_train(spike_times, alpha=alpha, total_duration=total_duration, dt=dt):
    times = np.arange(0, total_duration, dt)
    firing_rate = np.zeros_like(times)
    alpha_function = (alpha ** 2) * times * np.exp(-alpha * times) * (times >= 0)

    alpha_function = alpha_function / np.sum(alpha_function)

    # vectorized version
    for spike in spike_times:
        deltas = times - spike
        valid_mask = (deltas >= 0) & (deltas < len(alpha_function) * dt)
        firing_rate[valid_mask] += alpha_function[(deltas[valid_mask] // dt).astype(int)]

    firing_rate *= (1000 / dt)

    return firing_rate


def get_population_average_firing_rate(df, population):
    unique_neurons = df[df['Population'] == population]['Neuron ID'].unique()
    average_firing_rate = np.zeros(int(total_duration / dt))

    for neuron in unique_neurons:
        neuron_spike_times = df[(df['Neuron ID'] == neuron) & (df['Population'] == population)]['Spike Time']
        smoothed_fr = smooth_spike_train(neuron_spike_times)
        average_firing_rate += smoothed_fr

    average_firing_rate /= len(unique_neurons)
    
    return average_firing_rate

# get smooth firing rates averaged over population
smooth_stn_fr = get_population_average_firing_rate(spike_df, 'STN')
smooth_gpe_fr = get_population_average_firing_rate(spike_df, 'GPe')

plt.figure(figsize=(12, 6))
plt.plot(np.arange(0, total_duration, dt), smooth_stn_fr, label='STN', color='mediumspringgreen')
plt.plot(np.arange(0, total_duration, dt), smooth_gpe_fr, label='GPe', color='hotpink')

plt.legend(title='Population', loc='upper right')
plt.title('Smooth Firing Rates for STN and GPe Populations', fontsize=16)

set_time_xtics()
plt.xlabel('Time (ms)')
plt.ylabel('Firing Rate (Hz)')
plt.show()

3. Plot the power spectrum of the population firing rate in all three conditions. You can choose to present the power spectrum in the form of a spectrogram.¶

To compute it I've extracted smooth firing rates - as done above but they are not averaged - for each condition and used the Welch's method to extract the final power spectra. This shows us how the firing rates are distributed across frequencies.

from scipy.signal import welch

# Smooth the spike train for each population in each condition
smooth_stn_fr_1 = smooth_spike_train(healthy_df[healthy_df['Population'] == 'STN']['Spike Time'])
smooth_stn_fr_2 = smooth_spike_train(sick_df[sick_df['Population'] == 'STN']['Spike Time'])
smooth_stn_fr_3 = smooth_spike_train(dbs_df[dbs_df['Population'] == 'STN']['Spike Time'])

smooth_gpe_fr_1 = smooth_spike_train(healthy_df[healthy_df['Population'] == 'GPe']['Spike Time'])
smooth_gpe_fr_2 = smooth_spike_train(sick_df[sick_df['Population'] == 'GPe']['Spike Time'])
smooth_gpe_fr_3 = smooth_spike_train(dbs_df[dbs_df['Population'] == 'GPe']['Spike Time'])

# Compute power spectra using Welch's method
freqs_stn_1, pow_spec_stn_1 = welch(smooth_stn_fr_1, fs=1000/dt, nperseg=len(smooth_stn_fr_1))
freqs_stn_2, pow_spec_stn_2 = welch(smooth_stn_fr_2, fs=1000/dt, nperseg=len(smooth_stn_fr_2))
freqs_stn_3, pow_spec_stn_3 = welch(smooth_stn_fr_3, fs=1000/dt, nperseg=len(smooth_stn_fr_3))

freqs_gpe_1, pow_spec_gpe_1 = welch(smooth_gpe_fr_1, fs=1000/dt, nperseg=len(smooth_gpe_fr_1))
freqs_gpe_2, pow_spec_gpe_2 = welch(smooth_gpe_fr_2, fs=1000/dt, nperseg=len(smooth_gpe_fr_2))
freqs_gpe_3, pow_spec_gpe_3 = welch(smooth_gpe_fr_3, fs=1000/dt, nperseg=len(smooth_gpe_fr_3))

# Function to plot power spectrum
def plot_power_spectrum(freqs, powers, label, ax):
    ax.plot(freqs, powers, label=label)
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Power')
    ax.set_yscale('log')
    ax.set_ylim(10 ** (-8), 10 ** 10)

# Set up plot
fig, axs = plt.subplots(2, 1, figsize=(12, 6))
fig.subplots_adjust(hspace=0.3)
axs[0].set_title('STN Power Spectra', fontsize=14, bbox={'facecolor': 'mediumspringgreen', 'alpha': 0.5, 'boxstyle': 'round,pad=0.3'})
axs[1].set_title('GPe Power Spectra', fontsize=14, bbox={'facecolor': 'hotpink', 'alpha': 0.5, 'boxstyle': 'round,pad=0.3'})

# Plot power spectra for STN
plot_power_spectrum(freqs_stn_1, pow_spec_stn_1, 'Healthy', axs[0])
plot_power_spectrum(freqs_stn_2, pow_spec_stn_2, 'Sick', axs[0])
plot_power_spectrum(freqs_stn_3, pow_spec_stn_3, 'DBS', axs[0])
axs[0].set_xticklabels([])
axs[0].set_xlabel('')

# Plot power spectra for GPe
plot_power_spectrum(freqs_gpe_1, pow_spec_gpe_1, 'Healthy', axs[1])
plot_power_spectrum(freqs_gpe_2, pow_spec_gpe_2, 'Sick', axs[1])
plot_power_spectrum(freqs_gpe_3, pow_spec_gpe_3, 'DBS', axs[1])

axs[0].legend(title='Condition', loc='upper right')
plt.tight_layout()
plt.show()

4. Plot the auto-correlation function of the population firing rate for each population and each condition.¶

Plotting the auto-correlation of each population in each condition should tell us something about the periodicity and relationship between smoothed firing rates.

import numpy as np
import matplotlib.pyplot as plt

def compute_autocorrelation(firing_rate):
    autocorr = np.correlate(firing_rate, firing_rate, mode='full')
    # Normalize
    autocorr /= autocorr.max()
    # Return only the second half of the auto-correlation (non-negative lags)
    return autocorr[len(autocorr)//2:]

# Compute the auto-correlation for each condition
autocorr_stn_1 = compute_autocorrelation(smooth_stn_fr_1)
autocorr_stn_2 = compute_autocorrelation(smooth_stn_fr_2)
autocorr_stn_3 = compute_autocorrelation(smooth_stn_fr_3)

autocorr_gpe_1 = compute_autocorrelation(smooth_gpe_fr_1)
autocorr_gpe_2 = compute_autocorrelation(smooth_gpe_fr_2)
autocorr_gpe_3 = compute_autocorrelation(smooth_gpe_fr_3)

# Function to plot power spectrum
def plot_auto_corr_population(lags, autocorrelations, labels, ax):
    for autocorr, label in zip(autocorrelations, labels):
        ax.plot(lags, autocorr, label=label)
    ax.set_xlabel('Lag (ms)')
    ax.set_ylabel('Auto-correlation')
    ax.grid(axis='x')
    for pos in [start, sick, dbs, dbs_end, sick_end]:
        ax.axvline(x=pos, color='white', linewidth=2, zorder=0)

fig, axs = plt.subplots(2, 1, figsize=(12, 6))
fig.subplots_adjust(hspace=0.3)
axs[0].set_title('STN Auto-correlation', fontsize=14, bbox={'facecolor': 'mediumspringgreen', 'alpha': 0.5, 'boxstyle': 'round,pad=0.3'})
axs[1].set_title('GPe Auto-correlation', fontsize=14, bbox={'facecolor': 'hotpink', 'alpha': 0.5, 'boxstyle': 'round,pad=0.3'})

lags = np.arange(len(autocorr_stn_1))
conditions = ['Healthy', 'Sick', 'DBS']

plot_auto_corr_population(lags, [autocorr_stn_1, autocorr_stn_2, autocorr_stn_3], conditions, axs[0])
plot_auto_corr_population(lags, [autocorr_gpe_1, autocorr_gpe_2, autocorr_gpe_3], conditions, axs[1])

axs[0].legend(title='Condition', loc='upper right')
axs[0].set_xticklabels([])
axs[0].set_xlabel('')
plt.show()

5. Test different amplitudes of the DBS. How sensitive is the effect to changes in the DBS amplitude? Create figures to demonstrate your results.¶

To see how sensitive the change for different amplitudes is, I've tested it for 21 equally-spaced values. Value -180 pA is the default, simulations were done for values ranging from -280 to -80 (100 pA to each side from the default). Then, we simply took the firing rates of each population and plotted it with a lineplot to symbolize the change. Standard deviation is shown as well. To get a clearer view of what happens with slight nudges around the default, more precise simulations took place - they are shown in the rectangle with black borders - as a zoomed window.

import nest

def simulate_w_amp(amplitude=-180):
    nest.ResetKernel()

    # Set simulation kernel
    nest.SetKernelStatus({
      "local_num_threads": 1,
      "resolution": 0.1,
      "rng_seed": 1
    })

    # Create nodes
    stn = nest.Create("iaf_psc_alpha", 300, params={
      "C_m": 30,
      "E_L": -70,
      "V_reset": -70,
      "V_th": -60,
      "t_ref": 2,
      "tau_m": 20,
      "tau_syn_ex": 1,
      "tau_syn_in": 10,
    })
    gpe = nest.Create("iaf_psc_alpha", 500, params={
      "C_m": 30,
      "E_L": -70,
      "V_reset": -70,
      "V_th": -60,
      "t_ref": 2,
      "tau_m": 20,
      "tau_syn_ex": 1,
      "tau_syn_in": 10,
    })
    sr1 = nest.Create("spike_recorder", 1)
    sr2 = nest.Create("spike_recorder", 1)
    pg1 = nest.Create("poisson_generator", 1, params={
      "rate": 800,
      "start": 0,
      "stop": 1500,
    })
    vm1 = nest.Create("voltmeter", 1)
    vm2 = nest.Create("voltmeter", 1)
    pg2 = nest.Create("poisson_generator", 1, params={
      "rate": 800,
      "start": 0,
      "stop": 1500,
    })
    dc1 = nest.Create("dc_generator", 1, params={
      "amplitude": amplitude,
      "start": 500,
      "stop": 1500,
    })

    # Connect nodes
    nest.Connect(stn, gpe, conn_spec={
      "rule": "fixed_indegree",
      "indegree": 50,
    }, syn_spec={ 
      "weight": 30,
      "delay": 5,
    })
    nest.Connect(stn, stn, conn_spec={
      "rule": "fixed_indegree",
      "indegree": 50,
    }, syn_spec={ 
      "weight": 50,
      "delay": 2,
    })
    nest.Connect(stn, sr1)
    nest.Connect(gpe, sr2)
    nest.Connect(gpe, stn, conn_spec={
      "rule": "fixed_indegree",
      "indegree": 30,
    }, syn_spec={ 
      "weight": -15,
      "delay": 5,
    })
    nest.Connect(gpe, gpe, conn_spec={
      "rule": "fixed_indegree",
      "indegree": 30,
    }, syn_spec={ 
      "weight": -5,
      "delay": 2,
    })
    nest.Connect(pg1, gpe, syn_spec={ 
      "weight": 40,
    })
    nest.Connect(pg1, stn, syn_spec={ 
      "weight": 65,
    })
    nest.Connect(vm1, gpe)
    nest.Connect(vm2, stn)
    nest.Connect(pg2, gpe, syn_spec={ 
      "weight": -5,
    })
    nest.Connect(dc1, stn)

    # Run simulation
    nest.Simulate(2500)

    response = {
      "events": [sr1.events, sr2.events, vm1.events, vm2.events, ]
    }

    return response

def get_data(response, amp, data):
    events = response["events"]

    ids_STN, ids_GPe = events[0]["senders"], events[1]["senders"]
    spikes_t_STN, spikes_t_GPe = events[0]["times"], events[1]["times"]

    for neuron_id, spike_t in zip(ids_STN, spikes_t_STN):
        data.append({'Neuron ID': f'{neuron_id}_{-amp}', 'Spike Time': spike_t, 'Population': 'STN', 'Amplitude': amp})

    for neuron_id, spike_t in zip(ids_GPe, spikes_t_GPe):
        data.append({'Neuron ID': f'{neuron_id}_{-amp}', 'Spike Time': spike_t, 'Population': 'GPe', 'Amplitude': amp})

def get_firing_rates_of_dbs(df):
    df = df[(df['Spike Time'] > 500) & (df['Spike Time'] <= 1500)]
    df.loc[:, 'Firing Rate'] = df.groupby('Neuron ID')['Spike Time'].transform('count')
    return df

tested_amps = np.linspace(-280, -80, 21)

data = []

for amp in tested_amps:
    get_data(simulate_w_amp(amp), amp, data)

df = pd.DataFrame(data)

with open('different_21_amps_df.pkl', 'rb') as file:
    df = pickle.load(file)

fr_df = get_firing_rates_of_dbs(df)

with open('different_amps_zoom_df.pkl', 'rb') as file:
    df = pickle.load(file)

fr_zoom_df = get_firing_rates_of_dbs(df)

from mpl_toolkits.axes_grid1.inset_locator import inset_axes

fig, ax = plt.subplots(figsize=(12, 6))
sns.lineplot(data=fr_df, x='Amplitude', y='Firing Rate', hue='Population', palette=['hotpink', 'mediumspringgreen'], legend=True,\
             marker='X', linestyle='--', markerfacecolor='black', err_style='band', errorbar='sd', ax=ax)

ax_inset = inset_axes(ax, width="30%", height="30%", bbox_to_anchor=(-0.55, -0.32, 1, 1), bbox_transform=ax.transAxes)
sns.lineplot(data=fr_zoom_df, x='Amplitude', y='Firing Rate', hue='Population', palette=['mediumspringgreen', 'hotpink'],\
             marker='X', linestyle='--', markerfacecolor='black', err_style='band', errorbar='sd', ax=ax_inset)

ax_inset.set_xlabel('')
ax_inset.set_ylabel('')
ax_inset.get_legend().remove()
for spine in ax_inset.spines.values():
    spine.set_edgecolor('black')

ax.set_title('Firing Rates of STN and GPe Populations with Different Amplitudes of DBS', fontsize=16)
ax.set_xlabel('Amplitude (pA)')
ax.set_ylabel('Firing Rate (Hz)')
ax.set_xticks([x for x in tested_amps if x % 20 == 0])
plt.show()