Inconsistent behavior on modeling p2y2 Receptor

NMODL and the Channel Builder.
Post Reply
Archit
Posts: 2
Joined: Sun Jun 19, 2016 2:19 pm

Inconsistent behavior on modeling p2y2 Receptor

Post by Archit » 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)

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
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:

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 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:
  • 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.
I would be grateful for any insight as to why this is happening. (Also, I was not sure if this was the right place to put this question - if its not, please feel free to move it to the correct section, or point me to do so - this is my first post on this forum).

ted
Site Admin
Posts: 5115
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Inconsistent behavior on modeling p2y2 Receptor

Post by ted » Mon Nov 06, 2017 5:41 pm

I suspect the problem is in the treatment of R and LR. They are initialized to values that are far from their quasi steady state, and in the first time step of the solution R collapses to a value two orders of magnitude smaller than its initial value, while LR jumps from 0 to nearly 2e4. Thereafter they both decrease very slowly. That would be OK except for one thing: this transition is governed by an algebraic equation, not an ODE. Stability can be a problem for systems in which the dynamics of one or more states are approximated by an algebraic equation, as is done in your implementation. Suggest you replace the algebraic formulation with an ODE, even one in which the time constant (or transition rate constant) is completely arbitrary (you'll probably find that a wide range of time constants produces simulation results that are essentially identical; it's not the exact value that matters, it's the fact that it's much faster than the time constants that govern the other states, and the fact that it exists). This won't be a serious problem for performance, because this model already has 5 other states governed by ODEs.

Now, why does this instability happen only on the first run? Probably "some variable that matters" is not initialized to the same value on each run; if it started at the same value on each run, then each run would produce identical results (in this case, it would blow up). Suggest you add a printf statement to the INITIAL block that reports the value of every variable that appears on the left hand side of an assignment statement (NMODL's printf statement syntax is a subset of C's), then see if that reveals a difference between the values reported before the first and second runs. This is potentially important, because you really don't want the results of a previous run to affect the next run--bad for reproducibility and debugging.

Archit
Posts: 2
Joined: Sun Jun 19, 2016 2:19 pm

Re: Inconsistent behavior on modeling p2y2 Receptor

Post by Archit » Mon Nov 13, 2017 2:32 pm

Thank you Ted for your quick and insightful response, and really sorry for my late (and somewhat long) one.

As you suggested, I tried converting the algebraic approximations to calculate R and LR to ODE - this does not solve the problem. I even removed the algebraic equations entirely - note that the mechanism per se does not depend on the values of L, LR, Rp, LRp, Ri - and found that the instable behavior still remains.

I then tried the second suggestion in your post - to print out the initial values of the variables by putting a printf statement in NMODL. Doing so reveals that indeed there is a variable that is not set at start of the first run, but is set after the first run. This variable is L (the simulations were done at constant atpo for now). Putting the statement L = atpo in the INITIAL block removes the instability, irrespective of whether the algebraic approximations are used or not. The algebraic approx. also behave as expected, including the jump at the initial step - as noted in your earlier post. (Hence I will be removing them - its best not to introduce a potential cause of instability which may become a headache later on)

Regarding the last solution, I was curious as to why it works, since I already have the statement L = atpo in the BREAKPOINT, before the SOLVE statement. Hence I printed the values of all variables just after it enters the BREAKPOINT, after the L = atpo statement in BREAKPOINT, and after the SOLVE statement in BREAKPOINT. Doing so reveals that when the BREAKPOINT is called the first time, before the L = atpo statement is executed and the SOLVE statement is even reached - the state variables show changes. It was changes here that made the mechanism unstable. As L was 0 before L = atpo statement, the state variables just blew up on this first change - making the first run unstable. As L gets set on the second run, the instable behavior vanishes for all subsequent runs.
EDIT: I tried with changing atpo concentration. I see that the dependence between the runs is broken only when the L = atpo statement is introduced in the INITIAL block, as expected from above.

But is this behavior expected? As far as I had understood, on start of any run, neuron first calls the INITIAL block of all mechanisms, and then calls BREAKPOINT block for each mechanism at certain time points in each time step. Since the BREAKPOINT contains the statement L = atpo before the SOLVE statement, shouldn't changes in the state variables happen after the L = atpo statement has been executed - thus removing the possibility of such an instability? Or am I missing something obvious here?

ted
Site Admin
Posts: 5115
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Inconsistent behavior on modeling p2y2 Receptor

Post by ted » Tue Nov 14, 2017 9:30 am

I'll have to get back to you about this after the end of the SFN meeting.

Post Reply