Bursting and oscillations in RD1 Retina driven by AII Amacrine Neuron (Choi et al. 2014)

 Download zip file 
Help downloading and running models
Accession:156781
"In many forms of retinal degeneration, photoreceptors die but inner retinal circuits remain intact. In the rd1 mouse, an established model for blinding retinal diseases, spontaneous activity in the coupled network of AII amacrine and ON cone bipolar cells leads to rhythmic bursting of ganglion cells. Since such activity could impair retinal and/or cortical responses to restored photoreceptor function, understanding its nature is important for developing treatments of retinal pathologies. Here we analyzed a compartmental model of the wild-type mouse AII amacrine cell to predict that the cell's intrinsic membrane properties, specifically, interacting fast Na and slow, M-type K conductances, would allow its membrane potential to oscillate when light-evoked excitatory synaptic inputs were withdrawn following photoreceptor degeneration. ..."
Reference:
1 . Choi H, Zhang L, Cembrowski MS, Sabottke CF, Markowitz AL, Butts DA, Kath WL, Singer JH, Riecke H (2014) Intrinsic bursting of AII amacrine cells underlies oscillations in the rd1 mouse retina. J Neurophysiol 112:1491-504 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Retina;
Cell Type(s):
Channel(s): I Na,t; I A; I M;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: MATLAB;
Model Concept(s): Bursting; Oscillations; Action Potentials; Bifurcation;
Implementer(s):
Search NeuronDB for information about:  I Na,t; I A; I M;
/
ChoiEtAl2014
readme.txt
A2ONCB.m
A2ONCB_FHN.m
A2ONCB_Fig5a.m
                            
function A2ONCB
%%% AII-ONCB network for a single AII and a single ONCB
%%% Cell 1 is AII; Cell
% Units: mF, mV,ms,S,ohm,mA, cm

clc
clear all

global numcell armsection R capacitance                  % Number of cells and compartments, R and capacitance
global gbar_na gbar_M gbar_A G_gap Gpas epas ena ek      % conductances
global vhalfh_na vhalfm_na kh_na km_na                   % Na current activation/inactivation
global vhalfm_M km_M                                     % M-type K (slow) current activation
global vhalfh_A vhalfm_A vhalfw_A kh_A km_A kw_A f       % A-type K (fast) current activation/inactivation
global mtau_na htau_na mtau_M mtau_A                     % Time constants
global IinjectAII IinjectONCB                                              % Injected current
global epas2

% set the values of injected currents with InjectVector (a list of values indicates switching to other values during the run)

% Note that the plots only include data after a transient given by cuttime.

% Set the conductance of the gap junction between AII and ONCB using the switch ONCB:
% ONCB=1 for Fig.14 and for Fig.5B top. ONCB=2,3,4 for remaining panels of Fig.5B


%%%%%%%%%% -------------------------- Switches -------------------------------------%%%%%%%%

plotswitch = 14;  % 14 for Fig 14, 5 for Fig 5
scan = 0; % 0 if looking at voltage traces; 1 if scanning Vm (Fig2A); 2 if scanning VhlafM (Fig3B)
ONCB =1; % ONCB type (coupling strength)

RD = 1; % 1 if RD; 0 if wildtype
ONCBnum = 1;% 5; %

duration=1000; % duration of a run for a given parameter value

IinjectAII = 0;
IinjectONCB = 0;

% Current injection into ONCB
if plotswitch == 5
    InjectVector = [0];
elseif plotswitch == 14
    InjectVector = [0:2e-9:20*1e-9];
end

Vector = InjectVector;
tol_steady = 0.1; %(ms)

lookcell = 2; % 1 for cell 1 (AII), 2 for cell 2 (ONCB)
numcell = 2;

plotsection =1; %1 is Soma; 2 is cable; 3 is IS

v_init =  -62;%(mV)
v_init2 = -50;%(mV)

% set coupling between AII and ONCB
if ONCB == 1
    gjCondR=  750*1e-12 *ONCBnum; %(S)
elseif ONCB == 2
    gjCondR=  300*1e-12 *ONCBnum; %(S)
elseif ONCB == 3;
    gjCondR=  150*1e-12 *ONCBnum;
else
    gjCondR=  100*1e-12 *ONCBnum; %(S)
end

gjCondIS= 1e-5*1e-12; %(S)
gjCondDC= 1e-5*1e-12; %(S)

%%%%%%%%%% -------------------------- Parameters Unchanged --------------------------------%%%%%%%%
armsection = 1; % # of sections in the cable, 11 in Mark's code

somadiam=25*1e-4;  %(cm)
somalength=25*1e-4;  %(cm)
armdiam=0.3*1e-4;  %(cm)
armlength=32*1e-4;  %(cm)
handdiam=2*1e-4;  %(cm)
handlength=2*1e-4;  %(cm)

SA_soma= 1963.5*1e-8; %(cm^2)
SA_arm=30.1593*1e-8; %(cm^2)
SA_hand=12.5664*1e-8; %(cm^2)

SA_ONCB= 439.8*1e-8; %(cm^2)

global_ra = 150; %(ohm-cm)

ena = 50; %(mV)
ek = -77; %(mV)

Cm = 1e-3; %(mf/cm^2)


%%
%%% --------- Na and K currents activation/inactivation -------- %%%
vhalfh_na = -49.5; %(mV)
vhalfm_na = -48; %(mV)
kh_na = 2; %(mV)
km_na = 5; %(mV)
mtau_na = 0.01; %(ms)
htau_na = 0.5; %(ms)

vhalfm_M = -40; %(mV)
km_M = 4; %(mV)
mtau_M = 50;%(ms)

vhalfh_A = -40.5; %(mV)
vhalfm_A = -10; %(mV)
vhalfw_A= -45; %(mV)
kh_A = 2; %(mV)
km_A = 7; %(mV)
kw_A = 15; %(mV)
f=0.83;
mtau_A = 1; %(ms)

gbar_na_soma = 0;
gbar_na_arm = 0;
gbar_na_IS = 0.2; %(mho/cm^2)

gbar_M_soma = 0;
gbar_M_arm = 0;
gbar_M_IS = 0.03; %(mho/cm2)

gbar_A_soma = 0.004; %(mho/cm2)
gbar_A_arm = 0;
gbar_A_IS = 0.08; %(mho/cm2)



if numcell == 2
    
    gbar_na_soma2 = 0;
    gbar_na_arm2 = 0;
    gbar_na_IS2 = 0; %(mho/cm^2)
    
    gbar_M_soma2 = 0;
    gbar_M_arm2 = 0;
    gbar_M_IS2 = 0;  %(mho/cm2)
    
    gbar_A_soma2 = 0;  %(mho/cm2)
    gbar_A_arm2 = 0;
    gbar_A_IS2 = 0;  %(mho/cm2)
    
end


%%%%%%%% ----------------------------- Leak current for AII ------------------------------------ %%%%%%%

if RD == 0
    Rm = 40000; %(ohm*cm^2)
    epas = -40; %(mV)
    gpas=1/Rm; %(S/cm^2)
    
    Rm2 = 12000;
    epas2 = -30;%-50;
    gpas2=1/Rm2;
    
elseif RD == 1
    
    Rm = 40000; %(ohm*cm^2)
    epas = -65;%-30; %(mV)
    gpas=1/Rm; %(S/cm^2)
    
    Rm2 = 12000;
    epas2 = -35;%-50;
    gpas2=1/Rm2;
    
end

%%%%%%%% --------------------------  Resistance between sections -------------------------------- %%%%%%%
R=zeros(armsection+1,1);
for num=1:armsection+1
    if num==1
        R(num)=(global_ra/(pi*(somadiam/2)^2))*somalength/2 + (global_ra/(pi*(armdiam/2)^2))*armlength/armsection/2;
    elseif num==armsection+1
        R(num)=(global_ra/(pi*(handdiam/2)^2))*handlength/2 + (global_ra/(pi*(armdiam/2)^2))*armlength/armsection/2;
    else
        R(num)=(global_ra/(pi*(armdiam/2)^2))*armlength/armsection;
    end
end



%%%%%%%% ------------------------------ Conductances Set Up -------------------------%%%%%%%%
gbar_na=zeros(numcell,armsection+2);
gbar_M=zeros(numcell,armsection+2);
gbar_A=zeros(numcell,armsection+2);
capacitance=zeros(numcell,armsection+2);
G_gap=zeros(numcell,armsection+2);
Gpas=zeros(numcell,armsection+2);



for ii=1:numcell
    
    if ii == 2
        gbar_na_soma = gbar_na_soma2;
        gbar_na_arm = gbar_na_arm2;
        gbar_na_IS = gbar_na_IS2; %(mho/cm^2)
        
        gbar_M_soma = gbar_M_soma2;
        gbar_M_arm = gbar_M_arm2;
        gbar_M_IS = gbar_M_IS2;  %(mho/cm2)
        
        gbar_A_soma = gbar_A_soma2;  %(mho/cm2)
        gbar_A_arm = gbar_A_arm2;
        gbar_A_IS = gbar_A_IS2;  %(mho/cm2)
        
        gpas = gpas2;
        SA_soma = SA_ONCB;
    end
    
    
    for jj=1:armsection+2
        if jj==1
            gbar_na(ii,jj) = gbar_na_soma * SA_soma;
            gbar_M(ii,jj) = gbar_M_soma * SA_soma;
            gbar_A(ii,jj)= gbar_A_soma * SA_soma;
            capacitance(ii,jj) = Cm*SA_soma;
            G_gap(ii,jj) = gjCondR;
            Gpas(ii,jj) = gpas * SA_soma;
            
            
        elseif jj==armsection+2
            gbar_na(ii,jj) = gbar_na_IS * SA_hand;
            gbar_M(ii,jj) = gbar_M_IS * SA_hand;
            gbar_A(ii,jj) = gbar_A_IS * SA_hand;
            capacitance(ii,jj) = Cm*SA_hand;
            G_gap(ii,jj) = 0;%gjCondIS;
            
            if ii == 1
                Gpas(ii,jj) = gpas * SA_hand;
            elseif ii == 2
                Gpas(ii,jj) = 0;
            end
            
            
        else
            gbar_na(ii,jj) = gbar_na_arm * SA_arm/armsection;
            gbar_M(ii,jj) = gbar_M_arm * SA_arm/armsection;
            gbar_A(ii,jj) = gbar_A_arm * SA_arm/armsection;
            capacitance(ii,jj) = Cm*SA_arm/armsection;
            G_gap(ii,jj) = 0;%gjCondDC;
            
            if ii == 1
                Gpas(ii,jj) = gpas * SA_arm/armsection;
            elseif ii == 2
                Gpas(ii,jj) = 0;
            end
            
        end
        
    end
end


%%%%%% ------------------- Initial conditions --------------- %%%%%%%%%%

%%--------------- Initial conditions for Na
minit_na = 1/(1+exp(-(v_init-vhalfm_na)/km_na));
hinit_na = 1/(1+exp((v_init-vhalfh_na)/kh_na));

%%--------------- Initial conditions for M-type K
minit_M = 1/(1 + exp(-(v_init-vhalfm_M)/km_M));
%%--------------- Initial conditions for A-type K
minit_A = 1/(1+exp(-(v_init-vhalfm_A)/km_A));
hinit_A = f*(1/(1 + exp((v_init-vhalfh_A)/kh_A)))+(1-f);

%%%%%% ------------------- Initial conditions for cell 2 --------------- %%%%%%%%%%
if numcell == 2
    minit_na2 = 1/(1+exp(-(v_init2-vhalfm_na)/km_na));
    hinit_na2 = 1/(1+exp((v_init2-vhalfh_na)/kh_na));
    
    minit_M2 = 1/(1 + exp(-(v_init2-vhalfm_M)/km_M));
    
    minit_A2 = 1/(1+exp(-(v_init2-vhalfm_A)/km_A));
    hinit_A2 = f*(1/(1 + exp((v_init2-vhalfh_A)/kh_A)))+(1-f);
end

%%%%%% ------------------- Initial conditions Set Up  ------------------ %%%%%%%%%%
initial = zeros(numcell, armsection+2, 7);
for ii = 1:numcell
    if ii==1
        for jj = 1:armsection+2
            initial(ii,jj,1) = minit_na;
            initial(ii,jj,2) = hinit_na;
            initial(ii,jj,3) = minit_M;
            initial(ii,jj,4) = minit_A;
            initial(ii,jj,5) = hinit_A;
            initial(ii,jj,6) = hinit_A;
            initial(ii,jj,7) = v_init;
        end
    elseif ii==2
        for jj = 1:armsection+2
            initial(ii,jj,1) = minit_na2;
            initial(ii,jj,2) = hinit_na2;
            initial(ii,jj,3) = minit_M2;
            initial(ii,jj,4) = minit_A2;
            initial(ii,jj,5) = hinit_A2;
            initial(ii,jj,6) = hinit_A2;
            initial(ii,jj,7) = v_init2;
        end
    end
end
initial = reshape(initial,numcell*(armsection+2)*7,1);


vmaxmat = zeros(1,length(Vector ));%[];
vminmat = zeros(1,length(Vector ));%[];
frequencymat = zeros(1,length(Vector ));%[];
ampmat = zeros(1,length(Vector));

Tend = 0;

for ii = 1:length(InjectVector)
    
    IinjectONCB = InjectVector(ii);
    
    time_start = 0;
    time_end = duration;%000;
    limit_time = 1100;
    cuttime = 500;
    
    
    
    freq = 0;
    amp = 0;
    
    while time_end < limit_time
        
        %%%%%%%%%%% ------------------- Solving Diff Eqs ------------------------ %%%%%%%%%%
        
        options = odeset('RelTol', 1e-8);
        [T,Y1] = ode45(@solve_A2ONCB, [time_start time_end], initial,options);
        Y = zeros(numcell,armsection+2,7,length(T));
        for i = 1:length(T)
            Y(1:numcell,1:armsection+2,1:7,i) = reshape(Y1(i,:),numcell, armsection+2,7);
        end
        initial_new = reshape(Y(:,:,:,end),numcell*(armsection+2)*7,1);
        initial = initial_new;
        
        plotcell =1;
        v= Y(plotcell,plotsection,7,:);
        v_fix=reshape(v,size(T));
        
        vIS= Y(plotcell,3,7,:);
        vIS_fix=reshape(vIS,size(T));
        
        %%%%% -------------- If there are more than one cell -------------- %%%%%%%
        if numcell == 2
            plotcell2 = 2;
            v2= Y(plotcell2,plotsection,7,:);
            v2_fix=reshape(v2,size(T));
        end
        
        
        v_AII = v_fix;%[v_AII;v_fix];
        v_ONCB = v2_fix;%[v_ONCB;v2_fix];
        
        %%% ----------------- Cutting of the transient ------------------%%%
        savenow = [];
        for i = 1:length(T)
            if T(i) >= cuttime
                savenow = [savenow i];
            else
                savenow = savenow;
            end
            start = min(savenow);
            
        end
        
        if lookcell == 1
            v = v_fix(start:end);
        elseif lookcell == 2
            v = v2_fix(start:end);
        end
        
        t = T(start:end);
        
        
        %%% ----------------- Finding local max/minima ----------------------%%%
        if lookcell == 1
            [vmax,imax,vmin,imin] = extrema(v_fix(start:end));
        elseif lookcell == 2
            [vmax,imax,vmin,imin] = extrema(v2_fix(start:end));
        end
        
        
        imax = sort(imax);
        imin = sort(imin);
        
        
        
        %% No spiking nor bursting
        if length(vmax)<=2 || length(vmax)<=2
            
            time_end = time_end+500;
            time_start = time_start+500;
            cuttime = time_start;
            
            if time_end>=limit_time
                tmax = T(imax);
                Vmax = vmax;
                tmin = T(imin);
                Vmin = vmin;
                if length(vmax)<1 || length(vmax)<1
                    freq = 0; %(Hz)
                    amp = 0;
                    t_cross = [];
                else
                    freq = 0; %(Hz)
                    amp = max(Vmax)-min(Vmin);
                    t_cross = [];
                end
                break
            end
        else
            
            %% Quadratic interpolation
            
            tmax = zeros(1,length(imax));
            Vmax = zeros(1,length(imax));
            iMax = zeros(1,length(imax));
            
            for index = 1:length(imax)
                iMax = imax(index);
                x = [t(iMax-1) t(iMax) t(iMax+1)];
                y = [v(iMax-1) v(iMax) v(iMax+1)];
                
                [p,S,mu] = polyfit(x,y,2);
                A_new = p(1); B_new = p(2); C_new = p(3); mu1 = mu(1); mu2 = mu(2);
                c=A_new*(mu1)^2/(mu2)^2-(B_new*(mu1/mu2))+C_new;
                b = (B_new*mu2 - 2*A_new*mu1)/(mu2)^2;
                a = A_new/(mu2)^2;
                
                
                tmax(index) = -b/(2*a);
                Vmax(index) = c-b^2/(4*a);
            end
            
            tmin = zeros(1,length(imin));
            Vmin = zeros(1,length(imin));
            iMin = zeros(1,length(imin));
            
            for index = 1:length(imin)
                iMin = imin(index);
                x = [t(iMin-1) t(iMin) t(iMin+1)];
                y = [v(iMin-1) v(iMin) v(iMin+1)];
                
                [p,S,mu] = polyfit(x,y,2);
                A_new = p(1); B_new = p(2); C_new = p(3); mu1 = mu(1); mu2 = mu(2);
                c=A_new*(mu1)^2/(mu2)^2-(B_new*(mu1/mu2))+C_new;
                b = (B_new*mu2 - 2*A_new*mu1)/(mu2)^2;
                a = A_new/(mu2)^2;
                
                tmin(index) = -b/(2*a);
                Vmin(index) = c-b^2/(4*a);
            end
            
            %%
            crossline = (max(Vmax)+min(Vmin))/2;
            count1 = 1;
            count2 = 1;
            Vmin_re = [];
            while count1 <= length(tmin)
                if Vmin(count1)<crossline
                    Vmin_re(count2) = Vmin(count1);
                    count2=count2+1;
                end
                count1=count1+1;
            end
            
            count1 = 1;
            count2 = 1;
            Vmax_re = [];
            while count1 <= length(tmax)
                if Vmax(count1)>crossline
                    Vmax_re(count2) = Vmax(count1);
                    count2=count2+1;
                end
                count1=count1+1;
            end
            
            Vmax = Vmax_re;
            Vmin = Vmin_re;
            
            %%
            t_cross = [];
            v_cross = (max(vmax)+min(vmin))/2;
            for iter = 1:length(v)-1
                if (v(iter)-v_cross)<=0 && (v(iter+1)-v_cross)>=0
                    t_cross = [t_cross ((v_cross-v(iter))*(t(iter+1)-t(iter))/(v(iter+1)-v(iter)))+t(iter)];
                else
                    t_cross = t_cross;
                end
            end
            
            V_cross = ones(1,length(t_cross))*v_cross;
            
            if length(t_cross)<=2
                time_end = time_end+500;
                time_start = time_start+500;
                cuttime = time_start;
                % cuttime = cuttime+500;
                
                if time_end>=limit_time
                    
                    if length(t_cross)<=1
                        freq = 0; %(Hz)
                        amp = Vmax(end)-Vmin(end);
                    else
                        period = (t_cross(2)-t_cross(1))*1e-3;
                        freq=1/period;
                        amp = Vmax(end)-Vmin(end);
                    end
                    
                    break
                end
                
            else
                
                dtcross = zeros(1,length(t_cross)-1);
                for step = 1:length(t_cross)-1
                    dtcross(step) = t_cross(step+1)-t_cross(step);
                end
                
                vari = zeros(1,length(dtcross)-1);
                for step = 1:length(dtcross)-1
                    vari(step) = abs(dtcross(step+1)-dtcross(step));
                    if vari(step) < tol_steady
                        % cuttime = tmax(step);
                        period = dtcross(step+1)*1e-3; %(s)
                        freq = 1/period; %(Hz)
                        amp = Vmax(step+1)-Vmin(step+1);
                        break
                    end
                    
                end
                
                
                
                %% Increase the simulataion time if :
                % 1) it does not go to a steady state
                % 2) it does not have more than one period in it
                % Need a limit
                
                if min(vari)>=tol_steady
                    time_end = time_end+500;
                    time_start = time_start+500;
                    cuttime = time_start;
                    
                    if time_end>=limit_time
                        period = dtcross(end)*1e-3; %(s)
                        freq = 1/period; %(Hz)
                        amp = Vmax(end)-Vmin(end);
                        break
                    end
                else
                    break
                end
                
            end
            
        end
    end
    
    frequencymat(ii) = freq;
    ampmat(ii) = amp;
    
    
    %%%%%%%%%%------------------------------------  Decide on Plots -------------------------------------%%%%%%%%%%
    if plotswitch == 14
        figure(1)
        hold on
        subplot(3,1,1), hold on, plot(T+Tend,IinjectONCB*ones(1,length(T)))
        xlabel('Time (ms)','FontWeight','bold'), ylabel('I_{inject} (mA)','FontWeight','bold')
        subplot(3,1,2), hold on, plot(T+Tend,v_AII)
        xlabel('Time (ms)','FontWeight','bold'), ylabel('V_{AII} (mV)','FontWeight','bold')
        subplot(3,1,3), hold on, plot(T+Tend,v_ONCB)
        xlabel('Time (ms)','FontWeight','bold'), ylabel('V_{ONCB} (mV)','FontWeight','bold')       
    elseif plotswitch == 5
        figure(2)
        hold on, plot(T+Tend,v_AII,'b', T+Tend,v_ONCB,'r'), ylim([-68 -40]), xlim([0 900])
        xlabel('Time (ms)'); ylabel(' Voltage (mV)');
        legend('AII','ONCB')
    end
    
    
    
    Tend=Tend+T(end);
end

end

%%%%%%%%%%%%----------- MATLAB reconstruction of Mark's AII model on NEURON --------------%%%%%%%%%%
% Units: mF, mV,ms,S,ohm,mA, cm

function dy = solve_A2ONCB(t,y1)

global numcell armsection R capacitance                  % Number of cells and compartments, R and capacitance
global gbar_na gbar_M gbar_A G_gap Gpas epas ena ek      % conductances

global vhalfh_na vhalfm_na kh_na km_na                   % Na current activation/inactivation
global vhalfm_M km_M                                     % M-type K (slow) current activation
global vhalfh_A vhalfm_A vhalfw_A kh_A km_A kw_A f       % A-type K (fast) current activation/inactivation
global mtau_na htau_na mtau_M mtau_A                     % Time constants
global IinjectAII IinjectONCB     
global epas2
% global inject_start inject_end

y = reshape(y1,numcell,armsection+2,7);
dy = zeros(numcell,armsection+2,7);


%%%%%%%%%%%% ---------------------------- Differential Equations ------------------------- %%%%%%%%%%%
for i=1:numcell

    for j=1:armsection+2
        
        if i ==1 && j == 1% && t>=inject_start && t<inject_end% i == 17
            ie =  IinjectAII;
        elseif i == 2 && j == 1
             ie = IinjectONCB;
        else
            ie = 0;
        end

        
        dy(i,j,1) = (1/mtau_na)*(1/(1+exp(-(y(i,j,7)-vhalfm_na)/km_na))-y(i,j,1));  % m_na
        dy(i,j,2) = (1/htau_na)*(1/(1+exp((y(i,j,7)-vhalfh_na)/kh_na))-y(i,j,2));  %h_na
        dy(i,j,3) = (1/mtau_M)*(1/(1 + exp(-(y(i,j,7)-vhalfm_M)/km_M))-y(i,j,3));  %m_M
        dy(i,j,4) = (1/mtau_A)*(1/(1 + exp(-(y(i,j,7)-vhalfm_A)/km_A))-y(i,j,4));  %m_A
        dy(i,j,5) = (1/(-20/(1+exp(-(y(i,j,7)+35)/6))+25))*(f*(1/(1 + exp((y(i,j,7)-vhalfh_A)/kh_A)))+(1-f)-y(i,j,5));  %h0_A
        dy(i,j,6) = (1/max(100,((y(i,j,7)+17)^2/4+26)))*(f*(1/(1 + exp((y(i,j,7)-vhalfh_A)/kh_A)))+(1-f)-y(i,j,6));  %h1_A
        
        if i == 1
            if j==1
                
                isection = 1/R(j)*(y(i,j+1,7)-y(i,j,7));
                
                
            elseif j==armsection+2
                
                
                isection = 1/R(j-1)*(y(i,j-1,7)-y(i,j,7));
                
            else
                
                
                isection = 1/R(j)*(y(i,j+1,7)-y(i,j,7))+1/R(j-1)*(y(i,j-1,7)-y(i,j,7));
                
            end
        elseif i == 2
            isection = 0;
        end
        
  
        if numcell == 1
            igap = 0;
        else
            if i==1
                igap = G_gap(i,j)*(y(i+1,j,7)-y(i,j,7));
            elseif i==numcell
                igap = G_gap(i,j)*(y(i-1,j,7)-y(i,j,7));
            else
                igap = G_gap(i,j)*(y(i+1,j,7)-y(i,j,7))+G_gap(i,j)*(y(i-1,j,7)-y(i,j,7));
            end
        end
        
        if i==2
            im = gbar_na(i,j)*y(i,j,1).^3.*y(i,j,2).*(y(i,j,7) - ena)+gbar_M(i,j)*y(i,j,3)*(y(i,j,7) - ek)+1/(1+exp(-(y(i,j,7)-vhalfw_A)/kw_A))*gbar_A(i,j)*y(i,j,4)*y(i,j,5)*(y(i,j,7) - ek) + (1-1/(1+exp(-(y(i,j,7)-vhalfw_A)/kw_A)))*gbar_A(i,j)*y(i,j,4)*y(i,j,6)*(y(i,j,7) - ek)+Gpas(i,j)*(y(i,j,7)-epas2);
          
        elseif i == 1
            
            im = gbar_na(i,j)*y(i,j,1).^3.*y(i,j,2).*(y(i,j,7) - ena)+gbar_M(i,j)*y(i,j,3)*(y(i,j,7) - ek)+1/(1+exp(-(y(i,j,7)-vhalfw_A)/kw_A))*gbar_A(i,j)*y(i,j,4)*y(i,j,5)*(y(i,j,7) - ek) + (1-1/(1+exp(-(y(i,j,7)-vhalfw_A)/kw_A)))*gbar_A(i,j)*y(i,j,4)*y(i,j,6)*(y(i,j,7) - ek)+Gpas(i,j)*(y(i,j,7)-epas);
          
        end
        
        dy(i,j,7) = (1/capacitance(i,j))*(-im + isection + igap +ie);

        
    end
end

dy = reshape(dy,numcell*(armsection+2)*7,1);

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xmax,imax,xmin,imin] = extrema(x)
%EXTREMA   Gets the global extrema points from a time series.
%   [XMAX,IMAX,XMIN,IMIN] = EXTREMA(X) returns the global minima and maxima 
%   points of the vector X ignoring NaN's, where
%    XMAX - maxima points in descending order
%    IMAX - indexes of the XMAX
%    XMIN - minima points in descending order
%    IMIN - indexes of the XMIN
%
%   DEFINITION (from http://en.wikipedia.org/wiki/Maxima_and_minima):
%   In mathematics, maxima and minima, also known as extrema, are points in
%   the domain of a function at which the function takes a largest value
%   (maximum) or smallest value (minimum), either within a given
%   neighbourhood (local extrema) or on the function domain in its entirety
%   (global extrema).
%
%   Example:
%      x = 2*pi*linspace(-1,1);
%      y = cos(x) - 0.5 + 0.5*rand(size(x)); y(40:45) = 1.85; y(50:53)=NaN;
%      [ymax,imax,ymin,imin] = extrema(y);
%      plot(x,y,x(imax),ymax,'g.',x(imin),ymin,'r.')
%
%   See also EXTREMA2, MAX, MIN

%   Written by
%   Lic. on Physics Carlos Adrián Vargas Aguilera
%   Physical Oceanography MS candidate
%   UNIVERSIDAD DE GUADALAJARA 
%   Mexico, 2004
%
%   nubeobscura@hotmail.com

% From       : http://www.mathworks.com/matlabcentral/fileexchange
% File ID    : 12275
% Submited at: 2006-09-14
% 2006-11-11 : English translation from spanish. 
% 2006-11-17 : Accept NaN's.
% 2007-04-09 : Change name to MAXIMA, and definition added.


xmax = [];
imax = [];
xmin = [];
imin = [];

% Vector input?
Nt = numel(x);
if Nt ~= length(x)
 error('Entry must be a vector.')
end

% NaN's:
inan = find(isnan(x));
indx = 1:Nt;
if ~isempty(inan)
 indx(inan) = [];
 x(inan) = [];
 Nt = length(x);
end

% Difference between subsequent elements:
dx = diff(x);

% Is an horizontal line?
if ~any(dx)
 return
end

% Flat peaks? Put the middle element:
a = find(dx~=0);              % Indexes where x changes
lm = find(diff(a)~=1) + 1;    % Indexes where a do not changes
d = a(lm) - a(lm-1);          % Number of elements in the flat peak
a(lm) = a(lm) - floor(d/2);   % Save middle elements
a(end+1) = Nt;

% Peaks?
xa  = x(a);             % Serie without flat peaks
b = (diff(xa) > 0);     % 1  =>  positive slopes (minima begin)  
                        % 0  =>  negative slopes (maxima begin)
xb  = diff(b);          % -1 =>  maxima indexes (but one) 
                        % +1 =>  minima indexes (but one)
imax = find(xb == -1) + 1; % maxima indexes
imin = find(xb == +1) + 1; % minima indexes
imax = a(imax);
imin = a(imin);

nmaxi = length(imax);
nmini = length(imin);                

%%%% ----- Commented Out on 052012 by Hannah ------------------%%%%%
% Maximum or minumim on a flat peak at the ends?
if (nmaxi==0) && (nmini==0)
 if x(1) > x(Nt)
  xmax = x(1);
  imax = indx(1);
  xmin = x(Nt);
  imin = indx(Nt);
 elseif x(1) < x(Nt)
  xmax = x(Nt);
  imax = indx(Nt);
  xmin = x(1);
  imin = indx(1);
 end
 return
end

%% Maximum or minumim at the ends?
% if (nmaxi==0) 
%  imax(1:2) = [1 Nt];
% elseif (nmini==0)
%  imin(1:2) = [1 Nt];
% else
%  if imax(1) < imin(1)
%   imin(2:nmini+1) = imin;
%   imin(1) = 1;
%  else
%   imax(2:nmaxi+1) = imax;
%   imax(1) = 1;
%  end
%  if imax(end) > imin(end)
%   imin(end+1) = Nt;
%  else
%   imax(end+1) = Nt;
%  end
% end
%%%%% -----------------------------------------------------------%%%%%
xmax = x(imax);
xmin = x(imin);

% NaN's:
if ~isempty(inan)
 imax = indx(imax);
 imin = indx(imin);
end

% Same size as x:
imax = reshape(imax,size(xmax));
imin = reshape(imin,size(xmin));

% Descending order:
[temp,inmax] = sort(-xmax); clear temp
xmax = xmax(inmax);
imax = imax(inmax);
[xmin,inmin] = sort(xmin);
imin = imin(inmin);


% Carlos Adrián Vargas Aguilera. nubeobscura@hotmail.com

end

Loading data, please wait...