import csv import pandas as pd from scipy.signal import savgol_filter from sage.plot.matrix_plot import MatrixPlot from sympy import Symbol, cos, sympify, pprint from sympy.abc import x import numpy as np from sage.symbolic.integration.integral import indefinite_integral #import matplotlib.pyplot as plt import re #from ring import SR #from expression import Expression #Real data from SK Voltage versus time # -50mV with open('trace1.txt', 'r') as f: trace1 = [RDF(l.strip()) for l in f.readlines()] #-40mV with open('trace2.txt', 'r') as f: trace2 = [RDF(l.strip()) for l in f.readlines()] #-30mV with open('trace3.txt', 'r') as f: trace3 = [RDF(l.strip()) for l in f.readlines()] #-20mV with open('trace4.txt', 'r') as f: trace4 = [RDF(l.strip()) for l in f.readlines()] with open('time.txt', 'r') as f: time = [RDF(l.strip()) for l in f.readlines()] #https://en.wikipedia.org/wiki/Taylor_series The exponential function of the SK current is written as a Taylor series expansion #Smoothing out the exponential function as if it were a polynomial using Savitzky- Golay filter #13868 total points with wanting to display 500: 13868/500 = 27.7 trace1_smooth = np.array(trace1) #Window size must be larger than the order of the polynomial window_size = 25 order = 1 if window_size % 2 == 0: window_size = window_size -1 #6842 points trace1_smooth = savitzky_golay(trace1_smooth,window_size,order) trace1_smooth = numpy.round(trace1_smooth, decimals=3) sp1_smooth = scatter_plot(zip(time, trace1_smooth), markersize=1, marker='o') trace2_smooth = np.array(trace2) trace2_smooth = savitzky_golay(trace2_smooth,window_size,order) trace2_smooth = numpy.round(trace2_smooth, decimals=3) sp2_smooth = scatter_plot(zip(time, trace2_smooth), markersize=1, marker='o') trace3_smooth = np.array(trace3) trace3_smooth = savitzky_golay(trace3_smooth,window_size,order) trace3_smooth = numpy.round(trace3_smooth, decimals=3) sp3_smooth = scatter_plot(zip(time, trace3_smooth), markersize=1, marker='o') trace4_smooth = np.array(trace4) trace4_smooth = savitzky_golay(trace4_smooth,window_size,order) trace4_smooth = numpy.round(trace4_smooth, decimals=3) sp4_smooth = scatter_plot(zip(time, trace4_smooth), markersize=1, marker='o') all_smoothed_plots = sp4_smooth + sp3_smooth + sp2_smooth + sp1_smooth all_smoothed_plots.show(title="Smoothed Data Plots", gridlines="minor") #END OF SMOOTHING PROCESS AND DISPLAY THE 4 SMOOTHED PLOTS #Now lets try to find an equation #http://tutorial.math.lamar.edu/Classes/DE/RepeatedRoots.aspx o = open('trace1_truncated_smoothed.txt','w') for x in trace1_smooth: o.write(str(x)+'\n') o.close() o = open('trace1_straightup.txt','w') for x in trace1: #o.write(str(5) o.write(str(x)+'\n') o.close() o = open('time_truncated.txt','w') for x in time: #o.write(str(5) o.write(str(x)+'\n') o.close() #trace1_smooth #np.array(trace1) sp1 = scatter_plot(zip(time, trace1), markersize=1, marker='o') #sp1.show() sp2 = scatter_plot(zip(time, trace2), markersize=1, marker='s') #sp2.show() sp3 = scatter_plot(zip(time, trace3), markersize=1, marker='x') #sp3.show() sp4 = scatter_plot(zip(time, trace4), markersize=1, marker='+') #sp4.show() sp_total = sp1 + sp2 + sp3 + sp4; #sp_total.title() #sp_total.show(title= "Raw Data Plots") #data1 = list(csv.reader(f) #data1 = [row for row in data1] #M = matrix(data1) #data1 #M #plot(I_SK_K(-90*10^-3,u),(u,0,4*10^-7),title="I SK Plot at -90 milliVolts") #table(data1) #matrix_plot(numpy.random.rand(10, 10)) #matrix_plot(data1, fontsize=10) #matrix_plot(data1).show(fontsize=10) #plot(I_SK_K(90*10^-3,u),(u,0,4*10^-7),title="Tracer1") #s = scatter_plot(data1, marker='s') #f = open('trace2.txt', 'r') #data2 = csv.reader(f) #data2 = [row for row in data2] #data2 #f = open('trace3.txt', 'r') #data3 = csv.reader(f) #data3 = [row for row in data3] #data3 #f = open('trace4.txt', 'r') #data4 = csv.reader(f) #data4 = [row for row in data4] #data4 # Calcium dependent SK potassium current centi = 10^-2 milli = 10^-3 mirco = 10^-6 nano = 10^-9 pico = 10^-12 u = var('u') V = var('V') #100 mV this is from coefficents sheet E_Ca = 100 * milli #500nM #K_pump = 500 * 10^-9 * 10^-6/ 10^-3 #Used in equation 11 K_pump = 500 * nano #500nM*micro m/ms #Used in equation 11 M_pump = 500 * nano * mirco / milli #From coefficents sheet 2 ms/cm^2 #g_SK = 2 * 10^-3/100^2 g_SK = 2 * milli / centi^2 #From coefficent sheet 125.8nM #K_1 = 125.8 * 10^ -9 K_1 = 125.8 * nano #from coefficent sheet E_K = -90.0 * milli #Cystosolic buffering constant F_Ca = .01 # H is Faraday's Number H = .0193 #radius of neuron r = 20 * mirco #Equation 11 I_pump = (M_pump*u)/(u+K_pump) #L type calcium current #a_c = -0.0032*(V+50)/(exp(-(V+50)/5)-1) #b_c = exp(-(V+55)/(40)) #Take from equations sheet #g_L_Ca = .08 * milli / centi^2 #G_L_Ca = g_L_Ca*(a_c/(a_c + b_c))^4 #Equation 8 #I_L_Ca = G_L_Ca*(E_Ca - V) #Equation 10 change_in_calcium_with_respect_to_time=2*F_Ca/r *(I_L_Ca/H - I_pump) calcium_function_done_by_integration=indefinite_integral(change_in_calcium_with_respect_to_time,u) #need to solve this differential equation stringnewline ="\n\n\n" #string1= "This is du/dt:" #print string1 #change_in_calcium_with_respect_to_time #view (change_in_calcium_with_respect_to_time) #print stringnewline #string2= "This is u:" #print string2 #view(calcium_function_done_by_integration) #print stringnewline #string4="Write calcium (u) in terms of Voltage:" #As of now we have integrated du/dt and gotten the expression u which is in terms of Voltage and calcium (u) #Now we are taking the expression and expressing the calciumn in terms of voltage #print string4 #eq1 = calcium_function_done_by_integration == u #S_testing = solve(eq1, u, solution_dict=True) #S = solve(eq1, u) #calcium_in_term_of_voltage = S; calcium_in_term_of_voltage #calcium_in_term_of_voltage_testing = S_testing[0][u]; calcium_in_term_of_voltage_testing print stringnewline #string3="Let voltage be clamped at -20mV and calciumn will be:" #print string3 #calcium_in_term_of_voltage_testing.subs(V == -20*10^-3).n() #calcium_in_term_of_voltage[0][u].n() #calcium_in_term_of_voltage=calcium_in_term_of_voltage._sympy_() #u=x^2+7*x*x.subs_expr(x^2 == w) #u #Removing the u terming and taking the righthandside of the equation to get Voltage to minus 20 #Calcium_at_minus_20=calcium_in_term_of_voltage[0].rhs();Calcium_at_minus_20 #sin(x).subs(x=5).n() this is a symbolic expression evaluated numerically #Calcium_at_minus_20=calcium_in_term_of_voltage[0].rhs().subs(V == -20*10^-3); Calcium_at_minus_20 #Calcium_at_minus_20.expand_sum().numerical_approx(200) #view(Calcium_at_minus_20) #Calcium_at_minus_20.n() #view(calcium_in_term_of_voltage) #Calcium_at_minus_20=solve(calcium_in_term_of_voltage);Calcium_at_minus_20 #eqn = x^3 + 2/3 >= x - pi #eqn.lhs() #calcium_in_terms_of_voltage=calcium_in_term_of_voltage.left_hand_side() #Calcium_at_minus_20 = calcium_in_terms_of_voltage.subs(V == -20*10^-3) #print Calcium_at_minus_20 #calcium_function_done_by_integration_minus20=calcium_function_done_by_integration.subs(V == -20*10^-3) #calcium_function_done_by_integration_minus20 #eq1 = calcium_function_done_by_integration_minus20 == u #S = solve(eq1,u); S #T = solve(x*3 + 6 == x, x); T #calcium_function_done_by_integration_minus20 #eq1= calcium_function_done_by_integration_minus20 #calcium_at_minus_20=solve([u == calcium_function_done_by_integration_minus20],u) #print calcium_at_minus_20 print stringnewline print 1505003284098320 #Under equation 9 G_SK(u) =g_SK*(u^4/(u^4+K_1^4)) view(G_SK) I_SK_K(V, u) = G_SK*(E_K - V) view(I_SK_K) #This is a plot of the calcium #u #plot(u, (u,0,1000),axes_labels=['$Calcium$','$F(Calcium)$'], title="Calcium plot") #This is a plot of the I_pump #I_pump #plot(I_pump,(I_pump,0,1*10^-9),axes_labels=['$I pump$','$F( I Pump)$'], title="I pump plot") #g = var('g') #diff(sin(g), g) #G_SK #plot(G_SK,(u,0,4*10^-7),axes_labels=['$Calcium$','$Conductance$'],title="Small Conductance Potassium Channel's Conductance") ##range of the plot for calcium is from 0 to 100 micro Moles #I_SK_K #scale_V = -3 #scale_u = -7 #plot3d(I_SK_K,(V,89 * 10^scale_V,90*10^scale_V),(u,0,10^scale_u), title="I SK_K Plot") #-(2.00000000000000e-7)*(1.00000000000000*V + 0.0900000000000000)*u^4/(u^4 + 2.50450881409600e-28) #x1 = numpy.array([-20,-30,-40,-50]) #A=numpy.array([[x1, 1] for j in range(2)]) #n = numpy.array(range(2*4), dtype=int) #n.reshape((4,2)) #n.shape #a = numpy.array([[-20,0],[-30,0],[-40,0],[-50,0]]) a = numpy.array([[-20,1e-1],[-30,1e-8],[-40,1e-9],[-50,1e-9]]) for j in range (len(a)): a[j][0] = a[j][0]*milli*1 a #linestyling = numpy.array([['--'],['-.'],['-'],[':']]) #color = numpy.array([green],[purple],[yellow],[orange]) #a #a.shape #print "hell" #x = 0 #Voltage = -20 #print numpy.shape(A), numpy.ndim(A) #for j in range(numpy.ndim(A)): # A[j][0] = Voltage # Voltage = Voltage - 10 # print A[j] #print a[0][0],a[1][0],a[2][0],a[3][0] #I_SK_K(V,u) print stringnewline plot_curves_of_SK = Graphics() for j in range(len(a)): #print j #a[j] = (Voltage, I_SK_K(Voltage*10^-3,u),(u,0,4*10^-7)) Voltage = a[j][0] Calcium = 1*10^-7 #Storing current a[j][1]=I_SK_K(Voltage,Calcium) #title = "I SK Plot at " + str(Voltage) + " millivolts" #print title #Holding calcium constant at 2*10^-7, we will store the SK current #temp_color = color[j] #tempstyling = linestyling[j] #I_SK_K(Voltage,Calcium) plot_curves_of_SK = plot(I_SK_K(Voltage,u),(u,0,1*mirco),axes_labels=['$Calcium$','$Current in Nanoamps$'],title="I SK Plot from -20mV to -50mV") + plot_curves_of_SK print "Indeed, in the presence of 100 nm free Ca2+, the amplitude of the current was −288±72, −312±180, −527±149 pA for SK1, SK2 and SK3, respectively (n=3). \n In the presence of 300 nm Ca2+, the amplitude was −505±267, −483±226, −463±233 pA for SK1, SK2 and SK3, respectively (n=3). \n In the presence of 1 μm free intracellular Ca2+, currents had an average value of −2418±902, −2215±690 and −1783±997 pA, for SK1, SK2 and SK3, respectively." #x = x+1 #a plot_curves_of_SK.show() print stringnewline #Voltage versus SK current #plot(I_SK_K(Voltage*10^-3,2*10^-7),(u,0,4*10^-7),axes_labels=['$Voltage$','$SK Current$'],title="Voltage versus SK") scatter_plot_of_Voltage_versus_SK = Graphics() scatter_plot_of_Voltage_versus_SK=scatter_plot(a, marker='s',axes_labels=['$Voltage Step$','$SK Current$']) #x = numpy.array([float(0),0,0,0]) #y = numpy.array([float(0),0,0,0]) #for j in range (len(a)): # x[j] = a[j][0] # y[j] = a[j][1] #x #y #plt.scatter(x,y, alpha=.075, facecolor='black', marker='s') #plt.axis([numpy.amin(x),numpy.amax(x),numpy.amin(y),numpy.amax(y)]) #plt.show() scatter_plot_of_Voltage_versus_SK.show() print stringnewline #plot(I_SK_K(90*10^-3,u),(u,0,4*10^-7),title="I SK Plot at 90 milliVolts") #plot(I_SK_K(-90*10^-3,u),(u,0,4*10^-7),title="I SK Plot at -90 milliVolts") #plot3d(I_SK_K,(V,-90 * 10^scale_V,90*10^scale_V),(u,0,10^scale_u), title="I SK_K Plot")
1505003284098320
u ↦ u4+2.50450881409600×10−2820u4
(V,u) ↦ −u4+2.50450881409600×10−2820.0000000000000(1.00000000000000V+0.0900000000000000)u4
array([[ -2.00000000e-02, 1.00000000e-01],
[ -3.00000000e-02, 1.00000000e-08],
[ -4.00000000e-02, 1.00000000e-09],
[ -5.00000000e-02, 1.00000000e-09]])
Indeed, in the presence of 100 nm free Ca2+, the amplitude of the current was −288±72, −312±180, −527±149 pA for SK1, SK2 and SK3, respectively (n=3).
In the presence of 300 nm Ca2+, the amplitude was −505±267, −483±226, −463±233 pA for SK1, SK2 and SK3, respectively (n=3).
In the presence of 1 μm free intracellular Ca2+, currents had an average value of −2418±902, −2215±690 and −1783±997 pA, for SK1, SK2 and SK3, respectively.
%html (double click to edit)
(double click to edit)
︠c0069e73-ec4b-4d95-9961-f4c73ad7de43︠