I understand that the OP may no longer be following the thread, so my previous post is also a general question to the community at large.
On Tuesday, August 5, 2014 8:01:32 PM UTC-4, yaois...@gmail.com wrote: > > Frederick, thanks for following up with your post. > As a Julia newbie, I was just wondering where in the function you called > Sundials? I've been trying to write my own Adams Bashforth function, and > having dealt with an endless line of bugs, I think I may opt to learn how > to use an existing package instead. > > On Thursday, June 19, 2014 9:22:01 PM UTC-4, Frederick wrote: >> >> Success. To answer my own question, yes, Sundials does the extra work >> and no, the step size handed to Sundials doesn't matter. (That was a >> problem in my code.) >> >> For anyone that's interested, here is the julia code I managed to match >> the MATLAB code above. Maybe could be cleaner in some places, but at least >> it gets the job done. >> >> Again, Thanks to all for the help. (code highlighting fails halfway >> through...) >> >> >> function dydt_hh(t,statevar, statevardot) >> (V,m,h,n) = >> tuple([i for i in statevar[1:4]]...) ; >> >> (F, R, T, RTF, >> Nao, Ko, Nai, Ki, >> Cm, >> GNa, GK, Gl, ENa, EK, El, Istim) = >> tuple([i for i in statevar[5:20]]...) ; >> >> #%% compute ionic currents >> gNa = GNa*m^3*h; >> INa = gNa*(V-ENa); >> gK = GK*n^4; >> IK = gK*(V-EK); >> Il = Gl*(V-El) ; >> Iion = INa + IK + Il ; >> >> #% compute rate constants to update gates >> alphan = 0.01*(V+50)/(1-exp(-(V+50)/10)); >> betan = 0.125*exp(-(V+60)/80); >> >> alpham = 0.1*(V+35)/(1-exp(-(V+35)/10)); >> betam = 4*exp(-(V+60)/18) ; >> >> alphah = 0.07*exp(-(V+60)/20); >> betah = 1/(exp(-(V+30)/10)+1); >> >> statevardot[2] = alpham - (alpham+betam)*m ; # dm >> statevardot[3] = alphah - (alphah+betah)*h ; # dh >> statevardot[4] = alphan - (alphan+betan)*n ; # dn >> statevardot[1] = -(Istim+Iion)/Cm; # dV >> end >> >> >> # %% Hodgkin-Huxley model >> # t time ms >> # V membrane potantial mV >> # INa,IK,Il,Iion ionic current uA/cm2 >> # Cm capacitance uF/cm2 >> >> #%%%%%%%%%%%% Step 1: Define all constants >> >> # Physical constants >> #global F R T RTF >> F = 96.5; #% Faraday constant, coulombs/mmol >> R = 8.314; #% gas constant, J/K >> T_celsius = 6.3; #% Temperature in celsius >> T = 273 + T_celsius ; #% absolute temperature, K >> >> RTF = R*T/F ; >> >> # default concentrations for squid axon in sea water - mmol/l >> #global Nao Ko Nai Ki >> Nao = 491 ; >> Ko = 20 ; >> Nai = 50 ; >> Ki = 400 ; >> >> #% Cell constant >> #global Cm >> Cm = 1 ; #% membrane capacitance, uF/cm^2; >> >> #% Maximum channel conductances -- mS/cm^2 >> #global GNa GK Gl ENa EK El >> GNa = 120; >> GK = 36; >> Gl = 0.3; >> >> #% Nernst potentials -- mV >> ENa = RTF*log(Nao/Nai); >> EK = RTF*log(Ko/Ki); >> El = -49; >> >> #%%%%%%%%%%% Step 2: Define simulation and stimulus parameters >> >> t_end = 50 ; #% end of simulation, ms >> >> stim_delay = 5 ; >> stim_dur = 3 ; >> stim_amp = -50 ; >> >> stim_start = stim_delay ; >> stim_end = stim_delay + stim_dur ; >> >> # % % Intervals defined as follows >> # % % 1) start of simulation (t = 0) to beginning of stimulus (t = >> stim_start) >> # % % 2) beginning (t = stim_start) to end of stimulus (t = stim_end) >> # % % 3) end of stimulus (t = stim_end) to end of simulation (t = t_end) >> stimints = 3 ; >> intervals = zeros(3,2) >> intervals[1,:] = [0,stim_start] ; # first row >> intervals[2,:] = [stim_start,stim_end]<span style="color: #000 >> ... > >
