#First, we import different libraries into the script
import pandas as pd #for data handling
import matplotlib.pyplot as plt #for plotting
import numpy as np #for array manipulations

#defining the folder path and file name
path = 'C:\\Users\\Carl Zeiss\\Desktop\\N2a Ca Fura-2\\' #we use \\ to define the folders in the path
fn = 'KCl addition 10to50 in 10mM n2a undiff.csv' #this is the file name
save_folder = 'C:\\Users\\Carl Zeiss\\Desktop\\N2a Ca Fura-2\\' #we define where to save our plots and analysis

#loading the csv file
df = pd.read_csv(path+fn, delimiter = ';', )
c_names = df.columns #this takes all column names
calcium_340_col = [a for a in c_names if 'IntensityMean::340' in a] #this defines all columns containing 340 nm Ca traces
calcium_380_col = [a for a in c_names if 'IntensityMean::380' in a] #this defines all columns containing 380 nm Ca traces

calcium_340 = [np.array(df[a][1:]) for a in calcium_340_col]  #then we data only from the columns we defined above
calcium_380 = [np.array(df[a][1:]) for a in calcium_380_col]
time = df.iloc[1:,0] #taking the time
time = np.array(time).astype(float)

from scipy import signal #for signal processing

# # #first, we detrend the signals; this is useful when the baseline is not set to 0
ca_340 = [signal.detrend(a) for a in calcium_340]
ca_380 = [signal.detrend(a) for a in calcium_380]
# #if we are interested in the initial baseline to be set at 0, we can subtract some average value of the first 10 frames; this is important when the initial recording period is much shorter than the later stationary state;
# #to see how it looks like before subtraction, you can comment out the following line by putting # before it.
ca_380 = [a - np.mean(a[0:10]) for a in ca_380]

fig, ax = plt.subplots(1,2, figsize = (30,10), sharey = True )
for i, (tr1,tr2) in enumerate(zip(ca_340, ca_380)):
    ax[0].plot(time, tr1, linewidth = 3, color = 'green', alpha = 1-i/len(calcium_340), label = '{}'.format(i+1))
    ax[1].plot(time, tr2 , linewidth = 2.5, color = 'red', alpha = 1-i/len(calcium_340))
    ax[0].legend(title = 'Cell #', title_fontsize = 20, loc = 'lower left', prop = {'size': 15})
    
ax[1].set_title('Fura-2 380 nm signal', size = 30)
ax[0].set_title('Fura-2 340 nm signal', size = 30)
ax[0].set_xlabel('Time [s]', size = 30)
ax[0].set_ylabel('dF/F [a.u.]', size = 30)
ax[1].set_xlabel('Time [s]', size = 30)

for a in ax:
    a.tick_params(axis = 'both', which = 'major', labelsize = 24)

ratios = [tr_340/tr_380 for tr_340, tr_380 in zip(calcium_340,calcium_380)]
# #Detrending ratio for clarity
# ratios = [signal.detrend(ratio) for ratio in ratios] # this calculates the baseline from the longest steady segment
# ratios = [ratio - np.mean(ratio[:10]) for ratio in ratios] # because the baseline is shorter before than after, we subtract the mean of the first 10 frames to land at 0
ratios_filt = [signal.savgol_filter(ratio, 21, 3) for ratio in ratios] #additionally, we can filter the curve to more reliably define the maximum
# #Savitzky-Golay filter is polynomial sliding-window filter; here, we filter the curve with 3rd-order polynomial with 21-frame window
# #More on this filter: https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter

# #Plotting the 340/380 ratio
fig, ax = plt.subplots(1,1, figsize = (8,6))
for ratio, ratio_filt in zip(ratios, ratios_filt):
    ax.plot(time, ratio, linewidth = 2, color = 'k', alpha = 0.7, label = 'Raw trace') #plotting the raw, detrended ratio
    ax.plot(time, ratio_filt, linewidth = 1.5, color = 'deeppink', alpha = 0.8, label = 'Filtered trace') #filtered ratio; note the reduction in noise
    ax.set_title ('340/380 ratio', size = 18)
    
ax.set_xlabel('Time [s]', size = 18)
ax.legend(prop = {'size': 15})

# #next, we find the local maximum that defines the peak of Ca increase
# #in the function, we define parameters that describe our peak well:
# # - the height has to be at least 2.3 std above the median (this is the height threshold)
# # - the half-width of the peak has to be at least 100 frames (this is emipirically estimated and can vary)
# #Very similar algortihm can be used for ephys data, to detect action potentials
peak = [signal.find_peaks(ratio_filt, height = np.median(ratio) + 2.3*np.std(ratio), width = 100)[0] for ratio_filt in ratios_filt]

peak_ratio = [ratio_filt[pk] for ratio_filt,pk in zip(ratios, peak)] #taking the peak value of the ratio and storing in peak_ratio variable

for ratio, pk, number in zip(ratios_filt, peak, peak_ratio):
    number = round(number[0],2)
    ax.vlines(x = time[pk], ymin = np.min(ratio), ymax = ratio[pk], linestyle = '--', color = 'k', linewidth = 1.2)
    ax.hlines(xmin = 0, xmax = 300, y = np.min(ratio), linestyle = '--', color = 'k', linewidth = 1.2)
    ax.text(x = time[pk] + 10, y = ratio[pk] - ratio[pk]/4, s = 'R = {}'.format(number), size = 14)

ax.tick_params(axis = 'both', which = 'major', labelsize = 15)
# print(peak_ratio[0]) #printing the peak ratio value

plt.savefig(save_folder + 'Ratio_figure.png', dpi = 300)

kd = 370*(10**(-9)) #in M
Rmin = 0.77
Rmax = 0.996
F380max = 7.02
F380min = 6.99
Q = F380max/F380min#formula: F380max/F380min
R = peak_ratio[0][0]

ca_conc = kd*(R-Rmin)/(Rmax-R)*Q
print('Ca increase: {} nM'.format(ca_conc*1e9))

Ca increase: 66.9832891319402 nM

Evaluating [Ca2+] changes during KCl application in N2a cells using Fura-2¶

We first import all the necessary libraries and load our calcium traces¶

Then, we plot the 340nm and 380nm Fura-2 traces for each cell¶