Dopamine-modulated medium spiny neuron, reduced model (Humphries et al. 2009)

 Download zip file 
Help downloading and running models
Accession:128818
We extended Izhikevich's reduced model of the striatal medium spiny neuron (MSN) to account for dopaminergic modulation of its intrinsic ion channels and synaptic inputs. We tuned our D1 and D2 receptor MSN models using data from a recent (Moyer et al, 2007) large-scale compartmental model. Our new models capture the input-output relationships for both current injection and spiking input with remarkable accuracy, despite the order of magnitude decrease in system size. They also capture the paired pulse facilitation shown by MSNs. Our dopamine models predict that synaptic effects dominate intrinsic effects for all levels of D1 and D2 receptor activation. Our analytical work on these models predicts that the MSN is never bistable. Nonetheless, these MSN models can produce a spontaneously bimodal membrane potential similar to that recently observed in vitro following application of NMDA agonists. We demonstrate that this bimodality is created by modelling the agonist effects as slow, irregular and massive jumps in NMDA conductance and, rather than a form of bistability, is due to the voltage-dependent blockade of NMDA receptors
Reference:
1 . Humphries MD, Lepora N, Wood R, Gurney K (2009) Capturing dopaminergic modulation and bimodal membrane behaviour of striatal medium spiny neurons in accurate, reduced models. Front Comput Neurosci 3:26 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Neostriatum medium spiny direct pathway GABA cell;
Channel(s):
Gap Junctions:
Receptor(s): D1; D2; GabaA; AMPA; NMDA;
Gene(s):
Transmitter(s): Dopamine;
Simulation Environment: MATLAB;
Model Concept(s): Action Potential Initiation; Parameter Fitting; Simplified Models; Parkinson's; Bifurcation;
Implementer(s): Humphries, Mark D [m.d.humphries at shef.ac.uk];
Search NeuronDB for information about:  Neostriatum medium spiny direct pathway GABA cell; D1; D2; GabaA; AMPA; NMDA; Dopamine;
function [best_coeffs,best_fun,AICs,SS,best_residuals] = fitcurves(xdata,ydata,fits,LBC,UBC,IVC,options,varargin)

% FITCURVES test fit of curves to data
%   [C,F,AS,SS,R] = FITCURVES(X,Y,FITS,LB,UB,IV,OPTS) fits the set of functions requested
%   in FITS to find the best function that describes Y = F(X). The best function
%   is the one with the lowest AIC score.
%
%   [...] = FITCURVES(...,'S') will return the best function as the lowest
%   sum-of-squares, ignoring the number of free parameters - useful when
%   all you want is the very best fit!
%   
%   The set of functions is defined by the array FITS, whose numeric entries 
%   request the following:
%       1 = linear      (2 par)         a + b*x
%       2 = exponential (2 par)         a*exp(b*x)  
%       3 = exponential (3 par)         a + b*exp(c*x)
%       4 = double exponential (4 par)  a*exp(b*x) + c*exp(d*x)
%       5 = double exponential (5 par)  a + b*exp(c*x) + d*exp(e*x)
%       6 = triple exponential (6 par)  a*exp(b*x) + c*exp(d*x) + e*exp(f*x)
%       7 = power law  (2 par)          a*x^b
%       8 = power law (3 par)           a + b*x^c
%       9 = inverse second order (3 par)a + b/x + c/x^2
%      10 = Rayleigh (2 par)            a*x*exp(-b*sqrt(x))
%      11 = Gamma function (3 par)  1/(b^a * Gamma(a)) * x^(a-1) * exp(x/b) * c [x>0; set c=1 to get gamma PDF]    
%      12 = piece-wise linear#1 (5 par) 1: a+ b*x (x <= c) 2: d + e*x (x >= c)
%      13 = truncated power law (3 par) a*x^b * exp(x*c)    
%      14 = Pareto function (1 par)     1 - 1/(x^a)     
%      15 = modded power law I (2 par)  a*(1-x)^-b)
%      16 = modded power law II (2 par) a*(1+x)^b
%      17 = piece-wise power law (5 par) 1: a*x^b (x <= c) 2: d*x^e (x > c)
%      18 = inverse cubic polynomical (4 par) a/x^3 + b/x^2 + c/x + d  
%      19 = quadratic (3 par)           a + bx + cx^2 
%      20 = exponential rise-to-max (2par) a*(1-exp(-b*x))
%      21 = exponential rise-to-max (3par) a + b*(1-exp(-c*x))
%      22 = dbl exponential max (4par)  a*(1-exp(-b*x)) + c*(1-exp(-d*x))
%      23 = inverse first-order (2 par) a + b/x 
%      24 = stretched exponential decay (3 par)  a*exp(-x/b)^c
%      25 = cubic (4par):               a + bx + cx^2 + dx^3 
%      26 = quartic (5 par):            a + bx + cx^2 + dx^3 + e*x^4 
%      27 = hyperbola, double (4par):   a*x/(b+x) + c*x/(d+x)
%      28 = hyperbola, double (5par):   a*x/(b+x) + c*x/(d+x) + e*x
%      29 = quintic (6 par):            a + bx + cx^2 + dx^3 + e*x^4 + f*x^5
%      30 = exp+linear (4 par):         a + b*exp(c*x) + d*x
%      31 = rational (4 par):           (a + bx)/(1 + cx + dx^2)
%      32 = rational (5 par):           (a + bx + cx^2)/(1 + dx +ex^2)
%      33 = Gaussian (3 par):           a*exp(-(x-b)^2/2c^2) 
%      34 = double Gaussian (6 par):    a*exp(-(x-b)^2/2c^2) + d*exp(-(x-e)^2/2f^2)
%      35 = log Gaussian (3 par):       a*exp(-(ln(x)-b)^2/2c^2)     
%   
%   The range of values over which the coefficients for the requested functions 
%   are tested can be limited by lower (LB) and upper (UB) bound values.
%   Similarly, the initial values for the search can be given in the cell
%   array IV. The cell arrays LB, UB, and IV must have the same number of elements as FITS.
%   In each cell, there must be as many values as there are coefficients.
%   Any entry can be omitted by putting [] (including the entire array). OPTS is the
%   standard MatLab OPTIONS structure; set [] to omit (e.g. if using flag 'S').
%
%   Returns: C a cell array of the coefficients of the best-fit function; F
%   the inline function that was the best-fit; AS the array of AIC scores for every tested function; SS the
%   sum of squares all fitted functions; R the residuals for the *best*
%   fit function.
%
%   Note#1: TO FIX: random allocation of initial values always in [0 1] if
%   not specified; should be within [LB, UB] if these specified!!
%
%   Note#2: To get e.g. exponential decay, set bounds of LB = {[-Inf -Inf]}, UB =
%   {[Inf 0]}; 
%
%   Mark Humphries 22/6/2009

if any(fits > 35)
    error('Non-existent curve chosen')
end

rtn = 'AIC';
if nargin > 7 & strfind(varargin{1},'S')
    rtn = 'SS';
end

% linear fit
linfun = inline('x(1) + x(2) .* xdata','x','xdata');
% exponential fit (2 par)
expAfun = inline('x(1) .* exp(x(2).*xdata)','x','xdata');
% exponential fit (3 par)
expBfun = inline('x(1) + x(2) .* exp(x(3).*xdata)','x','xdata');
% double exponential fit (4 par)
dblexpAfun = inline('x(1) .* exp(x(2).*xdata) + x(3) .* exp(x(4).*xdata)','x','xdata');
% double exponential fit (5 par)
dblexpBfun = inline('x(1) + x(2) .* exp(x(3).*xdata) + x(4) .* exp(x(5).*xdata)','x','xdata');
% triple exponential fit
triexpfun = inline('x(1) .* exp(x(2).*xdata) + x(3) .* exp(x(4).*xdata) + x(5) .* exp(x(6).*xdata)','x','xdata');
% power law (2 par)
pwrAfun = inline('x(1) .* xdata.^x(2)','x','xdata');
% power law (3 par)
pwrBfun = inline('x(1) + x(2) .* xdata.^x(3)','x','xdata');
% inverse 2nd order fit
inv2fun = inline('x(1) + x(2)./xdata + x(3)./xdata.^2','x','xdata');
% Rayleigh function
ralfun = inline('x(1) .*xdata .* exp(-x(2).* sqrt(xdata))','x','xdata');
% Gamma function
gamdistfun = inline('1/(x(2).^x(1) .* gamma(x(1))) .* xdata.^(x(1)-1) .* exp(-xdata ./ x(2)) .* x(3)','x','xdata');
% piecewise linear fit 
piecelinfun = inline('(x(1) + x(2) .* xdata) .* (xdata<=x(3)) + (x(4) + x(5) .* xdata) .* (xdata > x(3)) ','x','xdata');
% truncated power law
trunpwrfun = inline('x(1) .* xdata.^x(2) .* exp(x(3) .* xdata)','x','xdata');
% Pareto 
Paretofun = inline('1 - 1 ./ xdata.^x(1)','x','xdata');
% modified power law I
modpwrAfun = inline('x(1) .* (1 - xdata.^-x(2))','x','xdata');
% modified power law II
modpwrBfun = inline('x(1) .* (1 + xdata.^x(2))','x','xdata');
% piecewise power law
piecepwrfun = inline('(x(1) .* xdata .^ x(2)) .* (xdata<=x(3)) + (x(4) .* xdata .^ x(5)) .* (xdata > x(3)) ','x','xdata');
% inverse cubic
invcubicfun = inline('x(1) ./ xdata.^3 + x(2) ./ xdata.^2 + x(3) ./ xdata + x(4)','x','xdata');
% quadratic
quadfun = inline('x(1) + x(2) .* xdata + x(3).*xdata.^2','x','xdata');
% exponential rise-to-max (2par) a*(1-exp(-b*x))
expmaxAfun = inline('x(1) .* (1-exp(-x(2).*xdata))','x','xdata');
% exponential rise-to-max (3par) a + b*(1-exp(-c*x))
expmaxBfun = inline('x(1) + x(2).* (1-exp(-x(3).*xdata))','x','xdata');
% dbl exponential max (4par)  a*(1-exp(-b*x)) + c*(1-exp(-d*x))
dblexpmaxAfun = inline('x(1) .* (1-exp(-x(2).*xdata))+x(3) .* (1-exp(-x(4).*xdata))','x','xdata');
% inverse first-order
invfun = inline('x(1) + x(2)./xdata','x','xdata');
% stretched exponential decay (3 par)
strexpfun = inline('x(1).*exp(-(xdata./x(2)).^x(3))','x','xdata');
% cubic
cubicfun =  inline('x(1) + x(2) .* xdata + x(3).*xdata.^2 + x(4).*xdata.^3','x','xdata');
% quartic 
quarticfun =  inline('x(1) + x(2) .* xdata + x(3).*xdata.^2 + x(4).*xdata.^3 + x(5).*xdata.^4','x','xdata');
% hyperbola, double (4 par)
hyprdblAfun =  inline('x(1)*xdata ./ (x(2) + xdata) + x(3).*xdata ./(x(4)+ xdata)','x','xdata');
% hyperbola, double (5 par)
hyprdblBfun =  inline('x(1)*xdata ./ (x(2) + xdata) + x(3).*xdata ./(x(4)+ xdata) + x(5) .* xdata','x','xdata');
% quintic
quinticfun =  inline('x(1) + x(2) .* xdata + x(3).*xdata.^2 + x(4).*xdata.^3 + x(5).*xdata.^4 + x(6).*xdata.^5','x','xdata');
% exp and linear (4 par)
explinfun = inline('x(1) + x(2) .* exp(x(3).*xdata) + x(4).*xdata','x','xdata');
% rational (4 par)
rtnl4fun = inline('(x(1) + x(2) .* xdata) ./ (1 + x(3).*xdata + x(4).*xdata.^2)','x','xdata');
% rational (5 par)
rtnl5fun = inline('(x(1) + x(2) .* xdata + x(3).*xdata.^2) ./ (1 + x(4).*xdata + x(5).*xdata.^2)','x','xdata');
% Gaussian (3 par)
Gaussfun = inline('x(1).*exp(-(xdata-x(2)).^2./x(3).^2)','x','xdata');
% double Gaussian (6 par)
dblGaussfun = inline('x(1).*exp(-(xdata-x(2)).^2./x(3).^2) + x(4).*exp(-(xdata-x(5)).^2./x(6).^2)','x','xdata');
% Gaussian (3 par)
logGaussfun = inline('x(1).*exp(-(log(xdata)-x(2)).^2./x(3).^2)','x','xdata');


nfits = length(fits);
SS = zeros(nfits,1);
coeffs = cell(nfits,1);
residuals = cell(nfits,1);

if isempty(LBC) LBC = cell(nfits,1); end
if isempty(UBC) UBC = cell(nfits,1); end
if isempty(IVC) IVC = cell(nfits,1); end

for loop = 1:nfits
    if ~isempty(LBC{loop})
        LB = LBC{loop};
        UB = UBC{loop};
    else
        LB = [];
        UB = [];
    end
    if ~isempty(IVC{loop}) IV = IVC{loop}; end

    switch fits(loop)
        case 1
            if isempty(IVC{loop}) IV = rand(2,1); end
            
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(linfun,IV, xdata, ydata,LB,UB,options);   
        case 2
            if isempty(IVC{loop}) IV = rand(2,1); end

            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(expAfun,IV, xdata, ydata,LB,UB,options);   
        case 3
            if isempty(IVC{loop}) IV = rand(3,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(expBfun,IV, xdata, ydata,LB,UB,options);

        case 4 
            if isempty(IVC{loop}) IV = rand(4,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(dblexpAfun,IV, xdata, ydata,LB,UB,options);   

        case 5
            if isempty(IVC{loop}) IV = rand(5,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(dblexpBfun,IV, xdata, ydata,LB,UB,options);   
        case 6
            if isempty(IVC{loop}) IV = rand(6,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(triexpfun,IV, xdata, ydata,LB,UB,options);   

        case 7
            if isempty(IVC{loop}) IV = rand(2,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(pwrAfun,IV, xdata, ydata,LB,UB,options);   

        case 8
            if isempty(IVC{loop}) IV = rand(3,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(pwrBfun,IV, xdata, ydata,LB,UB,options);   
        case 9
            if isempty(IVC{loop}) IV = rand(3,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(inv2fun,IV, xdata, ydata,LB,UB,options);   
        case 10
            if isempty(IVC{loop}) IV = rand(2,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(ralfun,IV, xdata, ydata,LB,UB,options);    
        case 11            
            if isempty(IVC{loop}) IV = rand(3,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(gamdistfun,IV, xdata, ydata,LB,UB,options);          
        case 12
            if isempty(IVC{loop}) IV = rand(5,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(piecelinfun,IV, xdata, ydata,LB,UB,options);
%             max_ix = find(ydata == max(ydata));
%             x1 = xdata(1:max_ix); y1 = ydata(1:max_ix);
%             x2 = xdata(max_ix:end); y2 = ydata(max_ix:end);
%             [coeffs1,SS1,residuals1] = lsqcurvefit(linfun,IV(1:2), x1, y1,LB,UB);          
%             [coeffs2,SS2,residuals2] = lsqcurvefit(linfun,IV(3:4), x2, y2,LB,UB);          
%             coeffs{loop} = [coeffs1 coeffs2];
%             SS(loop) = SS1 + SS2;
%             residuals{loop} = [residuals1; residuals2];
        case 13
            if isempty(IVC{loop}) IV = rand(3,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(trunpwrfun,IV, xdata, ydata,LB,UB,options);              
        case 14
            if isempty(IVC{loop}) IV = rand(1,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(Paretofun,IV, xdata, ydata,LB,UB,options);              
        case 15
            if isempty(IVC{loop}) IV = rand(2,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(modpwrAfun,IV, xdata, ydata,LB,UB,options);              
        case 16
            if isempty(IVC{loop}) IV = rand(2,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(modpwrBfun,IV, xdata, ydata,LB,UB,options);              
        case 17
            if isempty(IVC{loop}) IV = rand(5,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(piecepwrfun,IV, xdata, ydata,LB,UB,options);              
        case 18
            if isempty(IVC{loop}) IV = rand(4,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(invcubicfun,IV, xdata, ydata,LB,UB,options);              
        case 19
            if isempty(IVC{loop}) IV = rand(3,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(quadfun,IV, xdata, ydata,LB,UB,options);              
        case 20
            if isempty(IVC{loop}) IV = rand(2,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(expmaxAfun,IV, xdata, ydata,LB,UB,options);              
        case 21
            if isempty(IVC{loop}) IV = rand(3,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(expmaxBfun,IV, xdata, ydata,LB,UB,options);              
        case 22
            if isempty(IVC{loop}) IV = rand(4,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(dblexpmaxAfun,IV, xdata, ydata,LB,UB,options);              
        case 23
            if isempty(IVC{loop}) IV = rand(2,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(invfun,IV, xdata, ydata,LB,UB,options);                  
        case 24
            if isempty(IVC{loop}) IV = rand(3,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(strexpfun,IV, xdata, ydata,LB,UB,options);                      
        case 25
            if isempty(IVC{loop}) IV = rand(4,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(cubicfun,IV, xdata, ydata,LB,UB,options);                      
        case 26
            if isempty(IVC{loop}) IV = rand(5,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(quarticfun,IV, xdata, ydata,LB,UB,options);                      
        case 27
            if isempty(IVC{loop}) IV = rand(4,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(hyprdblAfun,IV, xdata, ydata,LB,UB,options);                      
        case 28
            if isempty(IVC{loop}) IV = rand(5,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(hyprdblBfun,IV, xdata, ydata,LB,UB,options);                      
        case 29
            if isempty(IVC{loop}) IV = rand(6,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(quinticfun,IV, xdata, ydata,LB,UB,options);                      
        case 30
            if isempty(IVC{loop}) IV = rand(4,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(explinfun,IV, xdata, ydata,LB,UB,options);                      
        case 31
            if isempty(IVC{loop}) IV = rand(4,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(rtnl4fun,IV, xdata, ydata,LB,UB,options);                      
        case 32
            if isempty(IVC{loop}) IV = rand(5,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(rtnl5fun,IV, xdata, ydata,LB,UB,options);                      
        case 33
            if isempty(IVC{loop}) IV = rand(3,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(Gaussfun,IV, xdata, ydata,LB,UB,options);                      
        case 34
            if isempty(IVC{loop}) IV = rand(6,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(dblGaussfun,IV, xdata, ydata,LB,UB,options);                      
        case 35
            if isempty(IVC{loop}) IV = rand(3,1); end
            [coeffs{loop},SS(loop),residuals{loop}] = lsqcurvefit(logGaussfun,IV, xdata, ydata,LB,UB,options);                      
   
    end
end

AICs = zeros(1,nfits);
for j = 1:nfits
    AICs(j) = AIC(SS(j),length(ydata),length(coeffs{j}));
end

switch rtn 
    case {'AIC'}
    %%% then find smallest AIC
    [AICsort,I] = sort(AICs);
    case {'SS'}
    %% then find smallest sum-of-squares
    [SSsort,I] = sort(SS);
    otherwise
        error('returned function criterion unknown')
end

best_model = fits(I(1));        
best_coeffs = coeffs{I(1)};
best_SS = SS(I(1));
best_residuals = residuals{I(1)};

switch best_model
    case 1
        best_fun = linfun;
    case 2
        best_fun = expAfun;
    case 3
        best_fun = expBfun;
    case 4
        best_fun = dblexpAfun;
    case 5
        best_fun = dblexpBfun;      
    case 6
        best_fun = triexpfun;
    case 7
        best_fun = pwrAfun;
    case 8
        best_fun = pwrBfun;
    case 9
        best_fun = inv2fun;
    case 10
        best_fun = ralfun;
    case 11
        best_fun = gamdistfun;    
    case 12
        best_fun = piecelinfun;    
    case 13
        best_fun = trunpwrfun;
    case 14
        best_fun = Paretofun;
    case 15
        best_fun = modpwrAfun;
    case 16
        best_fun = modpwrBfun;
    case 17
        best_fun = piecepwrfun;   
     case 18
        best_fun = invcubicfun;   
    case 19
        best_fun = quadfun;
    case 20
        best_fun = expmaxAfun;
    case 21
        best_fun = expmaxBfun;
    case 22
        best_fun = dblexpmaxAfun;       
    case 23
        best_fun = invfun;
    case 24
        best_fun = strexpfun;
    case 25
        best_fun = cubicfun;      
    case 26
        best_fun = quarticfun;   
    case 27
        best_fun = hyprdblAfun;   
    case 28
        best_fun = hyprdblBfun;      
    case 29
        best_fun = quinticfun;    
    case 30
        best_fun = explinfun;           
    case 31
        best_fun = rtnl4fun;
    case 32
        best_fun = rtnl5fun;       
    case 33
        best_fun = Gaussfun;
    case 34
        best_fun = dblGaussfun;
    case 35
        best_fun = logGaussfun;
                
end


Loading data, please wait...