% % % solve the 2 compartment lgmd model % for a constant current injection % % The model contains % % 1) I_Na, I_DR in the axon compartment with variables m_naa, h_nai, m_dra % % 2) I_Ca, I_AHP in the axon compartment with variables m_caa, i_ahp, ca_c % % % the state variable vector is: % % y = [va vd h_nai m_dra ca_c] % % where va = axon membrane potential, vd = dendritic membrane potential, % h_nai = sodium conductance inactivation variable, m_dra = delayed % rectifier activation variable, ca_c = calcium concentration. % All other variables are set to their instantaneous steady-state value. % % % usage three_cmpt(T,[t_s t_e Iv cmpt],param_struct) e.g., lgmd_mod(500,[50 450 1],p_s); % % where param_struct is a structure containing the parameters of the % model. Current injection parameters are determined as follows: % % T = duration of simulation (ms) % t_s, t_e = start and end time of current stimulation, % respectively % Iv = current injection value (muA/cm2) % cmpt = compartment (1 = axon 2 = dendrite 3 = mid compartment)* % % * YES, the wiring *IS* weird -- 2 <-> 3 <-> 1 is the scheme. % % The parameter structure (p_s) is organized as follows: % % p_s.v_ss = steady-state vm for initialization purposes % p_s.... = other variables of the model % % figure produces Vs and Vd vs. time % function [t, y] = three_cmpt(ps_mod) %initial steady state values v_ss = ps_mod.v_ss; cacss = ca_c_ss(v_ss, ps_mod.alphaCa, ps_mod.tauCa, ps_mod.VCa, ps_mod.gCa, ps_mod.gCaV12, ps_mod.gCaSl); %cacss = 0; % force start [ca] = 0 y0 = [v_ss v_ss v_ss m_naa_ss(v_ss) h_nai_ss(v_ss) m_dra_ss(v_ss) cacss m_ih_ss(v_ss) m_caa_ss(v_ss, ps_mod.gCaV12, ps_mod.gCaSl), 0, 0, 0]; % an interesting plot if ( 0 == 1 ) v = -100:1:100; for i=1:length(v) cass_a(i) = ca_c_ss(v(i), ps_mod.alphaCa, ps_mod.tauCa, ps_mod.VCa, ps_mod.gCa, ps_mod.gCaV12, ps_mod.gCaSl); cass_b(i) = ca_c_ss(v(i), 2*ps_mod.alphaCa, ps_mod.tauCa, ps_mod.VCa, ps_mod.gCa, ps_mod.gCaV12, ps_mod.gCaSl); cass_c(i) = ca_c_ss(v(i), .5*ps_mod.alphaCa, ps_mod.tauCa, ps_mod.VCa, ps_mod.gCa, ps_mod.gCaV12, ps_mod.gCaSl); cass_d(i) = ca_c_ss(v(i), ps_mod.alphaCa, ps_mod.tauCa, ps_mod.VCa, ps_mod.gCa, ps_mod.gCaV12, 2*ps_mod.gCaSl); cass_e(i) = ca_c_ss(v(i), ps_mod.alphaCa, ps_mod.tauCa, ps_mod.VCa, ps_mod.gCa, ps_mod.gCaV12, .5*ps_mod.gCaSl); % and "m_ahp(ca)" mahp_a(i) = cass_a(i)/(cass_a(i) + ps_mod.kD); mahp_b(i) = cass_a(i)/(cass_a(i) + 2*ps_mod.kD); mahp_c(i) = cass_a(i)/(cass_a(i) + 4*ps_mod.kD); mahp_d(i) = cass_b(i)/(cass_b(i) + ps_mod.kD); mahp_e(i) = cass_c(i)/(cass_c(i) + ps_mod.kD); end subplot(2,1,1); plot (v,cass_a, 'k'); hold on; plot (v,cass_b, 'rx'); plot (v,cass_c, 'gx'); plot (v,cass_d, 'm'); plot (v,cass_e, 'b'); subplot(2,1,2); plot (v,mahp_a, 'k'); hold on; plot (v,mahp_b, 'rx'); plot (v,mahp_c, 'gx'); plot (v,mahp_d, 'm'); plot (v,mahp_e, 'b'); pause end % Setup the current vectors -- if necessary if (isfield(ps_mod, 'I_of_t') == 0) ps_mod.I_of_t = zeros(3,length(0:ps_mod.dt:ps_mod.duration)); ps_mod.I_of_t(3,round(ps_mod.I_inj_start/ps_mod.dt):round(ps_mod.I_inj_end/ps_mod.dt)) = ps_mod.I_inj_nA(3)*.001/ps_mod.Am; % convert to uA/cm2 ps_mod.I_of_t(2,round(ps_mod.I_inj_start/ps_mod.dt):round(ps_mod.I_inj_end/ps_mod.dt)) = ps_mod.I_inj_nA(2)*.001/ps_mod.Ad; % convert to uA/cm2 ps_mod.I_of_t(1,round(ps_mod.I_inj_start/ps_mod.dt):round(ps_mod.I_inj_end/ps_mod.dt)) = ps_mod.I_inj_nA(1)*.001/ps_mod.Aa; end if (0 == 1) figure; plot(0:ps_mod.dt:ps_mod.duration, ps_mod.I_of_t(1,:), 'b-'); hold on; plot(0:ps_mod.dt:ps_mod.duration, ps_mod.I_of_t(2,:), 'r-'); plot(0:ps_mod.dt:ps_mod.duration, ps_mod.I_of_t(3,:), 'k-'); pause; end % Synaptic noise g_syn_of_t_exc_spont = get_synaptic_noise(ps_mod.n_syns_per_facet_exc,... ps_mod.spontaneous_firing_rate_exc,ps_mod.n_minis_avg_exc, ps_mod.tau_syn_exc, ... ps_mod.g_max_spont_exc, ps_mod.duration, ps_mod.dt); g_syn_of_t_inh_spont = get_synaptic_noise(ps_mod.n_syns_per_facet_inh,... ps_mod.spontaneous_firing_rate_inh,ps_mod.n_minis_avg_inh, ps_mod.tau_syn_inh, ... ps_mod.g_max_spont_inh, ps_mod.duration, ps_mod.dt); ps_mod.g_syn_of_t_exc = ps_mod.g_syn_of_t_exc + g_syn_of_t_exc_spont; ps_mod.g_syn_of_t_inh = ps_mod.g_syn_of_t_inh + g_syn_of_t_inh_spont; % For debug if ( 0 == 1 ) figure; plot(0:ps_mod.dt:ps_mod.duration, ps_mod.g_syn_of_t_exc', 'r'); hold on; end % Note: ode45 does not converge, probably because the equations are stiff. % ode15s and ode23 seem to give similar solutions disp('starting simulation ...'); tic; [t,y] = ode23(@red,0:ps_mod.dt:ps_mod.duration,y0,[],ps_mod); %options = [] lensec = toc; disp([num2str(lensec) ' seconds elapsed; ' num2str(1000*(lensec/(ps_mod.duration/ps_mod.dt))) ' ms/step']); % save? if (isfield(ps_mod, 'savefile') ~= 0) save(ps_mod.savefile, 't', 'y', 'ps_mod'); end return % % reduced wang model % function dy = red(t,y,ps_mod) Cm = ps_mod.Cm; %capacitance, muF/cm2 dy = zeros(12,1); % --- current injection Ia = ps_mod.I_of_t(1,1+round(t/ps_mod.dt)); Id = ps_mod.I_of_t(2,1+round(t/ps_mod.dt)); Im = ps_mod.I_of_t(3,1+round(t/ps_mod.dt)); % ---- axon compartment v - leak, Na HH, Kdr HH IL_a = ps_mod.gLa*(y(1)-ps_mod.VL); INa = ps_mod.gNa*y(4)^3*y(5)*(y(1)-ps_mod.VNa); IKDR = ps_mod.gKDR*y(6)^4*(y(1)-ps_mod.VK); dy(1) = (-IL_a - INa -IKDR - (ps_mod.gc_ma)*(y(1)-y(3)) + Ia)/ps_mod.Cm; % --- dendritic compartment v - syn, inj, leak, H IL_d = ps_mod.gLd*(y(2)-ps_mod.VL); % NO I_H for now IH = 0; dy(2) = (-IL_d -IH - (ps_mod.gc_md)*(y(2)-y(3)) + Id)/ps_mod.Cm; % --- Middle compartment v - leak, IL_m = ps_mod.gLm*(y(3)-ps_mod.VL); dy(3) = (-IL_m -(ps_mod.gc_am)*(y(3)-y(1)) - (ps_mod.gc_dm)*(y(3)-y(2)) + Im)/ps_mod.Cm; % --- synaptic current Isyn_exc = ps_mod.g_syn_of_t_exc(1+round(t/ps_mod.dt))*(y(ps_mod.exc_syn_cmpt)-ps_mod.erev_syn_exc); dy(ps_mod.exc_syn_cmpt) = dy(ps_mod.exc_syn_cmpt) - Isyn_exc/ps_mod.Cm; Isyn_inh = ps_mod.g_syn_of_t_inh(1+round(t/ps_mod.dt))*(y(ps_mod.inh_syn_cmpt)-ps_mod.erev_syn_inh); dy(ps_mod.inh_syn_cmpt) = dy(ps_mod.inh_syn_cmpt) - Isyn_inh/ps_mod.Cm; % --- Calcium compartment ICa = ps_mod.gCa*y(9)*(y(ps_mod.ca_cmpt)-ps_mod.VCa); ca_av = ps_mod.sf_bapta*y(7); IKAHP = ps_mod.gKAHP*(ca_av/(ca_av+ps_mod.kD))*(y(ps_mod.ca_cmpt)-ps_mod.VK); dy(ps_mod.ca_cmpt) = dy(ps_mod.ca_cmpt) + (-ICa - IKAHP)/ps_mod.Cm; %m_naa phi_m = 4; %temperature factor dy(4) = phi_m*(m_naa_ss(y(1))-y(4))/(ps_mod.tau_naa_sf/(a_naa(y(1)) + b_naa(y(1)))); %h_nai phi_n = 4; %temperature factor dy(5) = phi_m*(h_nai_ss(y(1))-y(5))/(ps_mod.tau_nai_sf/(a_nai(y(1)) + b_nai(y(1)))); %m_dra phi_h = 4; %temperature factor dy(6) = phi_h*(m_dra_ss(y(1))-y(6))/(ps_mod.tau_kdr_sf/(a_dra(y(1)) + b_dra(y(1)))); %calcium concentration (FREELY AVAILABLE!) dy(7) = -ps_mod.alphaCa*ICa - y(7)/ps_mod.tauCa ; % IH gate dy(8) = phi_m*(m_ih_ss(y(2))-y(8))/tau_ih(y(2)); % Ca gate dy(9) = phi_m*(m_caa_ss(y(ps_mod.ca_cmpt), ps_mod.gCaV12, ps_mod.gCaSl)-y(9))/tau_caa(y(ps_mod.ca_cmpt), ps_mod.gCaTau); % bound bapta - DO NOT USE THIS dy(10) = 0; % Ca bound to AHP gate; DO NOT USE dy(11) = 0; % Ca bound to internal buffers - these should have a low concentration and get overwhelmed, but still delay onset of ICa - UNUSED dy(12) = 0; return %steady-state sodium activation function val = m_naa_ss(v) val = a_naa(v)./(a_naa(v) + b_naa(v)); % sodium activation forward rate constant function val = a_naa(v) if ( v == -33 ) %a_naa is not defined there val = 1; else val= -0.1*(v+33)./(exp(-0.1*(v+33))-1); end % sodium activation backward rate constant function val = b_naa(v) val = 4*exp(-(v+60)/12); %steady-state sodium inactivation function val = h_nai_ss(v) val = a_nai(v)./(a_nai(v) + b_nai(v)); %sodium inactivation variable forward rate constant function val = a_nai(v) val = 0.07*exp(-0.1*(v+50)); %sodium inactivation variable backward rate constant function val = b_nai(v) val = 1/(exp(-0.1*(v+20)) + 1); %steady-state delayed-rectifier activation function val = m_dra_ss(v) val = a_dra(v)./(a_dra(v) + b_dra(v)); %delayed rectifier activation variable forward rate constant function val = a_dra(v) if ( v == -34 ) val = 0.1; else val = -0.01*(v+34)./(exp(-0.1*(v+34))-1); end; %delayed rectifier activation variable backward rate constant function val = b_dra(v) val = 0.125*exp(-(v+44)/25); %steady-state calcium conductance activation function val = m_caa_ss(v, v12, sl) val = 1./(1+exp((v-v12)/sl)); % Voltage % very fast, but not instantaneous function val = tau_caa(v, tau) %val = tau; val = (0.25)*tau + (0.75*tau)./(1+exp((v+10)/-5)); % Voltage %steady-state calcium concentration function val = ca_c_ss(v, alphaCa, tauCa, VCa, gCa, Cav12, CaSl) %steady-state calcium current ICa_ss = gCa*m_caa_ss(v, Cav12, CaSl)*(v-VCa); %steady-state Ca concentration val = -1*tauCa*alphaCa*ICa_ss; % steady-state Ih activation function val = m_ih_ss(v) val = 1./(1+exp((v+92)/27.4)); % time constant for ih function val = tau_ih(v) val = 1./(.0018*exp((v+86.6)/-47.5) + .0097*exp((v+85.2)/17.7));