Hi All
I am trying to implement a neural model in R from a paper but I have run
into a bit of a problem. The model is getting so complex that I can't
really test it in bits any more. At the moment I get this error:
Error in eval(expr, envir, enclos) : object 'mCaT_soma_AB' not found
I think I have been staring at it too long because I can't see what I have
done wrong. Can anyone perhaps spot the error I have made. The full code is
posted below.
Many thanks
Jannetta
# TODO: Add comment
#
# Author: a9912577
###############################################################################
library(deSolve)
ST <- function(time, init, parms) {
with(as.list(c(init, parms)),{
#functions to calculate activation m and inactivation h of the currents
mNax <- function(v) 1/(1+exp(-(v+24.7)/5.29));
taumNa <- function(v) 1.32 - (1.26/(1+exp(-v+120)/25));
hNax <- function(v) 1/(1+exp((v+48.9)/5.18));
tauhNa <- function(v) (0.67/(1+exp(-(v+62.9)/10))) *
(1.5+1/(1+exp(v+34.9)/3.6));
mCaTx <- function(v) 1/(1+exp(-(v+25)/7.2))
taumCaT <- function(v) 55-(49.5/(1+exp(-(v+58)/17)))
hCaTx <- function(v) 1/(1+exp((v+36)/7))
tauhCaT_AB <-function(v) 87.5-(75/(1+exp(-(v+50)/16.9)))
tauhCaT_PD <-function(v) 350-(76/(1+exp(-(v+50)/16.9)))
mCaSx <- function(v) 1/(1+exp(-(v+22)/8.5))
taumCaS <- function(v) 16-(13.1/(1+exp(-(v+25.1)/26.5)))
mNapx <- function(v) 1 / (1+exp(-(v+26.8)/8.2))
taumNap <- function(v) 19.8-(10.7/(1+exp(-(v+26.5)/86.)))
hNapx <- function(v) 1/1+exp((v+48.5)/4.8)
tauhNap <- function(v) 666-(379/(1+exp(-(v+33.6)/11.7)))
mhx <-function(v) 1/(1+exp((v+70)/6))
taumh <- function(v) 272+(1499/(1+exp(-(v+42.2)/8.73)))
mKx <- function(v) 1/(1+exp(-(v+14.2)/11.8));
taumK <- function(v) 7.2-(6.4/(1+exp(-(v+28.3)/19.2)))
# AB soma
iCaT_soma_AB <- gCaT_soma_AB * mCaT_soma_AB ^ 3 * hCaT_soma_AB * (v_soma -
ECaT_soma_AB)
iCaS_soma_AB <- gCaS_soma_AB * mCaS_soma_AB ^ 3 * (v_soma - ECaS_soma_AB)
iNap_soma_AB <- gNap_soma_AB * mNap_soma_AB ^ 3 * hNap_soma_AB * (v_soma -
ENap_soma_AB)
ih_soma_AB <- gh_soma_AB * mh_soma_AB ^ 3 * hh_soma_AB * (v_soma -
Eh_soma_AB)
iK_soma_AB <- gK_soma_AB * mK_soma_AB ^ 4 * mK_soma_AB * (v_soma -
EK_soma_AB)
iKCa_soma_AB <- gKCa_soma_AB * mKCa_soma_AB ^ 4 * (v_soma - EKCa_soma_AB)
# Total current for Calcium
totalICa <- iCaT_soma_AB + iCaS_soma_AB
# Differential equations
dCaConc_soma <- (-F_AB * totalICa - CaConc_soma + C0_AB)/tauCa_AB
mKCax_AB <- function(v, CaConc_soma)
(CaConc_soma/(CaConc_soma+30))*(1/(1+exp(-(v+51)/4)))
mKCax_PD <- function(v, CaConc_soma)
(CaConc_soma/(CaConc_soma+30))*(1/(1+epx(-(v+51)/8)))
taumKCa <-function(v) 90.3 - (75.09 / (1+exp(-(v+46)/22.7)))
mAx <- function(v) 1/(1+exp(-(v+27)/8.7))
taumA <- function(v) 11.6-(10.4/(1+exp(-(v+32.9)/15.2)))
hAx <- function(v) 1 / (1+exp((v+46.9)/4.9))
tauhA <- function(v) 38.6 - (29.2/(1+exp(-(v+38.9)/26.5)))
mProcx <- function(v) 1 / (1+exp(-(v+12)/3.05))
taumProc <- 0.5
# Currents as product of maximal conducatance(g), activation(m) and
inactivation(h)
# Driving force (v-E) where E is the reversal potential of the particular
ion
# AB axon
iNa_axon_AB <- gNa_axon_AB * mNa_axon ^ 3 * hNa_axon * (v - ENa_axon_AB)
iK_axon_AB <- gK_axon_AB * mK_axon ^ 4 * (v - EK_axon_AB)
iLeak_axon_AB <- gLeak_axon_AB * (v - ELeak_axon_AB)
dv <- (0 - iNa_axon_AB - iK_axon_AB - iLeak_axon_AB) / C_axon_AB
dmNa_axon_AB <- (mNax(v) - mNa_axon_AB)/taumNa(v)
dhNa_axon_AB <- (hNax(v) - hNa_axon_AB)/tauhNa(v)
dmK_axon_AB <- (mKx(v) - mK_axon_AB)/taumK(v)
dv_soma <- (I - iCaT_soma_AB - iCaS_soma_AB - iNap_soma_AB - ih_soma_AB -
iK_soma_AB - iKCa_soma_AB)
dmCaT_soma_AB <- (mCaTx(v_soma) - mCaT_soma_AB)/taumCaT(v_soma)
dhCaT_soma_AB <- (hCaTx(v_soma) - hCaT_soma_AB)/tauhCaT_AB(v_soma)
dmCaS_soma_AB <- (mCaSx(v_soma) - mCaS_soma_AB)/taumCaS(v_soma)
dmNap_soma_AB <- (mNapx(v_soma) - mNap_soma_AB)/taumNap(v_soma)
dhNap_soma_AB <- (hNapx(v_soma) - hNap_soma_AB)/tauhNap(v_soma)
dmh_soma_AB <- (mhx(v_soma) - mh_soma_AB)/taumh(v_soma)
dmK_soma_AB <- (mKx(v_soma) - mK_soma_AB)/taumK(v_soma)
dmKCa_soma_AB <- (mKCax_AB(v_soma,CaConc_soma) -
mKCa_soma_AB)/taumKCa(v_soma)
list(c(dv,dv_soma,dmNa_axon_AB, dhNa_axon_AB, dmK_axon_AB, dmCaT_soma_AB,
dhCaT_soma_AB, dmCaS_soma_AB, dmNap_soma_AB, dhNap_soma_AB, dmh_soma_AB,
dmK_soma_AB, dmKCa_soma_AB,dCaConc_soma))
})}
## Set initial state
init = c(dv=-55,dv_soma=-55,dmNa_axon_AB=0, dhNa_axon_AB=0, dmK_axon_AB=0,
dmCaT_soma_AB=0, dhCaT_soma_AB=0, dmCaS_soma_AB=0, dmNap_soma_AB=0,
dhNap_soma_AB=0, dmh_soma_AB=0, dmK_soma_AB=0,
dmKCa_soma_AB=0,dCaConc_soma=0)
## Set parameters
F_AB=0.418
C0_AB=0.5
tauCa_AB=303
ENa_axon_AB=50
EK_axon_AB=-80
ELeak_axon_AB=-60
gNa_axon_AB=0.300
gK_axon_AB=0.0525
gLeak_axon_AB=20
C_axon_AB=0.0015
ECaT_soma_AB=55.2e-3
ECaS_soma_AB=9e-3
ENap_soma_AB=50
Eh_soma_AB=-20
EK_soma_AB=-80
EKCa_soma_AB=-80
EA_soma_AB=-8
EP_soma_AB=0
ELeak_soma_AB=-50
gCaT_soma_AB=55.2e-3
gCaS_soma_AB=9e-3
gNap_soma_AB=2.7e-3
gh_soma_AB=0.00054
gK_soma_AB=0.0525
gKCa_soma_AB=0.600
gA_soma_AB=0.0216
gP_soma_AB=0.570
gLeak_soma_AB=0.000045
gAxial_soma_AB=0.003
gGap_AB=0.75e-3
C_soma_AB=9e-3
F_PD=0.515
C0_PD=0.5
tauCa_PD=300
I=6.5
parms =
c(ENa_axon_AB,EK_axon_AB,ELeak_axon_AB,gNa_axon_AB,gK_axon_AB,gLeak_axon_AB,C_axon_AB,ECaT_soma_AB,ENap_soma_AB,Eh_soma_AB,EK_soma_AB,EKCa_soma_AB,EA_soma_AB,EP_soma_AB,ELeak_soma_AB,gCaT_soma_AB,gCaS_soma_AB,gNap_soma_AB,gh_soma_AB,gK_soma_AB,gKCa_soma_AB,gA_soma_AB,gP_soma_AB,gLeak_soma_AB,gAxial_soma_AB,C0_AB,C0_PD)
## Set integrations times
times = seq(from=0, to=100, by = 0.25);
out<-ode(y=init, times=times, func=ST, parms=parms)
plot(out)
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
