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 > ...
