#Fitz-Nagumo Model var("V","w") #declares V and w as variables to be used in equations t=srange(0,200,.1) #creates a time range for which the simulation will occur. soln=desolve_odeint([V*(1-V)*(V-.2)-w+.6, .01*(V-(.5*w))], ics= [.8,.4], times=t, dvars=[V,w]) #creates a simulation for the two functions for voltage and current, producing values to store in soln. initial condiitons of .8 and .4 for time t, with the first equation being solved in respect to V and the second in respect to w. #****NOTE: We changed the value of Iext from the given value '0.1' to '0.6' to create a trajectory image that showed oscillatory behavior. With a value of 0.1, Iext is not high enough to cross the threshold to start the oscillatory behavior.
(V, w)
list_plot(soln, axes_labels=("V","w"), plotjoined=true) #plots the two functions' solutions against each other to create a trajectory
-def eulersthresh(init,init2, stepsize,numtime): #sets up a function that takes in two inital values, the stepsize, and the numtime numsteps= floor(numtime/stepsize) #finds the amount of steps that will be taken for the model, floor allows the whole number to be divided by a decimal listt = [init] #creates a list for values of voltage, and sets the first place as initial value inputed list2 = [init2] #creates a list to store the values of current and enters init2 as the first value functionEM(N,w)=N*(1-N)*(N-0.2)-w +0.16 #establishes the function for voltage with N=voltage and w=current functionE(N,w)= .01*(N-.5*w) #establishes the function for current vals = srange(0, numsteps+1,1) #creates the range of values for which the following for loop will run t=srange(0, numsteps*stepsize+stepsize,stepsize) #creates the time values that will be used for the plot for v in vals: #creates a for loop that will run for the length of vals if listt[v]>=1: #creates an if statement where if the present value in listt is greater or equal to one listt.append(init) #add the inital input value to the list of voltage values to bring the function back down list2.append(init2) #add the inital current value to the values of current ans= (zip(t,listt)) #creates the pairs of time and voltage to be plotted else: #if the if statement is not true g=listt[v]+ (stepsize*functionEM(listt[v],list2[v])) #runs euler's method for the voltage function and sets the next value of the function equal to g listt.append(g) #adds g to the list of voltage values h= list2[v] + (stepsize*functionE(listt[v],list2[v])) #runs euler's method for the current function and sets the next value of the function equal to h list2.append(h) #adds h to the list of current values ans= (zip(t,listt)) #creates the pairs of time and voltage to be plotted return ans #returns the ordered pairs that can be used to plot
list_plot(eulersthresh(0.4,0.8,0.1,500), plotjoined=true)
def eulerssupra(init,init2, stepsize,numtime): #sets up a function that takes in two inital values, the stepsize, and the numtime numsteps= floor(numtime/stepsize) listt = [init] list2 = [init2] functionEM(N,w)=N*(1-N)*(N-0.2)-w +0.6 #we changed the Iext value to 0.6 functionE(N,w)= .01*(N-.5*w) vals = srange(0, numsteps+1,1) t=srange(0, numsteps*stepsize+stepsize,stepsize) for v in vals: if listt[v]>=1: listt.append(init) list2.append(init2) ans= (zip(t,listt)) else: g=listt[v]+ (stepsize*functionEM(listt[v],list2[v])) listt.append(g) h= list2[v] + (stepsize*functionE(listt[v],list2[v])) list2.append(h) ans= (zip(t,listt)) return ans
list_plot(eulerssupra(0.4,0.8,.1,500), plotjoined=true)
def eulerssub(init,init2, stepsize,numtime): #sets up a function that takes in two inital values, the stepsize, and the numtime numsteps= floor(numtime/stepsize) listt = [init] list2 = [init2] functionEM(N,w)=N*(1-N)*(N-0.2)-w +0.1 #we changed the Iext value to .1 functionE(N,w)= .01*(N-.5*w) vals = srange(0, numsteps+1,1) t=srange(0, numsteps*stepsize+stepsize,stepsize) for v in vals: if listt[v]>=1: listt.append(init) list2.append(init2) ans= (zip(t,listt)) else: g=listt[v]+ (stepsize*functionEM(listt[v],list2[v])) listt.append(g) h= list2[v] + (stepsize*functionE(listt[v],list2[v])) list2.append(h) ans= (zip(t,listt)) return ans
list_plot(eulerssub(0.4,0.8,.1,500), plotjoined=true)
def eulersfreq(init,init2, stepsize,numtime): #sets up a function that takes in two inital values, the stepsize, and the numtime numsteps= floor(numtime/stepsize) listt = [init] list2 = [init2] functionEM(N,w)=N*(1-N)*(N-0.2)-w +0.6 #we changed the Iext value to 0.6 functionE(N,w)= .01*(N-.5*w) vals = srange(0, numsteps+1,1) t=srange(0, numsteps*stepsize+stepsize,stepsize) for v in vals: if listt[v]>=1: listt.append(init) list2.append(init2) ans= (zip(t,listt)) else: g=listt[v]+ (stepsize*functionEM(listt[v],list2[v])) listt.append(g) h= list2[v] + (stepsize*functionE(listt[v],list2[v])) list2.append(h) ans= listt return ans
#Frequency with Iext .6 var("N", "w") #establishes N and W as variables maxlist = [] #creates a list to hold the max values values= eulersfreq(0.4,0.8,.1,200) #runs the simulation for eulerssupra without ploting it, instead returning and setting it equal to values so it can be called later
(N, w)
length= srange(0,len(values),1) #creates a range for the for loop to run using the length of values for i in length: #creates a foor loop that will run for each value in values if values[i] > .9962: #checks to see if the current value is greater than .9962 Note:When .9962 is the approx. value returned when the max value is called. Below this is shown. Due to the amount of decimals behind each value, we picked a rounded value for it to be higher than. When using just max(values) we only recieved one value in return because of the variation in high decimal places. maxlist.append(i) #if the value is higher than that then the index is added into maxlist print maxlist #prints maxlist
[656, 657, 658, 1534, 1535, 1536]
#Note: There are more than two values, however, of the 6, 3 are concentrated around one area and the other three around another area. We concluded that one value from each area was a sufficient representation of the max point. values[1535] #checks to make sure that the values that were returned to us make sense
0.996282727155982
values[656] #checks to make sure that the values returned are similar
0.996264347916264
max(values) #We used this to change was value we used in the if statement to determine if it was a peak
0.996282727155982
#Frequency with Iext .7 def eulersfreq(init,init2, stepsize,numtime): #sets up a function that takes in two inital values, the stepsize, and the numtime numsteps= floor(numtime/stepsize) listt = [init] list2 = [init2] functionEM(N,w)=N*(1-N)*(N-0.2)-w +0.7 #we changed the Iext value to 0.7 functionE(N,w)= .01*(N-.5*w) vals = srange(0, numsteps+1,1) t=srange(0, numsteps*stepsize+stepsize,stepsize) for v in vals: if listt[v]>=1: listt.append(init) list2.append(init2) ans= (zip(t,listt)) else: g=listt[v]+ (stepsize*functionEM(listt[v],list2[v])) listt.append(g) h= list2[v] + (stepsize*functionE(listt[v],list2[v])) list2.append(h) ans= listt return ans
var("N", "w") #Declaring N and w as variables maxlist = [] #creates a list to hold the max values values= eulersfreq(0.4,0.8,.1,200) #runs the simulation for eulersfreq without ploting it, instead returning and setting it equal to values so it can be called later
(N, w)
length= srange(0,len(values),1) #creates a range for the for loop to run using the length of values for i in length: #creates a foor loop that will run for each value in values if values[i] > .9962: #checks to see if the current value is greater than .9962 Note:When .9962 is the approx. value returned when the max value is called. Below this is shown. Due to the amount of decimals behind each value, we picked a rounded value for it to be higher than. When using just max(values) we only recieved one value in return because of the variation in high decimal places. maxlist.append(i) #if the value is higher than that then the index is added into maxlist print maxlist #Print maxlist
[464, 465, 466, 931, 932, 933, 1398, 1399, 1400, 1865, 1866, 1867]
#frequency with Iext .85 def eulersfreq(init,init2, stepsize,numtime): #sets up a function that takes in two inital values, the stepsize, and the numtime numsteps= floor(numtime/stepsize) listt = [init] list2 = [init2] functionEM(N,w)=N*(1-N)*(N-0.2)-w +0.85 #we changed the Iext value to 0.85 functionE(N,w)= .01*(N-.5*w) vals = srange(0, numsteps+1,1) t=srange(0, numsteps*stepsize+stepsize,stepsize) for v in vals: if listt[v]>=1: listt.append(init) list2.append(init2) ans= (zip(t,listt)) else: g=listt[v]+ (stepsize*functionEM(listt[v],list2[v])) listt.append(g) h= list2[v] + (stepsize*functionE(listt[v],list2[v])) list2.append(h) ans= listt return ans
var("N", "w") maxlist = [] values= eulersfreq(0.4,0.8,.1,200)
(N, w)
length= srange(0,len(values),1) for i in length: if values[i] > .992: maxlist.append(i) print maxlist
[51, 52, 53, 105, 106, 107, 159, 160, 161, 213, 214, 215, 267, 268, 269, 321, 322, 323, 375, 376, 377, 429, 430, 431, 483, 484, 485, 537, 538, 539, 591, 592, 593, 645, 646, 647, 699, 700, 701, 753, 754, 755, 807, 808, 809, 861, 862, 863, 915, 916, 917, 969, 970, 971, 1023, 1024, 1025, 1077, 1078, 1079, 1131, 1132, 1133, 1185, 1186, 1187, 1239, 1240, 1241, 1293, 1294, 1295, 1347, 1348, 1349, 1401, 1402, 1403, 1455, 1456, 1457, 1509, 1510, 1511, 1563, 1564, 1565, 1617, 1618, 1619, 1671, 1672, 1673, 1725, 1726, 1727, 1779, 1780, 1781, 1833, 1834, 1835, 1887, 1888, 1889, 1941, 1942, 1943, 1995, 1996, 1997]
#frequency Iext = .9 def eulersfreq(init,init2, stepsize,numtime): #sets up a function that takes in two inital values, the stepsize, and the numtime numsteps= floor(numtime/stepsize) listt = [init] list2 = [init2] functionEM(N,w)=N*(1-N)*(N-0.2)-w +0.9 #we changed the Iext value to 0.9 functionE(N,w)= .01*(N-.5*w) vals = srange(0, numsteps+1,1) t=srange(0, numsteps*stepsize+stepsize,stepsize) for v in vals: if listt[v]>=1: listt.append(init) list2.append(init2) ans= (zip(t,listt)) else: g=listt[v]+ (stepsize*functionEM(listt[v],list2[v])) listt.append(g) h= list2[v] + (stepsize*functionE(listt[v],list2[v])) list2.append(h) ans= listt return ans
var("N", "w") maxlist = [] values= eulersfreq(0.4,0.8,.1,200)
(N, w)
length= srange(0,len(values),1) for i in length: if values[i] > .9962: maxlist.append(i) print maxlist
[35, 36, 72, 73, 109, 110, 146, 147, 183, 184, 220, 221, 257, 258, 294, 295, 331, 332, 368, 369, 405, 406, 442, 443, 479, 480, 516, 517, 553, 554, 590, 591, 627, 628, 664, 665, 701, 702, 738, 739, 775, 776, 812, 813, 849, 850, 886, 887, 923, 924, 960, 961, 997, 998, 1034, 1035, 1071, 1072, 1108, 1109, 1145, 1146, 1182, 1183, 1219, 1220, 1256, 1257, 1293, 1294, 1330, 1331, 1367, 1368, 1404, 1405, 1441, 1442, 1478, 1479, 1515, 1516, 1552, 1553, 1589, 1590, 1626, 1627, 1663, 1664, 1700, 1701, 1737, 1738, 1774, 1775, 1811, 1812, 1848, 1849, 1885, 1886, 1922, 1923, 1959, 1960, 1996, 1997]
lenvals= [656, 1536, 464, 931, 51, 105, 35, 72] #creates a list of the different max points pairs at four different frequencies g=1/(lenvals[1]-lenvals[0]) #creates the frequency values for each pair of max points for each Iext values h=1/(lenvals[3]-lenvals[2]) i=1/(lenvals[5]-lenvals[4]) j=1/(lenvals[7]-lenvals[6]) print g,h,i,j
1/880 1/467 1/54 1/37
freq=[g,h,i,j] #creates a list of the frequencies iext= [.6,.7,.85,.9] #creates a list of the Iext values points= [(.6,g),(.7,h),(.85,i),(.9,j)] #creates a list of pairs with Iext values as x and frequency values as y
list_plot(points, axes_labels=("Iext","frequency"), plotjoined=true) #plots the frequencies against the Iext values
