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 9908a46..9b4a1fc 100644 Binary files a/__pycache__/Functions.cpython-311.pyc and b/__pycache__/Functions.cpython-311.pyc differ diff --git a/__pycache__/b_factory.cpython-311.pyc b/__pycache__/b_factory.cpython-311.pyc new file mode 100644 index 0000000..5db60ec Binary files /dev/null and b/__pycache__/b_factory.cpython-311.pyc differ diff --git a/__pycache__/b_functions.cpython-311.pyc b/__pycache__/b_functions.cpython-311.pyc new file mode 100644 index 0000000..277b37a Binary files /dev/null and b/__pycache__/b_functions.cpython-311.pyc differ diff --git a/__pycache__/s_functions.cpython-311.pyc b/__pycache__/s_functions.cpython-311.pyc new file mode 100644 index 0000000..749b728 Binary files /dev/null and b/__pycache__/s_functions.cpython-311.pyc differ 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) + + + + + + + + +