Inconsistent behavior on modeling p2y2 Receptor
Posted: Mon Nov 06, 2017 12:03 am
I am trying to model th p2y2 receptor, as detailed in Lemon et al. 2003 paper (only the IP3 generation), and have written the following NMODL file for the same: (Modified version to remove IP3 dependence on Ca2+ - just putting Kc = 0 in the paper's model)
The model as such works fine, but the trouble that I am getting is that when it is run the first time (either through hoc interpreter or by loading a hoc file directly) the solution comes out to be unstable. Running it again just after that (without closing neuron) gives the expected solution. The hoc code that I am using is the following:
I am not sure what causes this behavior - I suspected the coupled differential equation in the NMODL, so I tried two ways to see if that is the case:
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.