i have been trying to model Ca diffusion in extracellular space containing neuron and glial cells in rat v1 area. i went through your mails on K accumulation in extracellular space to auatai and decided to use that approach. currently what i am doing is that i am assuming that there is extracellular space, neuron and glial cell spatially arranged in a grid with finite elements. The grid contains 256 elements (16x16). Neuron size is(10x10) , glial size is(3x3), Surrounding both of them and between them is extracellular space, with each element size as 2um*2um. The code is
Code: Select all
TITLE Calcium ion accumulation and diffusion
NEURON
{
SUFFIX cadifus
USEION ca READ cao,cai,ica WRITE cai,ica,cao
GLOBAL vol,Buffer0,cao_moles,elems
RANGE ipump,id,rho,cax,last_bpc,meanup0
}
DEFINE NANN 4
UNITS
{
(mol) = (1)
(molar) = (1/liter)
(mM) = (millimolar)
(um) = (micron)
(mA) = (milliamp)
FARADAY = (faraday) (10000 coulomb)
PI = (pi) (1)
}
PARAMETER
{
DFree = .6 (um2/ms)
diam (um)
ica (mA/cm2)
k1buf = 500.0 (/mM-ms)
k2buf = 0.5 (/ms)
k1=1.e10 (/mM-ms)
k2=50.e7 (/s)
k3=1.e10 (/s)
k4=5.e6 (/mM-ms)
cai0 = 0.000116 (mM)
cao0 = 1.2 (mM)
rho = 2.2727 (1) :intra/extracellular
area (um2)
flag
id
fact=13.8
}
CONSTANT
{
volo=3.168e-12 (liter) :volume of ecs surrounding neurom
evolo = 0.072e-12 (liter) :volume of each element
nvolo=7.2e-12 (litre) :vol of neuron
}
ASSIGNED
{
cax (mM)
last_bpc (1)
meanup0 (mM)
cai (mM)
vol[NANN] (1) : gets extra cm2 when multiplied by diam^2
ipump (mA/cm2)
last_ipump (mA/cm2)
gcai_init (mM)
cao_moles (mol)
elems[256] (mM) :THESE ARE THE 256 ELEMENTS OF THE 16x16 GRID
}
STATE
{
ca[NANN] (mM) <1.e-7> : ca[0] is equivalent to cai
cao (mM)
CaBuffer[NANN] (mM)
Buffer[NANN] (mM)
pump (mol/cm2) <1.e-3>
pumpca (mol/cm2) <1.e-15>
}
LOCAL totpump,kd,totbuf
INITIAL
{
totpump=0.2
pump=totpump/(1+1.e-18*k4*cao/k3)
pumpca=1.4e-22
ipump=0
flag = 0
totbuf=1.2
kd=k2buf/k1buf
FROM i=0 TO NANN-1
{
ca[i] = cai
CaBuffer[i] =(totbuf*ca[i])/(kd+ca[i])
Buffer[i] = totbuf - CaBuffer[i]
}
itotca = 0
ototca = 0
gcai_init = 0
cao_moles = cao*volo :for 160 it is 15.36e-9 mmoles
cai = cai0
cao = cao0
cax = cao0 + rho*(cai0 + CaBuffer[0]*0.9722) + rho*(0.00004)*pumpca/diam :0.9722 is sum0f vol[i]/(pi/4) for 0,1,2 annulus
meanup0=1.2
}
BREAKPOINT
{
VERBATIM
float sample_time, samplet2;
int i, j, n1, n2, n3, n4, x, p, neigh, neigh1;
//GIVEN BELOW ARE THE ELEMENT IDs OF THE 55 ELEMENTS OCCUPIED BY total EXTRACELLULAR SPACE//
int list[55]= {0,1,2,3,4,5,6,7,8,9,10,11, 16, 27,28, 29,30,31,32,43,47,48,59,63,64,75,79,80,91,92,93,94,95,96,107,112, 123,128,139,144,155, 160,171,176,177,178,179,180,181,182,183,184,185,186,187} ;
int list0[44]= {0,1,2,3,4,5,6,7,8,9,10,11, 16, 27,32,43,48,59,64,75,80,91,96,107,112, 123,128,139,144,155, 160,171,176,177,178,179,180,181,182,183,184,185,186,187} ; //these are elements occupied by ECS surrounding neuron
double del=20.0,neigh_ca;
double dc=0.0232,totca; //232um2/s divided by e4 (as timestep is 0.1ms and not s)
SOLVE state METHOD sparse
last_ipump=ipump
last_bpc=(4*rho/PI)*((cai*PI/4)+(CaBuffer[0]*vol[0])+(CaBuffer[1]*vol[1])+(CaBuffer[2]*vol[2])) + (rho*(0.00004)*pumpca/diam)
VERBATIM
if(cai<0.000116)
{
k1=0.0;
cai=0.000116;
}
else k1=1.e10;
if (flag>1)
{
sample_time = (t-(((t)-(int)(t))));
samplet2 = (t-(((t/1000)-(int)(t/1000))*1000));
if((int)id==0)
{
for(i=0;i<256;i++)
if (elems[i] < 0.0) elems[i]=0.0;
if((int)(t*100)==5) // INITIALIZE EXTRACELLULAR CALCIUM AT TIME = 5 MILLI SECONDS
{
for(i=0 ;i<256;i++) elems[i] = 0.0;
for(i=0;i<55;i++) elems[list[i]]=1.2;// INITIALIZE CA CONC IN ECS ELEMENTS TO 1.2 mM
}// END OF INITIALIZATION
neigh = 44;
for(i=0;i<neigh;i++) elems[list0[i]]= cao; //EXTRACELLULAR CALCIUM UPDATE
//********************DIFFUSION BEGIN***********
for(i=0;i<55;i++)
{
x=list[i];
n1=x-1;n2=x+1;n3=x-16;n4=x+16; //n1-n4 ARE THE ELEMENT IDs OF THE FOUR POSSIBLE NEIGHBORS OF X
neigh1=0; //neigh IS A COUNTER THAT COUNTS THE VALID NEIGHBORS OF X
neigh_ca=0.0; //neigh_ca IS A COUNTER THAT ACCUMULATES THE CALCIUM OF VALID NEIGHBORS
for(j=0;(j<55)&&(neigh1<=3);j++)
{
p=list[j];
if ((n1==p)||(n2==p)||(n3==p)||(n4==p))
{
neigh1++;
neigh_ca=neigh_ca+elems[p];
}
}
elems[x]=elems[x]+(1/del)*dc*(neigh_ca - neigh1*elems[x]);
}
//**********************DIFFUSION END************
float sumup0=0.0;// meanup0=0.0 ;
for(i=0;i<44;i++) sumup0 = sumup0 + elems[list0[i]];
meanup0 = sumup0/44.0;
cax = meanup0+last_bpc;
} //id==0 loop closed
flag = 0;
}
else flag = 2;
ENDVERBATIM
}
LOCAL coord_done
INITIAL {
if (coord_done == 0) {
coord_done = 1
coord()
}
: note Buffer gets set to Buffer0 automatically
: and CaBuffer gets set to 0 (Default value of CaBuffer0) as well
FROM i=0 TO NANN-1 {
ca[i] = cai
}
}
LOCAL frat[NANN] : gets extra cm when multiplied by diam
PROCEDURE coord() {
LOCAL r, dr2
: cylindrical coordinate system with constant annuli thickness to
: center of cell. Note however that the first annulus is half thickness
: so that the concentration is second order correct spatially at
: the membrane or exact edge of the cell.
: note ca[0] is at edge of cell
: ca[NANN-1] is at center of cell
r = 1/2 :starts at edge (half diam)
dr2 = r/(NANN-1)/2 :half thickness of annulus
vol[0] = 0
frat[0] = 2*r
FROM i=0 TO NANN-2 {
vol[i] = vol[i] + PI*(r-dr2/2)*2*dr2 :interior half
r = r - dr2
frat[i+1] = 2*PI*r/(2*dr2) :exterior edge of annulus
: divided by distance between centers
r = r - dr2
vol[i+1] = PI*(r+dr2/2)*2*dr2 :outer half of annulus
}
}
LOCAL dsq, dsqvol : can't define local variable in KINETIC block or use
: in COMPARTMENT
KINETIC state {
COMPARTMENT i, diam*diam*vol[i]*1(um) {ca CaBuffer Buffer}
COMPARTMENT (1.e10)*area {pump pumpca}
COMPARTMENT (1.e15)*volo {cao}
~ ca[0] << (-(ica-last_ipump)*PI*fact*diam*frat[0]*1(um)/(2*FARADAY))
FROM i=0 TO NANN-2 {
~ ca[i] <-> ca[i+1] (DFree*frat[i+1]*1(um), DFree*frat[i+1]*1(um))
}
dsq = diam*diam*1(um)
FROM i=0 TO NANN-1 {
dsqvol = dsq*vol[i]
~ ca[i] + Buffer[i] <-> CaBuffer[i] (k1buf*dsqvol,k2buf*dsqvol)
}
~ca[0] + pump <-> pumpca ((1.e-11)*k1*area, (1.e7)*k2*area)
~pumpca <-> pump + cao ((1.e7)*k3*area, (1.e-11)*k4*area)
ipump = 2*FARADAY*(f_flux-b_flux)/area
cai = ca[0]
cao = cax-(4*rho/PI)*((cai*PI/4)+(CaBuffer[0]*vol[0])+(CaBuffer[1]*vol[1])+(CaBuffer[2]*vol[2])) - (rho*(0.00004)*pumpca/diam)
}at line 151 in file cadiv.mod: SOLVE state METHOD sparse
Error at section location RSTsingle[0].soma[0](0.5)
Convergence not achieved in maximum number of iterations
/usr/local/nrn/i686/bin/nrniv: scopmath library error
near line 0
{counter1()}
^
fadvance()
advance()
step()
continuerun(30001)
and others
though the approach works absolutely fine if there is no glial cell. (there is no additional extracellular space apart from which is directly surrounding the neuron)
the code also works fine if there is no diffusion (which amounts to no additional extracellular space).
the problem is there in how i am writing cao.
for writing cao i did the following calculation based on your reply for K accumulation:
Since in calcium there are significant differences. More so, because we want diffusion to be included and hence cax will not be fixed. We also included concentrations of buffers and pumps attached to calcium in the calculations.
In calcium we will calculate cai from the ODE and calculate cao through algebraic equations. So if
rho = intracellular volume / extracellular volume, and
cai0, cao0, CaBuffer0[0],CaBuffer0[1], CaBuffer0[2], pumpca0 are initial concentrations then:
cai*intracellular_v + cao*extracellular_v + CaBuffer[0]*Volume[0] + CaBuffer[1]*Volume[1] + CaBuffer[2]*Volume[2] + pumpCa*area = cai0*intracellular_v + cao0*extracellular_v + CaBuffer0[0]*Volume[0] + CaBuffer0[1]*Volume[1] + CaBuffer0[2]*Volume[2] + pumpCa0*area
Since Volume = diam*diam*L*vol
where vol is the volume of the annulus i for a cylinder with diam=1 um and L=1um
intracellular_v = (PI/4)*diam*diam*L
area = PI*diam*L , here pumpca has units of mol/cm2 and area is um2 hence a factor of 1e-8 is required to convert the area to cm2. Also a factor of 1e3 is required to convert the moles to milimoles.
Replacing these values and dividing by intracellular_v
=> cai + cao/rho + (CaBuffer[0]*vol[0] + CaBuffer[1]*vol[1] + CaBuffer[2]*vol[3]) / (PI/4) + (1e-8)*(1e3)*pupmca*4/diam = cai0 + cao0/rho + (CaBuffer0[0]*vol[0] + CaBuffer0[1]*vol[1] + CaBuffer0[2]*vol[3]) / (PI/4) + (1e-8)*(1e3)*pupmca0*4/diam
=> cao/rho = cai0 + cao0/rho + (CaBuffer0[0]*vol[0] + CaBuffer0[1]*vol[1] + CaBuffer0[2]*vol[2]) / (PI/4) + (1e-8)*(1e3)*pupmca0*4/diam - [cai + (CaBuffer[0]*vol[0] + CaBuffer[1]*vol[1] + CaBuffer[2]*vol[3]) / (PI/4) + (1e-8)*(1e3)*pupmca*4/diam]
=> cao = [ cao0 + cai0*rho + (rho*4/PI) * (CaBuffer0[0]*vol[0] + CaBuffer0[1]*vol[1] + CaBuffer0[2]*vol[2]) + rho* (1e-8)*(1e3)*pupmca0*4/diam ] - rho*[cai + (CaBuffer[0]*vol[0] + CaBuffer[1]*vol[1] + CaBuffer[2]*vol[3]) / (PI/4) + (1e-8)*(1e3)*pupmca*4/diam]
=> cao = cax - rho*[cai + (CaBuffer[0]*vol[0] + CaBuffer[1]*vol[1] + CaBuffer[2]*vol[2]) / (PI/4) + (1e-8)*(1e3)*pupmca*4/diam] .................(1)
where
cax = cao0 + cai0*rho + (rho*4/PI) * (CaBuffer0[0]*vol[0] + CaBuffer0[1]*vol[1] + CaBuffer0[2]*vol[2]) + rho* (1e-8)*(1e3)*pupmca0*4/diam
So the method which we use in ca mechanism is as follows:
(a) we calculate cai, CaBuffer[0], CaBuffer[1], CaBuffer[2], pumpca based on ODE
(b) we calculate cao based on these values and cax value, using eq (1) above
(c) We update elements value to cao : elems[list] = cao
(d) carry out diffusion on the elements value
(e) calculate mean of the element value
(f) calculate new cax using mean element value as cao0 and using fresh concentrations of cai, buffers and pumps
(g) calculate cai, CaBuffer[0], CaBuffer[1], CaBuffer[2], pumpca as in step (a)
THe code above is based on this approach.(at least i think so)
SO where am i going wrong?
thanks a lot for your help
with regards
kukkoo161