% spinemodel_fig6.m % % Biophysical dendritic spine model which was used to generate Figure 6 of % O'Donnell et al, J Neurosci (2011). Please send any comments or questions % to cian@salk.edu. % % - Detailed biophysical model of dendritic spine as metaplastic device. % - Three functional compartments: spinehead, dendrite, connected by spine neck. % - [Ca+] tracked in spinehead. % Influx from NMDA receptors and CaV channels. % Outflux from diffusion to spine neck, internal buffering and membrane pumps % (last two lumped together). % Amplitude and timecourse of Ca transients depend on spine geometery and % synapse properties. % - Synaptic plasticity depends on [Ca]. % - includes neck capacitance, and multiple CaV channels %%%%%%%%%%%%%%%%%%%%% % SETUP %%%%%%%%%%%%%%%%%%%%% % Run parameters dt = 0.1e-3; % s duration = 1500e-3; % s; 8000s for figure starttime = cputime; % to calculate total CPU simulation time later headvolvec = 1e-18*[0.01:0.01:0.3]'; % List of potential spine head volumes headradvec = (3.*headvolvec./(4*pi)).^(1/3); % corresponding list of spine head radii % Electrical parameters Cm = 0.75; % uF/cm^2; Membrane capacitance Ra = 200E-2; % Ohm.m; Axial resistivity gbar_leak = 2; % S/m^2 = pS/um^2; Membrane leak conductance density e_leak = -70; % mV; Leak conductance reversal potential e_CaR = 10; % mV; Nernst eCa is higher (+50mV), but this approximation reduces GHK error (Koch book, 1999) Vsp0 = -50; % Initial spine head voltage Vd0 = Vsp0; % Initial dendritic voltage % Spine L = 0.5e-6; % m; neck length diam = 1*0.1e-6; % m; neck diameter Rneck = Ra*L*4/(pi*diam^2); % resulting neck resistance %headrad = 1.5299e-07; % m; Specify initial spine head radius headrad = headradvec(2); % Choose initial spine head radius from list min_headrad = 0.1e-6; % Lower limit on spine head radius Chead = 4*pi*headrad^2*1E4*Cm*1E-6; % Spine head capacitance glhead = gbar_leak*4*pi*headrad^2; % Head leak conductance alpha = 2e-10; % Ratio of headvol/gbar_ampa (assume always proportional) % Dendrite dend_area = 0.5E4; % um^2; gives input resistance of 200MOhm if gleak = 1pS/um^2 dend_diam = 2; % um dend_length = dend_area/(pi*dend_diam); % um dend_vol = pi*dend_diam^2*dend_length/4; % um^3 Cdend = dend_area*1E-8*Cm*1E-6; % dendritic capacitance Rdend = dend_area/gbar_leak; % Dendritic input resistance Istim = 200e-12; % A; stimulus current to dendrite (with default parameters Istim give V=-70, and every 100pA depolarises by 10mV. % Synapse ampa_rise = 0.18e-3; % s; Rise time for synaptic AMPAR conductance ampa_decay = 1.8e-3; % s; Decay time for synaptic AMPAR conductance nmda_rise = 2e-3; % s; Rise time for synaptic NMDAR conductance nmda_decay = 89e-3; % s; Decay time for synaptic NMDAR conductance % set ampa and nmda g separately %gbar_ampa = 0.6E-9; % S gbar_nmda = 0.090E-9; % S % set approx like Noguchi et al (2005) gbar_ampa = (4/3)*pi*headrad^3/alpha; % S %gbar_nmda = 0.035E-9 + gbar_ampa/8; % S; undercompensating; % S esyn = 0; % mV; reversal potential of synapse % Ca parameters Ca0 = 0.05E-6; % M; resting Ca concentration Casp0 = Ca0; Cad0 = Ca0; Dca = 2.2E-10; % m^2/s; Ca diffusion constant beta_sp = 1.3E-4; % /s; spine extrusion rate beta_neck = beta_sp; % /s; neck extrusion rate; assumed same as spine beta_d = 4*beta_sp; % /s; dendrite extrusion rate. 4x spine rate (to reproduce Sabatini et al, 2002) kf = 100E6; % /M/s; buffer forward (binding) rate kb = 1*5E-6*kf; % /s; buffer backward (unbinding) rate Bt_sp = 100E-6; % M; total buffer concentration in spine Bt_neck = 100E-6; % M; assumed same as spine head concentration Bt_d = 5*Bt_sp; % total buffer in dendrite (Sabatini et al, 2002) % R-type Ca current gbar_CaRspine = 75; % S/m^2 = pS/um^2 gbar_CaRdend = 10; % S/m^2 = pS/um^2 tau_mR = 0.0036; % s; activation time constant tau_hR = 0.2; % s % inactivation time constant VmR_half = -14; % mV; activation half-voltage kmR = 6.7; % 1/mV; activation slope VhR_half = -65; % mV; inactivation half-voltage khR = -11.8; % 1/mV; inactivation slope msp_R0 = 1/(1+exp(-(Vsp0-VmR_half)/kmR)); md_R0 = 1/(1+exp(-(Vd0-VmR_half)/kmR)); hsp_R0 = 1/(1+exp(-(Vsp0-VhR_half)/khR)); hd_R0 = 1/(1+exp(-(Vd0-VhR_half)/khR)); % T-type Ca current (CaV 3.1) from Traboulsie et al J Physiol 2007 gbar_CaT_spine = 75; % S/m^2 = pS/um^2 gbar_CaT_dend = 10; % S/m^2 = pS/um^2 tau_mT_0 = -0.855809*1e-3; % ms; zero activation time constant tau_mT_A = 1.493527*1e-3; % ms; prefactor on activation function tau_mT_k = 27.414182; % 1/mv; scale on activation exponential tau_hT_0 = 9.987873*1e-3; % ms % inactivation time constant tau_hT_A = 0.002883*1e-3; % ms; prefactor on inactivation function tau_hT_k = 5.598574; % 1/mv; scale on inactivation exponential VmT_half = -42.921064; % mV; activation half-voltage kmT = 5.163208; % 1/mV; activation slope VhT_half = -72.907420; % mV; inactivation half-voltage khT = -4.575763; % 1/mV; inactivation slope msp_T0 = 1 /(1+exp(-(Vsp0-VmT_half)/kmT)); md_T0 = 1 /(1+exp(-(Vd0-VmT_half)/kmT)); hsp_T0 = 1 /(1+exp(-(Vsp0-VhT_half)/khT)); hd_T0 = 1 /(1+exp(-(Vd0-VhT_half)/khT)); % L-type Ca current (CaV 1.2/1.3) from Carlin et al E J Neuro, 2000 gbar_CaL_spine = 75; % S/m^2 = pS/um^2 gbar_CaL_dend = 10; % S/m^2 = pS/um^2 tau_mL_0 = 5*1e-3; % ms; zero activation time constant tau_mL_A = 20*1e-3; % ms; prefactor on activation function tau_mLVhalf = -15; % mV; half voltage for activation tau curve tau_mL_k = 5; % 1/mv; scale on activation exponential tau_hL_0 = 20*1e-3; % ms % inactivation time constant tau_hL_A = 50*1e-3; % ms; prefactor on activation function tau_hLVhalf = -30; % mV; half voltage for inactivation tau curve tau_hL_k = 7; % mv; scale on activation exponential VmL_half = -20; % mV; activation half-voltage kmL = 6; % mV; activation slope VhL_half = -70; % mV; inactivation half-voltage khL = -6.4; % 1/mV; inactivation slope msp_L0 = 1 /(1+exp(-(Vsp0-VmL_half)/kmL)); md_L0 = 1 /(1+exp(-(Vd0-VmL_half)/kmL)); hsp_L0 = 1 /(1+exp(-(Vsp0-VhL_half)/khL)); hd_L0 = 1 /(1+exp(-(Vd0-VhL_half)/khL)); % N-type Ca current (CaV 2.2) from Huang et al, Neuroscience, 1998 gbar_CaN_spine = 75; % S/m^2 = pS/um^2 gbar_CaN_dend = 10; % S/m^2 = pS/um^2 alpha_m_N_A = 0.1*1e-3; % s; activation forward rate prefactor alpha_m_N_Vhalf = 20; % mV; activation forward rate Vhalf alpha_m_N_k = 10; % mV; activation forward rate slope beta_m_N_A = 0.4*1e-3; % s; activation backward rate prefactor beta_m_N_Vhalf = -25; % mV; activation backward rate threshold beta_m_N_k = 18; % mV; activation backward rate exponent alpha_h_N_A = 0.01*1e-3; % s; activation forward rate prefactor alpha_h_N_Vhalf = -50; % mV; activation forward rate Vhalf alpha_h_N_k = 10; % mV; activation forward rate slope beta_h_N_A = 0.1*1e-3; % s; activation backward rate prefactor beta_h_N_Vhalf = -17; % mV; activation backward rate threshold beta_h_N_k = 17; % mV; activation backward rate exponent % Initialize CaV channel state variables to steady-state values alpha_m_0_sp = (alpha_m_N_A*(Vsp0-alpha_m_N_Vhalf)/(1-exp(-(Vsp0-alpha_m_N_Vhalf)/alpha_m_N_k))); beta_m_0_sp = beta_m_N_A*exp(-(Vsp0-beta_m_N_Vhalf)/beta_m_N_k); msp_N0 = alpha_m_0_sp/(alpha_m_0_sp+beta_m_0_sp); alpha_h_0_sp = alpha_h_N_A*exp(-(Vsp0-alpha_h_N_Vhalf)/alpha_h_N_k); beta_h_0_sp = beta_h_N_A/(1+exp(-(Vsp0-beta_h_N_Vhalf)/beta_h_N_k)); hsp_N0 = alpha_h_0_sp/(alpha_h_0_sp+beta_h_0_sp); alpha_m_0_d = (alpha_m_N_A*(Vd0-alpha_m_N_Vhalf)/(1-exp(-(Vd0-alpha_m_N_Vhalf)/alpha_m_N_k))); beta_m_0_d = beta_m_N_A*exp(-(Vd0-beta_m_N_Vhalf)/beta_m_N_k); md_N0 = alpha_m_0_d/(alpha_m_0_d+beta_m_0_d); alpha_h_0_d = alpha_h_N_A*exp(-(Vd0-alpha_h_N_Vhalf)/alpha_h_N_k); beta_h_0_d = beta_h_N_A/(1+exp(-(Vd0-beta_h_N_Vhalf)/beta_h_N_k)); hd_N0 = alpha_h_0_d/(alpha_h_0_d+beta_h_0_d); % Plasticity parameters eta_ltp = 1e-9; % rate of potentiation ltp_midCa = 5.5e-6; % Potentiation curve offset ltp_slope = 0.2e-6; % Potentiation curve slope eta_ltd = eta_ltp/2; % rate of depression ltd_midCa = 4e-6; % Depression curve offset ltd_slope = 0.2e-6; % Depression curve slope %%%%%%%%%%%%%%% % STIMULATION %%%%%%%%%%%%%%% delay = 50e-3; % s; Wait time at beginning of simulation % Regular stimulation % rate = 10; % Hz % isi = 1/rate; % s % spiketimes = (delay:isi:duration)'; % s % spiketimes = [spiketimes; duration]; % spiketimes = (floor(spiketimes)/dt)*dt; % Poisson stimulation % rate = 10; % Hz % num = rate*duration/800; % isi = -(1000/rate)*log(rand(num,1)); % spiketimes = cumsum(isi); % valid = find(spiketimes