Code: Select all
TITLE G Lemon et al (2003), part I
NEURON {
POINT_PROCESS p2y2lemon1
USEION atp READ atpo VALENCE 1
USEION ip3p2y READ ip3p2yi WRITE ip3p2yi VALENCE 1
RANGE Rs, Rsp, G, IP3, PIP2, R, LR, Rp, LRp
RANGE kr, kp, ke, Rt, kdeg, ka, kd, rr, alpha, Gt, PIP2t, L, d
GLOBAL Kc, N, K1, K2
}
UNITS {
(mM) = (milli/liter)
(uM) = (micro/liter)
}
PARAMETER {
:Receptor regulation
kr = 0.000175 (/ms)
kp = 3e-5 (/ms)
ke = 6e-7 (/ms)
K1 = 1e-2 (mM)
K2 = 1 (mM)
Rt = 1.7e4 (1)
:G protein cascade
kdeg = 0.00125 (/ms)
ka = 1.7e-5 (/ms)
kd = 0.00015 (/ms)
rr = 0.0015 (/ms)
alpha = 1e-7(/ms)
Kc = 0 (mM)
Gt = 1e5 (1)
PIP2t = 2e7 (1)
d = 1.234e-2 (1)
N = 6.022e23 (1)
}
ASSIGNED {
:Receptor regulation
L (mM)
atpo(mM)
R (1)
LR (1)
Rp (1)
LRp (1)
Ri (1)
:G protein cascade
rho (1)
rh (/ms)
:Conversion
ip3p2yi (mM)
area (micron2)
diam (micron)
}
STATE {
:Receptor regulation
Rs (1)
Rsp (1)
:G protein cascade
G (1)
IP3 (1)
PIP2(1)
}
INITIAL{
:Receptor States - No ligand at start
R = Rt
LR = 0 (1)
Rs = R + LR
Rp = 0 (1)
LRp = 0 (1)
Rsp = Rp + LRp
Ri = Rt - Rs - Rsp
:G protein
G = ka*d*Gt/(ka*d + kd)
rh = (alpha*G)
PIP2 = kp*rr*PIP2t/(kp*(rh + rr) + rh*rr)
IP3 = rh*rr*PIP2t/(kp*(rh + rr) + rh*rr)
ip3p2yi = (1e18)*IP3/(N*(area*diam/4))
}
BREAKPOINT{
L = atpo
SOLVE receptor METHOD cnexp
}
DERIVATIVE receptor{
:Time evolution for R + LR = Rs and Rp + LRp = Rsp assuming fast ligand binding kinetics
:kr is recycling rate, kp is phosphorylation rate, ke is internalization rate
Rs' = kr*(Rt - Rs - Rsp) - kp*(L*Rs/(K1 + L))
Rsp' = kp*(L*Rs/(K1 + L)) - ke*(L*Rsp/(K2 + L))
:Fast Ligand binding == Steady state solution for ligand binding
R = K1*Rs/(K1 + L)
LR = L*Rs/(K1 + L)
Rp = K2*Rsp/(K2 + L)
LRp = L*Rsp/(K2 + L)
Ri = Rt - Rs - Rsp
:G protein dynamics combined in one differential equation
:Assumes the G protein based reactions are fast and well below saturation
:Also assumes G protein dissociation and association are irreversible
:G denotes Ga.GTP, and Gt is total
:rho - Bound receptor to total receptor, ka - activation rate, kd - deactivation rate
:d - Ratio of activities of unbound and bound receptor (Bound/unbound to ligand)
rho = (L*Rs)/(Rt*(K1 + L))
G' = ka*(d + rho)*(Gt - G) - kd*G
:PLC dynamics - production of ip3p2y
:PLC activated when Ga.GTP and Ca2+ binds to PLC - Assuming again rapid kinetics for PLC bindings
:Hydrolysis rate is then rh
rh = (alpha*G)
:Thus IP3 conc. changes as
IP3' = rh*PIP2 - kdeg*IP3
:Need to replenish the source of IP3 - PIP2.
:Done via a reservoir that is replenished by degraded IP3 and gives out PIP2 - both at rate rr
PIP2' = -(rh + rr)*PIP2 - rr*IP3 + rr*PIP2t :PIP2t = PIP2 + IP3 + Reservoir
ip3p2yi = (1e18)*IP3/(N*(area*diam/4))
}
COMMENT
In the above, everything other than L is in terms of molecules, not in conc. Thus they are all dimensionless and refer directly to number of molecules
IP3 is converted back to conc. at the end of DERIVATIVE block
ENDCOMMENT
Code: Select all
load_file("nrngui.hoc")
create soma
objref p2y2
p2y2 = new p2y2lemon1(0.5)
objref g
g = new Graph()
addplot(g,0)
g.addvar("soma.ip3p2yi(0.5)")
tstop = 1000
init()
run()//This time the solution becomes unstable - I just get a straight line in the plot in the 3rd Quadrant
init()
run()//Now the solution is fine - I even get the correct profiles for IP3 and G protein as discussed in Lemon et al. 2003
- I wrote a separate NMODL file with only a pair of coupled differential equations, with the method still as cnexp. I don't get this behavior.
- I changed the method to derivimplicit - now the code promptly crashes on any run, unless cvode is used. If cvode is used, the solutions are still as expected from the paper. Since the receptor has to be put in another model, which uses a lot of pointers, I am not sure if using cvode is advisable.