# XPP scripts for supplementary Figure S5A of
# Maex R, Grinevich VP, Grinevich V, Budygin E, Bencherif M, Gutkin B (2014)
# Understanding the role a7 nicotinic receptors play in dopamine efflux
# in nucleus accumbens. ACS Chemical Neuroscience 5, 1032-1040.
#
# Here I want to implement the facilitation of the receptor
# at very low (partial) agonist concentration, as described
# by Smulders et al. (2005) and Prickaerts et al. (2012).
# The graph plots the fraction of open a7 channels across time.
# Naming conventions
#
# V_ = membrane voltage of (neuron or neuron population)
# R_ = release of (transmitter)
# C_ = concentration of (transmitter)
# I_ = membrane current of (channel or receptor)
# P_ = presynatic membrane current of (receptor)
# tau_ = time-constant of
# tmin_ = minimum value of time-constant of
# tmax_ = maximum value of time-constant of
# inf_ = steady-state value of
# act_ = level of activation of (receptor)
# des_ = level of desensitisation of (receptor)
# w_ = weight of
# EC50_ = halfmax concentration of steady-state response
# Ktau_ = halfmax concentration of time-constant
# stim_ = stimulation by
# glu = glutamate or glu-ergic neuron (population)
# dop = dopamine or dopamine-ergic neuron (population)
# gab = gaba or gaba-ergic neuron (population)
# ach = acetylcholine
# nic = nicotine
# a7 = alpha7-type nicACh receptor
# a4b2 = alpha4-beta2-type nicACh receptor
# bas = basal
#
# Units
# time seconds
# concentration micromolar
# Transfer functions (input-output)
# sigmoid
f(x) = 1/(1 + exp(-x))
# halfwave rectification
hwr(x) = heav(x) * x
# clipped at 0 and 1
clip(x) = (1 + (x-1) * (heav (1-x))) * heav(x)
# Hill equation
# Hill(x,K,n) = 1 / (1 + (K/x)^n)
Hill(x,K,n) = x^n / (x^n + K^n)
# competitive Hill
compHill(x,y,z,wx,wy,wz,Kx,Ky,Kz,n) = \
(wx * x^n/(x^n + Kx^n * (1 + (y/Ky)^n + (z/Kz)^n))) + \
(wy * y^n/(y^n + Ky^n * (1 + (x/Kx)^n + (z/Kz)^n))) + \
(wz * z^n/(z^n + Kz^n * (1 + (x/Kx)^n + (y/Ky)^n)))
# steady-state Hill function with cross-terms representing the effect of
# the simultaneous binding of two different compounds, for instance ACh and EVP
Hillcross(x,y,z,wx,wy,wz,Kx,Ky,Kz,n,fxy,fyz,fxz) = \
((wx * (x/Kx)^n) + \
(wy * (y/Ky)^n) + \
(wz * (z/Kz)^n) + \
(fxy * (x/Kx) * (y/Ky)) + \
(fxz * (x/Kx) * (z/Kz)) + \
(fyz * (y/Kz) * (z/Kz))) / (1 + (z/Kz)^n + (x/Kx)^n + (y/Ky)^n + \
2 * (z/Kz) + 2 * (x/Kx) + 2 * (y/Ky) + \
2 * (x/Kx) * (y/Ky) + \
2 * (x/Kx) * (z/Kz) + \
2 * (y/Ky) * (z/Kz) )
# rectangular pulse
rect_pulse(t,from,to,amp) = amp * (heav (t-from) - heav (t-to))
# THE MODEL
# We now have two neuron populations sharing inputs from
# glu neurons (mostly a7-driven) and through a4b2 receptors.
# The parameter coefficients r and s specify the balance
# between these inputs.
par s=0.8,r=0.8
# the dynamics of the dopaminergic neuron population V_dop
V_dop' = (-V_dop - I_gab + hwr (I_basDop + s*I_glu + r*I_b2)) / tau_Vdop
par tau_Vdop=0.02,I_basDop=0.1
aux aux_hwr=hwr (I_basDop + s*I_glu + r*I_b2)
# the dynamics of the gaba-ergic neuron population V_gab
V_gab' = (-V_gab + hwr ((1-s)*I_glu + (1-r)*I_b2)) / tau_Vgab
par tau_Vgab=0.02
# the (stationary) glutamatergic input I_glu to dopamine (and GABA)neurons
I_glu = w_glu * clip (V_glu + P_glu)
par w_glu=1,V_glu=0.1
# the (varying) gaba-ergic input I_gab to dopamine neurons
I_gab = w_gab * clip (V_gab + I_basGab)
aux auxI_gab=I_gab
par w_gab=1.5
par I_basGab=-0.05
aux auxI_gab=I_gab
# the presynaptic facilitation of glu input to dopamine (and GABA) neurons
P_glu = P_a7
# the dynamics of presynaptic a7 receptors (same single variable for
# glu inputs to dopamine cells in VTA and medium spiny neurons in striatum)
P_a7 = act_a7 * (1 - des_a7)
aux auxP_a7=P_a7
aux auxact_a7=act_a7
aux auxsen_a7=1-des_a7
# /* activation of alpha7 */
act_a7' = (- act_a7 + inf_actA7) / tau_actA7
inf_actA7 = Hillcross(C_ach + C_chol, C_nic, C_agA7, \
1.0, w_actNicA7, w_actAgA7, \
EC50_A7ach, EC50_A7nic, EC50_A7ag, Hill_actA7, \
0, 0, crossactA7)
par crossactA7=0.3
aux ainf_actA7=inf_actA7
par EC50_A7=80
par EC50_A7ach=68
par EC50_A7nic=13
par Hill_actA7=1.73
par EC50_A7ag=0.03
par w_actNicA7=0.8
par w_actAgA7=0
par tau_actA7=0.005
# /* desensitisation of alpha7 */
des_a7' = (- des_a7 + inf_desA7) / tau_desA7
inf_desA7 = compHill( 0, C_nic, C_agA7, \
1, 1, 0, \
1, IC50_A7nic, IC50_A7ag, Hill_desA7)
par crosssenA7=0
aux ainf_desA7=inf_desA7
par IC50_A7ach=1.3
par IC50_A7nic=1.3
par IC50_A7ag=0.002
par Hill_desA7=2
tau_desA7 = tmin_desA7 + \
tmax_desA7 * (1 - inf_desA7)
aux aux_tdA7=tau_desA7
par tmin_desA7=0.05
par tmax_desA7=120
# the dynamics of somatic a4b2 receptors on the soma/dendrite of dopamine neurons
I_b2 = act_b2 * (1 - des_b2)
aux auxI_b2=I_b2
aux auxact_b2=act_b2
aux auxsen_b2=1-des_b2
# /* activation of alpha4beta2 */
act_b2' = (- act_b2 + inf_actB2) / tau_actB2
inf_actB2 = compHill(C_ach, C_nic, C_agB2, \
1, w_actNicB2, w_actAgB2, \
EC50_B2ach, EC50_B2nic, EC50_B2ag, Hill_actB2)
par EC50_actB2=30
par EC50_B2ach=30
par EC50_B2ag=30
par EC50_B2nic=0.23
par Hill_actB2=1.05
par w_actNicB2=2
par w_actAgB2=1
par tau_actB2=0.005
# /* desensitisation */
des_b2' = (- des_b2 + inf_desB2) / tau_desB2
inf_desB2 = Hill (C_nic + C_agB2, EC50_desB2, Hill_desB2)
tau_desB2 = tmin_desB2 + \
tmax_desB2 * Hill (Ktau_desB2, C_nic + C_agB2, Htau_desB2)
aux aux_tdB2=tau_desB2
par EC50_desB2=0.061
par Hill_desB2=0.5
par tmin_desB2=0.5
par tmax_desB2=600
par Ktau_desB2=0.11
par Htau_desB2=3
# the dynamics of dopamine release AND RE-UPTAKE
R_dop' = (C_basDop * (1 + (V_dop-ss_Vdop) / ss_Vdop) - R_dop) / tau_Rdop - R_dop * Vmax_Rdop / (R_dop + EC50_Rdop)
par tau_Rdop=0.2,P_a7bas=0.01
# the assumed basal DA production, denotated C_basDop, in the absence of re-uptake
par C_basDop=0.1
# the following values are from Chen and Budygin 2007
par Vmax_Rdop=1.3
par EC50_Rdop=0.2
# the physiological (ach) and pharmacological (nic) stimuli
# /* acetylcoline HERE I USE CHOLINE BECAUSE THIS SELECTIVELY ACTIVATES THE A7s */
C_ach' = (-C_ach + stim_ach) / 1
stim_chol = bas_chol + rect_pulse (t, 1, 10, dose_chol) + rect_pulse (t, 31, 35, dose_chol) \
+ rect_pulse (t, 61, 65, dose_chol) \
+ rect_pulse (t, 91, 95, dose_chol) \
+ rect_pulse (t, 121, 125, dose_chol) \
+ rect_pulse (t, 151, 155, dose_chol) \
+ rect_pulse (t, 181, 185, dose_chol) \
+ rect_pulse (t, 211, 215, dose_chol) \
+ rect_pulse (t, 241, 245, dose_chol) \
+ rect_pulse (t, 271, 275, dose_chol) \
+ rect_pulse (t, 301, 305, dose_chol) \
+ rect_pulse (t, 331, 335, dose_chol) \
+ rect_pulse (t, 361, 365, dose_chol) \
+ rect_pulse (t, 391, 395, dose_chol) \
+ rect_pulse (t, 421, 425, dose_chol) \
+ rect_pulse (t, 481, 485, dose_chol) \
+ rect_pulse (t, 541, 545, dose_chol) \
+ rect_pulse (t, 571, 575, dose_chol) \
+ rect_pulse (t, 601, 605, dose_chol) \
+ rect_pulse (t, 631, 635, dose_chol)
C_chol' = (-C_chol + stim_chol) / 1
par bas_ach=0
par bas_chol=0
par dose_ach=0
par stim_ach=0
# FOR THE UPPER PLOT PULSES OF 1 microM ARE APPLIED
par dose_chol=1
# FOR THE LOWER PLOT PULSES OF 60 microM ARE APPLIED
# par dose_chol=60
# /* nicotine */
stim_nic = bas_nic + rect_pulse (t, 66, 67, dose_Nic)
par dose_Nic=0
# we now simulate an alpha function as a second-order ode
# the alpha-function is area normalized
# to calculate its peak, take amplitude and divide by 10*e
C_nic_nic' = (-C_nic_nic + stim_nic)/tau1_nic
C_nic' = (-C_nic + C_nic_nic)/tau2_nic
par tau1_nic=10,tau2_nic=10,bas_nic=0
# /* alpha7 agonist */
stim_a7 = rect_pulse (t, 170, 510, dose_Aga7)
par dose_Aga7=0.003
C_agA7A7' = (-C_agA7A7 + stim_a7) / tau1_agA7
C_agA7' = (-C_agA7 + C_agA7A7) / tau2_agA7
par tau1_agA7=10,tau2_agA7=10
# /* alpha4 beta2 agonist */
stim_b2 = rect_pulse (t, 60, 61, 0)
C_agB2' = (-C_agB2 + stim_b2) / tau_agB2
par tau_agB2=10
# /* calculating the steady-states (in reverse order) */
ss_Ib2 = compHill( bas_ach, 0, 0, \
1, w_actNicB2, w_actAgB2, \
EC50_B2ach, EC50_B2nic, EC50_B2ag, Hill_actB2)
ss_Pa7 = compHill( bas_ach + bas_chol, 0, 0, \
1, w_actNicA7, w_actAgA7, \
EC50_A7ach, EC50_A7nic, EC50_A7ag, Hill_actA7)
ss_Pglu = ss_Pa7
ss_Iglu = w_glu * clip (V_glu + ss_Pglu)
ss_Vgab = hwr ((1-s)*ss_Iglu + (1-r)*ss_Ib2)
ss_Igab = w_gab * clip (ss_Vgab + I_basGab)
ss_Vdop = - ss_Igab + hwr (I_basDop + s*ss_Iglu + r*ss_Ib2)
# aux auxss_Vdop=ss_Vdop
# the initial conditions
act_a7(0)=0
des_a7(0)=0
act_b2(0)=0
des_b2(0)=0
V_dop(0)=0
V_gab(0)=0
C_nic(0)=bas_nic
C_nic_nic(0)=bas_nic
# the numerical parameters
@ MAXSTOR=500000
@ BOUNDS=1e+8
@ DT=0.005
@ TOTAL=700
@ xlo=0
@ xhi=700
@ ylo=0
@ yhi=0.001
@ METHOD=rk
@ NJP=200
@ YP=auxP_a7
done