From 9685a9b310061ce1c2642fbfa723cea40742fba2 Mon Sep 17 00:00:00 2001 From: BrianNGitahi <102825698+BrianNGitahi@users.noreply.github.com> Date: Wed, 11 Oct 2023 18:28:40 -0700 Subject: [PATCH] no artifact --- Bleach_compare.py | 72 +++++++ Bleaching Analysis.py | 50 ++++- Bleaching w:out_bline.py | 178 +++++++++++++++++ Bleaching_bline.py | 82 +++++++- Frate.py | 63 ++++++ Functions.py | 46 ----- __pycache__/Functions.cpython-311.pyc | Bin 5115 -> 5115 bytes __pycache__/b_factory.cpython-311.pyc | Bin 0 -> 4039 bytes __pycache__/b_functions.cpython-311.pyc | Bin 0 -> 7584 bytes __pycache__/s_functions.cpython-311.pyc | Bin 0 -> 3452 bytes b_factory.py | 251 +++++++++++++++++++++++ b_functions.py | 253 ++++++++++++++++++++++++ f_rate.py | 67 +++++++ s_functions.py | 180 +++++++++++++++++ 14 files changed, 1182 insertions(+), 60 deletions(-) create mode 100644 Bleach_compare.py create mode 100644 Bleaching w:out_bline.py create mode 100644 Frate.py create mode 100644 __pycache__/b_factory.cpython-311.pyc create mode 100644 __pycache__/b_functions.cpython-311.pyc create mode 100644 __pycache__/s_functions.cpython-311.pyc create mode 100644 b_factory.py create mode 100644 b_functions.py create mode 100644 f_rate.py create mode 100644 s_functions.py diff --git a/Bleach_compare.py b/Bleach_compare.py new file mode 100644 index 0000000..14ff5c4 --- /dev/null +++ b/Bleach_compare.py @@ -0,0 +1,72 @@ +# This script is used to compare the bleaching effects established in the +# Bleaching Analysis and bleaching_bline scripts + + +# import necessary packages +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import pandas as pd +import scipy as sp + +# import functions from the simulation and bleaching libraries +from s_functions import simulate_neuron, simulate_nm_conc, simulate_flourescence_signal +from b_factory import bleach_1, bleach_2, bleach_3, bleach_4 + + +# ANALYSIS 1: comparing the different contributions at one timescale for the bleach factor + +# choose the timescale for the bleach factor +chosen_tau = 10e4 + +# generate an input nm conc array from firing neuron +firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# plot bleached signal for the bleach factor acting on different components of f + +# create timesteps array for the plot +t = np.linspace(0,nm_conc.size-1,nm_conc.size) + +b1 = bleach_1(K_D = 1000, tau=chosen_tau, F_max = 45, F_min = 10, nm_conc=nm_conc, bline_len=5000) +b2 = bleach_2(K_D = 1000, tau=chosen_tau, F_max = 45, F_min = 10, nm_conc=nm_conc, bline_len=5000) +b3 = bleach_3(K_D = 1000, tau=chosen_tau, F_max = 45, F_min = 10, nm_conc=nm_conc, bline_len=5000) +b4 = bleach_4(K_D = 1000, tau=chosen_tau, F_max = 45, F_min = 10, nm_conc=nm_conc, bline_len=5000) + + +# create a fit for the initial values of the bleach 2 and 4 result and subtract it -- revise this later +# it should be a fit of the f value +poly2 = np.polyfit(t,b2,2) +fit2 = np.polyval(poly2,t) + +poly4 = np.polyfit(t,b4,2) +fit4 = np.polyval(poly4,t) + + + +plt.plot(t,b1,label='bleach 1') +#plt.plot(t,b2,label='bleach 2') +# plt.plot(fit2, label='fit2') +plt.plot(b2-fit2,label='subtracted 2') + + +plt.plot(t,b3,label='bleach 3') +#plt.plot(t,b4,label='bleach 4') +plt.plot(b4-fit4,label='subtracted 4') + +plt.xlabel('time(ms)') +plt.ylabel('Delta F/F0') +plt.title('Flourescence intensity signal over time') +plt.legend() + +plt.show() + + +# ANALYSIS 2: + + + + +# ANALYSIS 3: +# check this bleaching effect 1 for different values of tau -- different bleach factor +#different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) diff --git a/Bleaching Analysis.py b/Bleaching Analysis.py index 60feb3a..cc22f8d 100644 --- a/Bleaching Analysis.py +++ b/Bleaching Analysis.py @@ -25,7 +25,7 @@ def bleach_1(K_D, tau, F_max, F_min, nm_conc): bleach = np.exp(-t/tau) # calculate bleached f values: derived in part from eq 2 in Neher/Augustine - f_b = bg_tissue + (K_D*F_min + bleach*nm_conc*F_max)/(K_D + bleach*nm_conc) + f_b = bg_tissue + (K_D*F_min + bleach*nm_conc*F_max)/(K_D + nm_conc) # calculate normalized signal: (assume f0 is the initial f value) delta_ft_f0 = (f_b-f_b[0])/f_b[0] @@ -37,12 +37,6 @@ def bleach_1(K_D, tau, F_max, F_min, nm_conc): plt.ylabel('Delta F/F0') plt.title('Flourescence intensity signal over time (bleach 1)') plt.legend() - #plt.show() - - # print the bleach factor - # plt.plot(bleach) - # plt.xlabel('timesteps(ms)') - # plt.title('Bleach factor as function of time') return delta_ft_f0 @@ -137,4 +131,44 @@ def bleach_3(K_D, tau, F_max, F_min, nm_conc): # plot bleached signal with different time constants for the bleach factor for i in range(len(different_taus)): bleach_3(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) -plt.show() \ No newline at end of file +plt.show() + + +# ANALYSIS 4: bleaching factor acting on all the contributions +def bleach_4(K_D, tau, F_max, F_min, nm_conc): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f_b= bleach*(bg_tissue + (K_D*F_min + nm_conc*F_max)/(K_D + nm_conc)) + + # calculate normalized signal: (assume f0 is the initial f value) + delta_ft_f0 = (f_b-f_b[0])/f_b[0] + + # plot the normalized signal delta f/f0 at the different t + plt.plot(t,delta_ft_f0, label = tau) + plt.xlabel('time(ms)') + plt.ylabel('Delta F/F0') + plt.title('Flourescence intensity signal over time (bleach 4)') + plt.legend() + + return delta_ft_f0 + + +# check this bleaching effect 1 for different values of tau -- different bleach factor +different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# generate an input nm conc array from firing neuron +firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# plot bleached signal with different time constants for the bleach factor +for i in range(len(different_taus)): + bleach_4(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +plt.show() diff --git a/Bleaching w:out_bline.py b/Bleaching w:out_bline.py new file mode 100644 index 0000000..f9760e0 --- /dev/null +++ b/Bleaching w:out_bline.py @@ -0,0 +1,178 @@ +# This script was the original analysis script and +# has since been updated -- cf b_compare or analysis + +# it has the bleaching analysis with the df/f computed without the +# moving baseline calculation of f0 + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import pandas as pd +import scipy as sp + +# import the necessary functions +from s_functions import simulate_neuron, simulate_nm_conc + + +# define the background fluorescence due to the tissue +bg_tissue = 1.5 + +# ANALYSIS 1: bleaching factor acting on the [NM] contribution to F +def bleach_1(K_D, tau, F_max, F_min, nm_conc): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f_b = bg_tissue + (K_D*F_min + bleach*nm_conc*F_max)/(K_D + nm_conc) + + # calculate normalized signal: (assume f0 is the initial f value) + delta_ft_f0 = (f_b-f_b[0])/f_b[0] + + + # plot the normalized signal delta f/f0 at the different t + plt.plot(t,delta_ft_f0, label = tau) + plt.xlabel('time(ms)') + plt.ylabel('Delta F/F0') + plt.title('Flourescence intensity signal over time (bleach 1)') + plt.legend() + + return delta_ft_f0 + + +# check this bleaching effect 1 for different values of tau -- different bleach factor +different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# generate an input nm conc array from firing neuron +firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# plot bleached signal with different time constants for the bleach factor +for i in range(len(different_taus)): + bleach_1(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +plt.show() + + +# ANALYSIS 2: bleaching factor acting on the contributions of the dye, F0 and [NM] +def bleach_2(K_D, tau, F_max, F_min, nm_conc): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f_b= bg_tissue+ bleach*(K_D*F_min + nm_conc*F_max)/(K_D + nm_conc) + + # calculate normalized signal: (assume f0 is the initial f value) + delta_ft_f0 = (f_b-f_b[0])/f_b[0] + + # plot the normalized signal delta f/f0 at the different t + plt.plot(t,delta_ft_f0, label = tau) + plt.xlabel('time(ms)') + plt.ylabel('Delta F/F0') + plt.title('Flourescence intensity signal over time (bleach 2)') + plt.legend() + + return delta_ft_f0 + + +# check this bleaching effect 1 for different values of tau -- different bleach factor +different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# generate an input nm conc array from firing neuron +firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# plot bleached signal with different time constants for the bleach factor +for i in range(len(different_taus)): + bleach_2(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +plt.show() + + +# ANALYIS 3: bleaching factor acting on the background fluorescence +def bleach_3(K_D, tau, F_max, F_min, nm_conc): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f_b= bleach*bg_tissue + (K_D*F_min + nm_conc*F_max)/(K_D + nm_conc) + + # calculate normalized signal: (assume f0 is the initial f value) + delta_ft_f0 = (f_b-f_b[0])/f_b[0] + + # plot the normalized signal delta f/f0 at the different t + plt.plot(t,delta_ft_f0, label = tau) + plt.xlabel('time(ms)') + plt.ylabel('Delta F/F0') + plt.title('Flourescence intensity signal over time (bleach 3)') + plt.legend() + + return delta_ft_f0 + + +# check this bleaching effect 1 for different values of tau -- different bleach factor +different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# generate an input nm conc array from firing neuron +firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# plot bleached signal with different time constants for the bleach factor +for i in range(len(different_taus)): + bleach_3(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +plt.show() + + +# ANALYSIS 4: bleaching factor acting on all the contributions +def bleach_4(K_D, tau, F_max, F_min, nm_conc): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f_b= bleach*(bg_tissue + (K_D*F_min + nm_conc*F_max)/(K_D + nm_conc)) + + # calculate normalized signal: (assume f0 is the initial f value) + delta_ft_f0 = (f_b-f_b[0])/f_b[0] + + # plot the normalized signal delta f/f0 at the different t + plt.plot(t,delta_ft_f0, label = tau) + plt.xlabel('time(ms)') + plt.ylabel('Delta F/F0') + plt.title('Flourescence intensity signal over time (bleach 4)') + plt.legend() + + return delta_ft_f0 + + +# check this bleaching effect 1 for different values of tau -- different bleach factor +different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# generate an input nm conc array from firing neuron +firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# plot bleached signal with different time constants for the bleach factor +for i in range(len(different_taus)): + bleach_4(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +plt.show() diff --git a/Bleaching_bline.py b/Bleaching_bline.py index 7f9d4c5..0e68d0c 100644 --- a/Bleaching_bline.py +++ b/Bleaching_bline.py @@ -28,7 +28,7 @@ def bleach_1(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): bleach = np.exp(-t/tau) # calculate bleached f values: derived in part from eq 2 in Neher/Augustine - f = bg_tissue + (K_D*F_min + bleach*nm_conc*F_max)/(K_D + bleach*nm_conc) + f = bg_tissue + (K_D*F_min + bleach*nm_conc*F_max)/(K_D + nm_conc) # define the signal array delta_ft_f0 = np.zeros(f.size) @@ -50,7 +50,7 @@ def bleach_1(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): # plot the normalized signal delta f/f0 at the different t - plt.plot(t,delta_ft_f0, label = tau) + plt.plot(t,delta_ft_f0, label = tau + 1) plt.xlabel('time(ms)') plt.ylabel('Delta F/F0') plt.title('Flourescence intensity signal over time (bleach 1)') @@ -72,6 +72,7 @@ def bleach_1(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): plt.show() + # ANALYSIS 2: bleaching factor acting on the contributions of the dye, F0 and [NM] def bleach_2(K_D, tau, F_max, F_min, nm_conc, bline_len =5000): @@ -105,7 +106,7 @@ def bleach_2(K_D, tau, F_max, F_min, nm_conc, bline_len =5000): delta_ft_f0[i] = (f[i]-f0)/(f0) # plot the normalized signal delta f/f0 at the different t - plt.plot(t,delta_ft_f0, label = tau) + plt.plot(t,delta_ft_f0, label = tau + 2) plt.xlabel('time(ms)') plt.ylabel('Delta F/F0') plt.title('Flourescence intensity signal over time (bleach 2)') @@ -127,8 +128,8 @@ def bleach_2(K_D, tau, F_max, F_min, nm_conc, bline_len =5000): plt.show() -# ANALYIS 3: bleaching factor acting on the background fluores -# cence + +# ANALYSIS 3: bleaching factor acting on the background fluorescence def bleach_3(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): # create timesteps array for the plot @@ -161,7 +162,7 @@ def bleach_3(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): delta_ft_f0[i] = (f[i]-f0)/(f0) # plot the normalized signal delta f/f0 at the different t - plt.plot(t,delta_ft_f0, label = tau) + plt.plot(t,delta_ft_f0, label = tau + 3) plt.xlabel('time(ms)') plt.ylabel('Delta F/F0') plt.title('Flourescence intensity signal over time (bleach 3)') @@ -180,4 +181,73 @@ def bleach_3(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): # plot bleached signal with different time constants for the bleach factor for i in range(len(different_taus)): bleach_3(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +plt.show() + + + + +# ANALYSIS 4: bleaching on all contributions +def bleach_4(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f= bleach*(bg_tissue + (K_D*F_min + nm_conc*F_max)/(K_D + nm_conc)) + + # define the signal array + delta_ft_f0 = np.zeros(f.size) + + # calculate f0 values and populate the signal array + for i in range(f.size): + + # calculate f0 by averaging the previous x number of f values + # if x is bigger than the current index then use all the prev f values + # where x is the length of the moving baseline + + if i < bline_len: + f0 = np.average(f[:i]) + else: + f0 = np.average(f[i-bline_len:i]) + + # calculate normalized signal using the calculated f0 + delta_ft_f0[i] = (f[i]-f0)/(f0) + + # plot the normalized signal delta f/f0 at the different t + plt.plot(t,delta_ft_f0, label = tau + 4) + plt.xlabel('time(ms)') + plt.ylabel('Delta F/F0') + plt.title('Flourescence intensity signal over time (bleach 4)') + plt.legend() + + return delta_ft_f0 + + +# check this bleaching effect 1 for different values of tau -- different bleach factor +different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# generate an input nm conc array from firing neuron +firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# plot bleached signal with different time constants for the bleach factor +for i in range(len(different_taus)): + bleach_4(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +plt.show() + + +# ANALYSIS 5: comparing the different contributions at one timescale for the bleach factor + +# choose the timescale for the bleach factor +chosen_tau = 10e4 + +bleach_1(K_D=1000,tau=chosen_tau, F_max=45, F_min=10, nm_conc=nm_conc) +bleach_2(K_D=1000,tau=chosen_tau, F_max=45, F_min=10, nm_conc=nm_conc) +bleach_3(K_D=1000,tau=chosen_tau, F_max=45, F_min=10, nm_conc=nm_conc) +bleach_4(K_D=1000,tau=chosen_tau, F_max=45, F_min=10, nm_conc=nm_conc) plt.show() \ No newline at end of file diff --git a/Frate.py b/Frate.py new file mode 100644 index 0000000..808a04b --- /dev/null +++ b/Frate.py @@ -0,0 +1,63 @@ +# This script is used for analysis of the effect of +# changing the firing rate of the neuron on the measured signal + +# import necessary pacakges +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import pandas as pd +import scipy as sp + +from Functions import simulate_neuron, simulate_nm_conc, simulate_flourescence_signal + + +# ANALYSIS 1: + +# Define function that makes generates several firing neurons with different firing rates +# and then compares the average df/f signal across them + +def signal_vs_activity(number_neurons,firing_rates): + + # create array to store the average signals generated by the neurons + average_signals = [] + + # for each firing rate, simulate a neuron, the nm dynamics from it and the signal + # then get the average signal + for i in range(firing_rates.size): + + # simulate neuron with different firing rate + guinea_neuron = simulate_neuron(70000,firing_rates[i]) + + # generate the nm_conc + guinea_nm_conc, guinea_b_conc, guinea_c_conc = simulate_nm_conc(guinea_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + + # then generate the signal + guinea_signal = simulate_flourescence_signal(K_D = 1000, F_max = 45, F_min = 10, nm_conc=guinea_nm_conc) + + # get the average of this signal and add it the original array + average_signals.append(np.average(guinea_signal)) + + print('simulated {} neuron(s)'.format(i+1)) + + + # plot the different firing rates vs their signals + plt.scatter(firing_rates,average_signals) + plt.xlabel('Firing rates(Hz)') + plt.ylabel('Average df/f signal') + plt.title('Signal vs activity plot') + plt.show() + + # plot at log scale + plt.loglog(firing_rates,average_signals, 'o') + plt.xlabel('Firing rates(Hz)') + plt.ylabel('Average df/f signal') + plt.title('Signal vs activity plot -Log') + plt.show() + + +# Check 1: +signal_vs_activity(2,np.array([1,2,4,10,50,70,100,200,300,400,500,600,700,800,900,1000,1500,2000,2500,3000,5000])) + +# results: looks plausible -- tho at a big scale it has a linear/exponential trend then it levels off +# -- it's linear + diff --git a/Functions.py b/Functions.py index 6317ec6..d093251 100644 --- a/Functions.py +++ b/Functions.py @@ -171,56 +171,10 @@ def simulate_flourescence_signal(K_D, F_max, F_min, nm_conc): #f_signal = simulate_flourescence_signal(K_D = 1000, F_max = 45, F_min = 10, nm_conc=nm_conc) -# ANALYSIS 1: -# Define function that makes generates several firing neurons with different firing rates -# and then compares the average df/f signal across them - -def signal_vs_activity(number_neurons,firing_rates): - - # create array to store the average signals generated by the neurons - average_signals = [] - - # for each firing rate, simulate a neuron, the nm dynamics from it and the signal - # then get the average signal - for i in range(firing_rates.size): - - # simulate neuron with different firing rate - guinea_neuron = simulate_neuron(70000,firing_rates[i]) - - # generate the nm_conc - guinea_nm_conc, guinea_b_conc, guinea_c_conc = simulate_nm_conc(guinea_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) - - # then generate the signal - guinea_signal = simulate_flourescence_signal(K_D = 1000, F_max = 45, F_min = 10, nm_conc=guinea_nm_conc) - - # get the average of this signal and add it the original array - average_signals.append(np.average(guinea_signal)) - - print('simulated {} neuron(s)'.format(i+1)) - - - # plot the different firing rates vs their signals - plt.scatter(firing_rates,average_signals) - plt.xlabel('Firing rates(Hz)') - plt.ylabel('Average df/f signal') - plt.title('Signal vs activity plot') - plt.show() - - # plot at log scale - plt.loglog(firing_rates,average_signals, 'o') - plt.xlabel('Firing rates(Hz)') - plt.ylabel('Average df/f signal') - plt.title('Signal vs activity plot -Log') - plt.show() -# Check 1: -#Signal_vs_activity(2,np.array([1,2,4,10,50,70,100,200,300,400,500,600,700,800,900,1000,1500,2000,2500,3000,5000])) - -# results: looks plausible -- tho at a big scale it has a linear/exponential trend then it levels off - diff --git a/__pycache__/Functions.cpython-311.pyc b/__pycache__/Functions.cpython-311.pyc index 9908a4672d6a124d8fb48191c9e7077dfdab7c49..9b4a1fc977b82a1b5f9d4a8d058f6c55690adf6f 100644 GIT binary patch delta 19 ZcmeyZ{#%`EIWI340}!;n+Q{`$7yv+j21oz^ delta 19 ZcmeyZ{#%`EIWI340}veS*vR!!7yv*01|k3e diff --git a/__pycache__/b_factory.cpython-311.pyc b/__pycache__/b_factory.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..5db60eca6fe78aab9812b356fff7633744666a2d GIT binary patch literal 4039 zcmeHK-D}%c6u+`9OR{Vyc9MpsUE8FkO-9=!OB(1X8*v?%5EeGdHfBhTY{{+~OZFtW zYk~}U*uX6rxDc2EVK^^cLSNFi^bZ(&jYb~Ag+7fvM|zIg3guY6kpt@@NpjW?zms%MEYAQ9suqI9t6$+ZvoB%ZxyMVsL!`eCY|ha;TBz1 z3R%IB5~@_taw=}!6DWyPPEBQM>-+QL?NS>;DR1r3p``DId*@e*EZy|nX()?cgJ^c! zo%$wWbEa2x&BAe*-bY~ZxM5XfnmihEli`s`&GHC6f>|EJU28?wVC#Jt4Vms)8aQiu z9PT1?77Re+A)Jgrc#&(G2drw%j>C8TcTR|IW4Nw{JCm_4Hw$WOV0j|Iy#%Onz^4)9-lr$@!Q* zlfMQFnU_r8O|3*nGCTw--nx&&rXq{xEj`*8Hux z*mE!9!sVM*$HT8G$5xo-Q){%{eZ}s$vhH5bJPp|GF`JEnFO1!>&}$7Yjo9Iy&2Z1$ z1c-ud<-vuK`&08%Rn$WFI;2EZ{m36EQx8T;BU`PZa`%F#6e&fv0>Sd}2UDe~XQ8%* zE^E|ET3=hqr7=6yvl;3k7QgQ*gbKGfz0z$dn

Yp-Y>gOQlGaMJG;Lqf6&3ZSg`S zvUsWl8@wMZ1*`4I-?`yAwgX4-?nVK&5eGbGBLX<_;3x_-M^FFjMgfgPCXFK`4CGa` zX=gy(K+6Z(Ie7e(SiTur*iiYrb=Hy=`zo0gX<4zmhV1b0T5A2q6U7ckZ6^9UXbBRu z0C`Jj@g9bja1AY54;;k1lrznythtMfTIPHSF-JLbf6jeNXc^MZ!w2%<4bcJ#i4A@) z+f1i5WEB>#R=!y&EEnw)!)yBbmrsgz*O<+W9iiplSD4LcIViX6M+!mIQHA-I&@!y` z0#bhdPtY=BU9ZrK<15S>y&AAj*2}OHF}o#ZGdKP(Xz8Og@(0rw!n9gY@-rG#@d$hK zf+`9+cI8DJQUoKP%^6ua2~{oWfx5=s`K*CyJtaf6>(cX>*AwZ2nlj{^s%zv2O!E*D zB!tw!PHJ}|mh&Q3k(1t>{QrWzrsTvzR{BUA0!#A8p#K7->ZT}Cs8VekC|tX?5cGGi zEp&RL`D)hgx_~3ZIO@XO4dN@99Q+ePjyyNGr33g5@9nj4ieryNZ4?6Hu^Vp CB3Q!! literal 0 HcmV?d00001 diff --git a/__pycache__/b_functions.cpython-311.pyc b/__pycache__/b_functions.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..277b37ae3c82dfc57fb1cfe8961c38e30626fd18 GIT binary patch literal 7584 zcmeHMUu+Y}8K3p8?e*Wp4oQd;zzOKUZ6FD7s8o04FQzFF>2!B+XO*+o-pyu%y>@5U zg`^wFSLG^GrK*vT&Z&3^eMl3*L+&kosVen_M{nhabS0#dPF0(?oX|-p)Q5hvYdaW@ z^eD7SIYzrP-_G~VH}lP(-#0V+d!NsZAWU~07E%Wh`htAaf~!hA`~`#1r-(;9l|nJf zOl>jROqsMj=AbNIEzcQafoDs(((af$?TLASPshA`J8z$-Vm=!ZsA+$#5)tu^-%-mH zLU-Y{LaYg7+`N-#@41$3z=5~QX$Fo*Xa@am-otzE`Iha#ht~=}p%ID&V!>Do$nAl* zRcLC|-S{`(|RuZY*XJyArw6YbfBqP52tPH(0 zE8Ak7U?Eupoz}d6dOiT3%$=S1q}O-illJb!-)T)z@1n@wO|r@+A~vH*kSOVwyU^S)wVe}~Otwwp0#_-+ zp5hWh>Jh6-Y2ip(i9B+T2`QBeNBg4vv&W;UOb!c5Qjn8ESdvviRwQ*MtVp8FrNWt? z2{=p&hL0pt0+*Z$4@6e|2BS!`f?<=h2Ah&(CCeoRgBGr54Oc>ptCFJR1j9Zn;EZC} zF(-?H;p9LQ2Lzo>sRolxWmLm)-E5iRn6Xm3Dyb>KaHIrLkoi@gL7$J08MMmf40|-5 z=B}HGBpXiIT9;5AR4>F+f^2wXs5mVss*qI-3Y1&zHK<9$!xMYslWKgj-=HMJHrcO` z9t=Mo9-ZxbUlFj@ng#dk}sH$PlzTXEjmzh>9F-_%>*T&LH? zKlt=LQJsxKSqu$f_9;EsQwjFu&qJ2Kr8scw)Qu|(SGJJXe#DFl)Rq@{eTA7{ zy?#f4K2*YshnGec_ZPsJ8~%bH46w8FKPtRa`1#W4O81>}Wu{kWdMivX5DysV{ErK> zIuoiep)wP)zBQfcsxVz;ri*-!m;EC;Gg@Ir%gm@L@bQKD3;7G5d)mrO+e7F|nA;{6 zmxIY+lgQ5>8~!_Cpy{uN%p6qsbCnqSO_e|&)`&wxQZ_&#>?EpexCau$Gbv#h(l`bW zAf7P87lYz&7fJZ*p!m<#FaL5eksXXM23wtW{RVwCo*)r6>>`&=bJ9Nmt2d8=x=jKJ z5eM#q1n(o!ei9uZQJ6&C5batI!q$BVqt(5K=}TUFkmvUxwCO>ZZ|OmH_6*w|WM|K? z?LnL`;6Yf+gP7}NllT%$_8FMWZA@D7{BGlT1e2!WPKcfcUmSt#I&>b#r{MDe*&@MNBbY_XkQULBFCZ^#-kI@fk!9sX{hp(7sMmjH?jU-OS-eYbfT17 zJiT;nC3j~|?>xDttY7|YP7j^cooBy8c(i*#=bJV`PXZ1Brq{>lZ*h!%ML-8GtzR5; zPCf_doW!Ja|M*fsXQni|WM3R#ajw}{eR@~z@YosEo1(h&+;<3ca4&Rdgu>*?E7E1y zL54mNx=vP9;|qDJcPGdh4>nN^vWZwJS*p&iXbbYj~9H^FKCS zlQ~t$&60ToVAwf~xf#RHOOulV7GyOJ7aEG;s~$V8J4f?uXt>~jnJ{l0O%j_VhAT&8 z&~V3r4-<+}rZU&e9kS_`%xPk{lT#T*fGZF#w~EZ|p1Gxd#J-iz@VS(522Vg)av`LA z360#MDT-VKP%UK?w4MzFboJRl2g;34Bez-v3W8Ikepp7QpMJJjhQ z%ds*ySq_VOSgb&d_T`<0j^a3!bT<`R3rayLp3$3nDos7Q>tMxoa0^k?6lKNv;SJi8 z|3UFkX<}ud(sgtV|LOSsYvsgbnHF_gtk5D5@N6J9PtUU-yXW0`_m=w&VnSo|-v8Q! z4@?t2pq^m@eOae3SLn-S^Le$}?``aQb0q$s`rnKE!Myu(m{V-wz2fwWZLMcrETgE7 Wq7@Vcf(+M7XVvF#9_FwK3I75_%UyE- literal 0 HcmV?d00001 diff --git a/__pycache__/s_functions.cpython-311.pyc b/__pycache__/s_functions.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..749b72886b956db8c454d38f71d7c06597fe9d2b GIT binary patch literal 3452 zcmaJ@Uu@gP89z!#iIQx|mJ~a#lTNmaTJ+-B2G*cyi*?4F#eZgFtFt|L7zA3R9oizP zkkk?@6sBPwpacdc11X{ceNbLf7wLdK3|Jrc*hkSI5Q9L00(1F*NL#G1 zTOr09>*~T=oEQ;c+6~W{~ZcNSj zX5HA148yxQ9b9okYDw9abeDPH4tDgOb1tXmJiiQbbQxN{yEXS4L(+}8ZCw`B>$|ck z+?D5+BxB#yigSxA8%uK=>zNfr*o6&1ob8H6Zl|Q?O<7Zo>%{{%VRVin=WV6~un~ET zKSNbTuR|yjJBs`@{y1KLXme8zH&tD1htQ>w2HV)TFY~8qzPf&vo_Lwgo~E;ZOWEmJ zCp}wTe|;GZC6DfZtUoC?XYJ8>XLR09EI5gU>N415qpj%D%jnu^bnUC<)^CMYUbUl| z6V+^7ba1hSi?GMQSPNf%1#Y0ck0j_2nSAS47GMof--|e3w$PrB`+huBSdpoO0PDdr z7_kDh#^o47m9R+|)hR!0VF~?yv%*=N7}%$K=pl;%9C9DPN;OU-ktsu{80(rs-BHr^ zXz45~O!tlej(f&LreiyCVNHno+u7~t)`igs*|)pf(VgCRs6_it3U{3o16I`Y*NXP@ z2;MPZu_Pi?Vpgyt4S@BSHLw$;oYKjM(8($AtWG|R%M#ohI*F7n3;fuXT&^IiMp4L1 z{GglgV3rdAf`5?B0|2*@%WGV8m=xAP!85MAH7pgJ-S(!?C5y)=YZyTs{WyqnzjeC&S18U66`#! zY7CuUXCB|E-8dVYsNJqFK3c7zF6OHP8=)EHJh7vxOt#2hlbPj+sCo$3e5Gz3>rw)R~iSc zShj_;`Up8JX4~y2pPzE>Q1SslD2M?-Fdag+Pr%+31oQwPT>w{z1a2$ft^7;ioEAXC z9=h-OVciN@VSpk2BJdo*`@#n>gZk$CdNJ5f{0sB^J;V!X=)ycdU5RY7Cp^IGc8Z0` zd-T8E|J~V7WY59Ei>?_&ZU#WA&;}(HA8^^cR?vVA1z}q%lp`c>WU`9Emxmy<6V6Dg z38@1@dHF5SK#uBDW_ENh^+f!DI8Y!)P) z@CWkQUF|;4LE2MX)-VN~q&P@1ZiqNn(%lF#;jN<{(uAo`k+o5vr6GB_hP-o^j2;gM z$U?$m;u0pS)fF zr87KHT_u@?J6x{)x>jzEoE)_9O&i~I@J$fD#^J-&+Cv+UJ9xZ>$7#Qnn6UArgC|>f zk|ycVF4{(;z3?ZgMSFLEDh`syh-b&BZ7D$i2)MUAsh$yV9YC( zi0Lx}xaKE2NNf>rz7^@mBqV?y{C=Mu76<6V8#J0z2$rwrHKkaP$USK%XUMy%P{?$m zCxDK}G|Hc2N*@5%_?R19%PqO=FLH|TfHtzKXJEv|pfA!0O+I;U9>kI&0fGVKF;P-A zMOM2V67P4++|E3W@Dxu5oP>`dz>)`}D}FfNA$XJwU`1Z)c&jhLlE?5EHK-`81`>6) zE*;10*p=!MiG)cS1(PF<@bR0?sgv-tS$k-{x^k8pZSZz#s+m7od8XOP8`U)cW%_F4 zcXld!vS=q~sUr2plWRX1XkRDkMTfWXjDu%dc;=hfNDGfp-0)0?=6V)f(|M>MJ+g2j z#RD&shKg{z6bv^|6y0ICJyN2{+i+|1SnIX=EUW}-=4H6e1&kuSz|oXOIYm$4b40hK z2lgHJfuf0}g7i~;9afTfZF~Y%JIFB18A`Ug|1*?ob^mARN3H&U|LFOuc4mwjYoq?^ dXUw|{-$s6AF>xl*M*Y?GtITK{eZQiv{|`L)GNu3k literal 0 HcmV?d00001 diff --git a/b_factory.py b/b_factory.py new file mode 100644 index 0000000..5c67dcc --- /dev/null +++ b/b_factory.py @@ -0,0 +1,251 @@ +# This is an identical script to the bleachign analysis +# The only difference is how we calculate the df/f values +# Here we use a moving baseline: F0 is the average of x previous f values + +# import necessary pacakges +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import pandas as pd +import scipy as sp + +# import the necessary functions +from s_functions import simulate_neuron, simulate_nm_conc + + +# define the background fluorescence due to the tissue +bg_tissue = 1.5 + +# ANALYSIS 1: bleaching factor acting on the [NM] contribution to F +def bleach_1(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f = bg_tissue + (K_D*F_min + bleach*nm_conc*F_max)/(K_D + nm_conc) + + # define the signal array + delta_ft_f0 = np.zeros(f.size) + + # calculate f0 values and populate the signal array + for i in range(f.size): + + # calculate f0 by averaging the previous x number of f values + # if x is bigger than the current index then use all the prev f values + # where x is the length of the moving baseline (first element if f[0]) + + if i == 0: + f0 =f[0] + elif i < bline_len: + f0 = np.average(f[:i]) + else: + f0 = np.average(f[i-bline_len:i]) + + + # calculate normalized signal using the calculated f0 + delta_ft_f0[i] = (f[i]-f0)/(f0) + + + # plot the normalized signal delta f/f0 at the different t + # plt.plot(t,delta_ft_f0, label = tau + 1) + # plt.xlabel('time(ms)') + # plt.ylabel('Delta F/F0') + # plt.title('Flourescence intensity signal over time (bleach 1)') + # plt.legend() + + return delta_ft_f0 + +# TESTING IF THE FUNCTION WORKS +# # check this bleaching effect 1 for different values of tau -- different bleach factor +# different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# # generate an input nm conc array from firing neuron +# firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +# nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# # plot bleached signal with different time constants for the bleach factor +# for i in range(len(different_taus)): +# bleach_1(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc, bline_len=5000) +# plt.show() + + + +# ANALYSIS 2: bleaching factor acting on the contributions of the dye, F0 and [NM] +def bleach_2(K_D, tau, F_max, F_min, nm_conc, bline_len =5000): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f= bg_tissue+ bleach*(K_D*F_min + nm_conc*F_max)/(K_D + nm_conc) + + # define the signal array + delta_ft_f0 = np.zeros(f.size) + + # calculate f0 values and populate the signal array + for i in range(f.size): + + # calculate f0 by averaging the previous x number of f values + # if x is bigger than the current index then use all the prev f values + # where x is the length of the moving baseline (first element if f[0]) + + if i == 0: + f0 =f[0] + elif i < bline_len: + f0 = np.average(f[:i]) + else: + f0 = np.average(f[i-bline_len:i]) + + # calculate normalized signal using the calculated f0 + delta_ft_f0[i] = (f[i]-f0)/(f0) + + # plot the normalized signal delta f/f0 at the different t + # plt.plot(t,delta_ft_f0, label = tau + 2) + # plt.xlabel('time(ms)') + # plt.ylabel('Delta F/F0') + # plt.title('Flourescence intensity signal over time (bleach 2)') + # plt.legend() + + return delta_ft_f0 + +# TESTING IF THE FUNCTION WORKS +# # check this bleaching effect 1 for different values of tau -- different bleach factor +# different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# # generate an input nm conc array from firing neuron +# firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +# nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# # plot bleached signal with different time constants for the bleach factor +# for i in range(len(different_taus)): +# bleach_2(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +# plt.show() + + + +# ANALYSIS 3: bleaching factor acting on the background fluorescence +def bleach_3(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f= bleach*bg_tissue + (K_D*F_min + nm_conc*F_max)/(K_D + nm_conc) + + # define the signal array + delta_ft_f0 = np.zeros(f.size) + + # calculate f0 values and populate the signal array + for i in range(f.size): + + # calculate f0 by averaging the previous x number of f values + # if x is bigger than the current index then use all the prev f values + # where x is the length of the moving baseline (first element if f[0]) + + if i == 0: + f0 =f[0] + elif i < bline_len: + f0 = np.average(f[:i]) + else: + f0 = np.average(f[i-bline_len:i]) + + # calculate normalized signal using the calculated f0 + delta_ft_f0[i] = (f[i]-f0)/(f0) + + # plot the normalized signal delta f/f0 at the different t + # plt.plot(t,delta_ft_f0, label = tau + 3) + # plt.xlabel('time(ms)') + # plt.ylabel('Delta F/F0') + # plt.title('Flourescence intensity signal over time (bleach 3)') + # plt.legend() + + return delta_ft_f0 + +# TESTING IF THE FUNCTION WORKS +# # check this bleaching effect 1 for different values of tau -- different bleach factor +# different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# # generate an input nm conc array from firing neuron +# firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +# nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# # plot bleached signal with different time constants for the bleach factor +# for i in range(len(different_taus)): +# bleach_3(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +# plt.show() + + + + +# ANALYSIS 4: bleaching on all contributions +def bleach_4(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f= bleach*(bg_tissue + (K_D*F_min + nm_conc*F_max)/(K_D + nm_conc)) + + # define the signal array + delta_ft_f0 = np.zeros(f.size) + + # calculate f0 values and populate the signal array + for i in range(f.size): + + # calculate f0 by averaging the previous x number of f values + # if x is bigger than the current index then use all the prev f values + # where x is the length of the moving baseline (first element if f[0]) + + if i == 0: + f0 =f[0] + elif i < bline_len: + f0 = np.average(f[:i]) + else: + f0 = np.average(f[i-bline_len:i]) + + # calculate normalized signal using the calculated f0 + delta_ft_f0[i] = (f[i]-f0)/(f0) + + # plot the normalized signal delta f/f0 at the different t + # plt.plot(t,delta_ft_f0, label = tau + 4) + # plt.xlabel('time(ms)') + # plt.ylabel('Delta F/F0') + # plt.title('Flourescence intensity signal over time (bleach 4)') + # plt.legend() + + return delta_ft_f0 + +# TESTING IF THE FUNCTION WORKS +# # check this bleaching effect 1 for different values of tau -- different bleach factor +# different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# # generate an input nm conc array from firing neuron +# firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +# nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# # plot bleached signal with different time constants for the bleach factor +# for i in range(len(different_taus)): +# bleach_4(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +# plt.show() + diff --git a/b_functions.py b/b_functions.py new file mode 100644 index 0000000..fcef438 --- /dev/null +++ b/b_functions.py @@ -0,0 +1,253 @@ +# This is an identical script to the bleachign analysis +# The only difference is how we calculate the df/f values +# Here we use a moving baseline: F0 is the average of x previous f values + +# import necessary pacakges +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import pandas as pd +import scipy as sp + +# import the necessary functions +from s_functions import simulate_neuron, simulate_nm_conc + + +# define the background fluorescence due to the tissue +bg_tissue = 1.5 + +# ANALYSIS 1: bleaching factor acting on the [NM] contribution to F +def bleach_1(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f = bg_tissue + (K_D*F_min + bleach*nm_conc*F_max)/(K_D + nm_conc) + + # define the signal array + delta_ft_f0 = np.zeros(f.size) + + # calculate f0 values and populate the signal array + for i in range(f.size): + + # calculate f0 by averaging the previous x number of f values + # if x is bigger than the current index then use all the prev f values + # where x is the length of the moving baseline + + if i < bline_len: + f0 = np.average(f[:i]) + else: + f0 = np.average(f[i-bline_len:i]) + + # calculate normalized signal using the calculated f0 + delta_ft_f0[i] = (f[i]-f0)/(f0) + + + # plot the normalized signal delta f/f0 at the different t + plt.plot(t,delta_ft_f0, label = tau + 1) + plt.xlabel('time(ms)') + plt.ylabel('Delta F/F0') + plt.title('Flourescence intensity signal over time (bleach 1)') + plt.legend() + + return delta_ft_f0 + + +# check this bleaching effect 1 for different values of tau -- different bleach factor +different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# generate an input nm conc array from firing neuron +firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# plot bleached signal with different time constants for the bleach factor +for i in range(len(different_taus)): + bleach_1(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc, bline_len=5000) +plt.show() + + + +# ANALYSIS 2: bleaching factor acting on the contributions of the dye, F0 and [NM] +def bleach_2(K_D, tau, F_max, F_min, nm_conc, bline_len =5000): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f= bg_tissue+ bleach*(K_D*F_min + nm_conc*F_max)/(K_D + nm_conc) + + # define the signal array + delta_ft_f0 = np.zeros(f.size) + + # calculate f0 values and populate the signal array + for i in range(f.size): + + # calculate f0 by averaging the previous x number of f values + # if x is bigger than the current index then use all the prev f values + # where x is the length of the moving baseline + + if i < bline_len: + f0 = np.average(f[:i]) + else: + f0 = np.average(f[i-bline_len:i]) + + # calculate normalized signal using the calculated f0 + delta_ft_f0[i] = (f[i]-f0)/(f0) + + # plot the normalized signal delta f/f0 at the different t + plt.plot(t,delta_ft_f0, label = tau + 2) + plt.xlabel('time(ms)') + plt.ylabel('Delta F/F0') + plt.title('Flourescence intensity signal over time (bleach 2)') + plt.legend() + + return delta_ft_f0 + + +# check this bleaching effect 1 for different values of tau -- different bleach factor +different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# generate an input nm conc array from firing neuron +firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# plot bleached signal with different time constants for the bleach factor +for i in range(len(different_taus)): + bleach_2(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +plt.show() + + + +# ANALYSIS 3: bleaching factor acting on the background fluorescence +def bleach_3(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f= bleach*bg_tissue + (K_D*F_min + nm_conc*F_max)/(K_D + nm_conc) + + # define the signal array + delta_ft_f0 = np.zeros(f.size) + + # calculate f0 values and populate the signal array + for i in range(f.size): + + # calculate f0 by averaging the previous x number of f values + # if x is bigger than the current index then use all the prev f values + # where x is the length of the moving baseline + + if i < bline_len: + f0 = np.average(f[:i]) + else: + f0 = np.average(f[i-bline_len:i]) + + # calculate normalized signal using the calculated f0 + delta_ft_f0[i] = (f[i]-f0)/(f0) + + # plot the normalized signal delta f/f0 at the different t + plt.plot(t,delta_ft_f0, label = tau + 3) + plt.xlabel('time(ms)') + plt.ylabel('Delta F/F0') + plt.title('Flourescence intensity signal over time (bleach 3)') + plt.legend() + + return delta_ft_f0 + + +# check this bleaching effect 1 for different values of tau -- different bleach factor +different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# generate an input nm conc array from firing neuron +firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# plot bleached signal with different time constants for the bleach factor +for i in range(len(different_taus)): + bleach_3(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +plt.show() + + + + +# ANALYSIS 4: bleaching on all contributions +def bleach_4(K_D, tau, F_max, F_min, nm_conc, bline_len=5000): + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # bleaching factor -- starts off as 1 then exponentially decreases + # we set tau to be a very large constant so this is a slow decrease + bleach = np.exp(-t/tau) + + # calculate bleached f values: derived in part from eq 2 in Neher/Augustine + f= bleach*(bg_tissue + (K_D*F_min + nm_conc*F_max)/(K_D + nm_conc)) + + # define the signal array + delta_ft_f0 = np.zeros(f.size) + + # calculate f0 values and populate the signal array + for i in range(f.size): + + # calculate f0 by averaging the previous x number of f values + # if x is bigger than the current index then use all the prev f values + # where x is the length of the moving baseline + + if i < bline_len: + f0 = np.average(f[:i]) + else: + f0 = np.average(f[i-bline_len:i]) + + # calculate normalized signal using the calculated f0 + delta_ft_f0[i] = (f[i]-f0)/(f0) + + # plot the normalized signal delta f/f0 at the different t + plt.plot(t,delta_ft_f0, label = tau + 4) + plt.xlabel('time(ms)') + plt.ylabel('Delta F/F0') + plt.title('Flourescence intensity signal over time (bleach 4)') + plt.legend() + + return delta_ft_f0 + + +# check this bleaching effect 1 for different values of tau -- different bleach factor +different_taus = np.array([10e6,10e5,10e4,10e3,10e2,10e1,10e0]) + +# generate an input nm conc array from firing neuron +firing_neuron = simulate_neuron(n_timesteps=70000,firing_rate=1) +nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# plot bleached signal with different time constants for the bleach factor +for i in range(len(different_taus)): + bleach_4(K_D = 1000, tau=different_taus[i], F_max = 45, F_min = 10, nm_conc=nm_conc) +plt.show() + + +# ANALYSIS 5: comparing the different contributions at one timescale for the bleach factor + +# choose the timescale for the bleach factor +chosen_tau = 10e4 + +bleach_1(K_D=1000,tau=chosen_tau, F_max=45, F_min=10, nm_conc=nm_conc) +bleach_2(K_D=1000,tau=chosen_tau, F_max=45, F_min=10, nm_conc=nm_conc) +bleach_3(K_D=1000,tau=chosen_tau, F_max=45, F_min=10, nm_conc=nm_conc) +bleach_4(K_D=1000,tau=chosen_tau, F_max=45, F_min=10, nm_conc=nm_conc) +plt.show() \ No newline at end of file diff --git a/f_rate.py b/f_rate.py new file mode 100644 index 0000000..dbb9eb2 --- /dev/null +++ b/f_rate.py @@ -0,0 +1,67 @@ +# This script is used for analysis of the effect of +# changing the firing rate of the neuron on the measured signal + +# import necessary pacakges +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import pandas as pd +import scipy as sp + +from Functions import simulate_neuron, simulate_nm_conc, simulate_flourescence_signal + + +# ANALYSIS 1: + +# Define function that makes generates several firing neurons with different firing rates +# and then compares the average df/f signal across them + +def signal_vs_activity(number_neurons,firing_rates): + + # create array to store the average signals generated by the neurons + average_signals = [] + + # for each firing rate, simulate a neuron, the nm dynamics from it and the signal + # then get the average signal + for i in range(firing_rates.size): + + # simulate neuron with different firing rate + guinea_neuron = simulate_neuron(70000,firing_rates[i]) + + # generate the nm_conc + guinea_nm_conc, guinea_b_conc, guinea_c_conc = simulate_nm_conc(guinea_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + + # then generate the signal + guinea_signal = simulate_flourescence_signal(K_D = 1000, F_max = 45, F_min = 10, nm_conc=guinea_nm_conc) + + # get the average of this signal and add it the original array + average_signals.append(np.average(guinea_signal)) + + print('simulated {} neuron(s)'.format(i+1)) + + + # plot the different firing rates vs their signals + plt.scatter(firing_rates,average_signals) + plt.xlabel('Firing rates(Hz)') + plt.ylabel('Average df/f signal') + plt.title('Signal vs activity plot') + plt.show() + + # plot at log scale + plt.loglog(firing_rates,average_signals, 'o') + plt.xlabel('Firing rates(Hz)') + plt.ylabel('Average df/f signal') + plt.title('Signal vs activity plot - Log') + plt.show() + + +# Check 1: +signal_vs_activity(2,np.array([1,2,4,10,50,70,100,200,300,400,500,600,700,800,900,1000,1500,2000,2500,3000,5000])) + +# results: looks plausible -- tho at a big scale it has a linear/exponential trend then it levels off +# -- it's linear + +# TO ADD: +# constrict to the region , 1000 hz +# then also get an r2 estimate + diff --git a/s_functions.py b/s_functions.py new file mode 100644 index 0000000..d093251 --- /dev/null +++ b/s_functions.py @@ -0,0 +1,180 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import pandas as pd +import scipy as sp + +# This script contains the functions in the Simulation notebook +# Its purpose is to run all the functions and produce a df/f plot + +# TO DO: make separate functions to plot the results of each of these functions +# print output line after function 2 to summarize the result: average conc for each of them maybe +# make the function 1 able to make multiple neurons at once! +# add the firing rate to the plots of flourescence + +# check the inputs to each of the functions and raise errors if something's wrong + +# Bonus: find way to show the figures all at once + + +# Function 1: firing neuron +def simulate_neuron(n_timesteps, firing_rate, number=1): + + # generate a base array with random numbers between 0-1 + x = np.random.rand(n_timesteps) + + # Then populate the bins with signals - firing & not firing -- with a specific probability that you choose + firing_neuron = x < 0.001*firing_rate + firing_neuron = firing_neuron.astype(int) + + + # # then make a plot of it! + # plt.plot(firing_neuron) + # plt.xlabel('timesteps') + # plt.ylabel('Neuron activity') + # plt.title('Neuron Activity over {} timesteps'.format(n_timesteps)) + # plt.show() + + + # check exactly how many spikes were produced: to see if it works + n_spikes = np.size(np.nonzero(firing_neuron)) + + # print simulated neuron summary: + print('Simulated neuron with {} spikes in {} timesteps ({} Hz).'.format(n_spikes, n_timesteps, firing_rate)) + + + return firing_neuron + +# output 1 +#firing_neuron = simulate_neuron(70000,1) + + +# Function 2: takes in an array of neuron activity and gives corresponding [NM] +def simulate_nm_conc(neuron_activity,nm_conc0, k_b,k_r,gamma): + + + # create array of [NM] with same size as neuron activity + nm_conc = np.zeros(neuron_activity.size) + + # define delta_nm, the increase in [NM] due to a spike + # this will be a constant value -- calculate the amplitude of the exponential function, A + delta_nm = 1 + + # first define tau the time constant + tau = (1+k_r+k_b)/gamma + + # create a for-loop where we update the value of [NM] for current timestep + for t in range(neuron_activity.size): + + # first timebin condition + if t == 0 : + nm_conc[t] = nm_conc0 + else: + nm_conc[t] = nm_conc[t-1] + + # update [NM] value + + # define d_nm/dt, the inifinitesimal decay in [NM] in one timestep + d_nm_dt = (nm_conc[t]-nm_conc0)/tau + + # if there's a spike add Delta_nm else subtract d_nm/dt + if neuron_activity[t]==1: nm_conc[t] = nm_conc[t] + delta_nm + else: nm_conc[t] = nm_conc[t] - d_nm_dt + + # plot the [NM] at all timesteps + n_timesteps = neuron_activity.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # Calculate the concentrations of the bound forms of the NM + # start with [NM B] the NM bound to the sensor + nm_b_conc = k_b*nm_conc + + # then [NM R], the NM bound to the receptor + nm_r_conc = k_r*nm_conc + + # # plot [NM], [NM B] and [NM R] simulataneously + # plt.plot(t,nm_conc, color = 'b', label='[NM]') + # plt.plot(t,nm_b_conc, color = 'g', label='[NM B]') + # plt.plot(t,nm_r_conc, color = 'r', label='[NM R]') + + # # label the axes and make legend + # plt.xlabel('time (ms)') + # plt.ylabel('Concentration') + # plt.title('NM concentration across {} ms'.format(n_timesteps)) + # plt.legend() + # plt.show() + + + # return the array of the [NM], [NM B], and [NM R] + return nm_conc, nm_b_conc, nm_r_conc + +# output 2 +#nm_conc, nm_b_conc, nm_r_conc = simulate_nm_conc(firing_neuron,nm_conc0=0,k_b=0.6, k_r=0.4,gamma=0.004) + +# Function 2.1: produces zoomed in plot +def plot_nm_conc(nm,start,stop,colour='b', plotlabel = ''): + + # define the timesteps to plot the [NM] + timesteps = stop - start + 1 + t = np.linspace(start,stop,timesteps) + + # get that section of the [NM] array + nm_section = nm[start:stop+1] + + # plot the [NM] + plt.plot(t,nm_section, color=colour, label=plotlabel) + plt.xlabel('time (ms)') + plt.ylabel('NM concentration') + plt.title('NM {} concentration from {} to {} ms'.format(plotlabel, start,stop)) + plt.show() + +# output 2.1 +# plot_nm_conc(nm_conc, start = 22000,stop = 26000) +# plot_nm_conc(nm_b_conc, start = 22000,stop = 26000, colour='g', plotlabel='B') +# plot_nm_conc(nm_r_conc, start = 22000,stop = 26000, colour='r', plotlabel='R') + + +# function 3: get that signal! +def simulate_flourescence_signal(K_D, F_max, F_min, nm_conc): + + # define K_D prime as + K_Dp = K_D*(F_min/F_max) + + # the initial/steady state concentration, [NM]i,0, of the neuromdultor + # CONFIRM VALUE FROM KENTA + nm_conc_0 = 0 + + # define the numerator and denominator + numerator = (K_Dp + nm_conc)/(K_D + nm_conc) + denominator = (K_Dp + nm_conc_0)/(K_D + nm_conc_0) + + # derive delta f/f0 by plugging in + delta_ft_f0 = (numerator/denominator) - 1 + + # create timesteps array for the plot + n_timesteps = nm_conc.size + t = np.linspace(0,n_timesteps-1,n_timesteps) + + # plot the normalized signal delta f/f0 at the different t + # plt.plot(t,delta_ft_f0) + # plt.xlabel('time(ms)') + # plt.ylabel('Delta F/F0') + # plt.title('Flourescence intensity signal over time') + # plt.show() + + print('completed f_signal simulation') + + return delta_ft_f0 + + +# output 3 +#f_signal = simulate_flourescence_signal(K_D = 1000, F_max = 45, F_min = 10, nm_conc=nm_conc) + + + + + + + + +