# ODE vector field: y' = f(t,y;p)
GABAB_signaling_vf <- function(t, state, parameters) {
## constants
## parameter values
kf_R0 = parameters[1];
kf_R1 = parameters[2];
kr_R1 = parameters[3];
kf_R2 = parameters[4];
kr_R2 = parameters[5];
kf_R3 = parameters[6];
kf_R4 = parameters[7];
kf_R5 = parameters[8];
kr_R5 = parameters[9];
kf_R6 = parameters[10];
kr_R6 = parameters[11];
kf_R7 = parameters[12];
kr_R7 = parameters[13];
kf_R8 = parameters[14];
kf_R9 = parameters[15];
kf_R10 = parameters[16];
kr_R10 = parameters[17];
kf_R11 = parameters[18];
kr_R11 = parameters[19];
kf_R12 = parameters[20];
kr_R12 = parameters[21];
kf_R13 = parameters[22];
kr_R13 = parameters[23];
kf_R14 = parameters[24];
kr_R14 = parameters[25];
k_gaba = parameters[26];
## state variables
gaba = state[1];
gabaGABABR = state[2];
gabaGABABRGi = state[3];
gabaGABABRGibg = state[4];
GABABR = state[5];
GABABRGi = state[6];
Gi = state[7];
GiaGDP = state[8];
GiaGTP = state[9];
GiaGTPRGS = state[10];
Gibg = state[11];
GIRK = state[12];
GIRKGibg = state[13];
GIRKGibg2 = state[14];
GIRKGibg3 = state[15];
GIRKGibg4 = state[16];
RGS = state[17];
VGCC = state[18];
VGCCGibg = state[19];
## expressions
R0 = k_gaba - kf_R0*gaba;
R1 = kf_R1*gaba*GABABR-kr_R1*gabaGABABR;
R2 = kf_R2*gabaGABABR*Gi-kr_R2*gabaGABABRGi;
R3 = kf_R3*gabaGABABRGi;
R4 = kf_R4*gabaGABABRGibg;
R5 = kf_R5*GABABR*Gi-kr_R5*GABABRGi;
R6 = kf_R6*gaba*GABABRGi-kr_R6*gabaGABABRGi;
R7 = kf_R7*GiaGTP*RGS-kr_R7*GiaGTPRGS;
R8 = kf_R8*GiaGTPRGS;
R9 = kf_R9*GiaGDP*Gibg;
R10 = kf_R10*GIRK*Gibg-kr_R10*GIRKGibg;
R11 = kf_R11*GIRKGibg*Gibg-kr_R11*GIRKGibg2;
R12 = kf_R12*GIRKGibg2*Gibg-kr_R12*GIRKGibg3;
R13 = kf_R13*GIRKGibg3*Gibg-kr_R13*GIRKGibg4;
R14 = kf_R14*VGCC*Gibg-kr_R14*VGCCGibg;
f_ <- numeric(19)
f_[1] <- +R0-R1-R6
f_[2] <- +R1-R2+R4
f_[3] <- +R2-R3+R6
f_[4] <- +R3-R4
f_[5] <- -R1-R5
f_[6] <- +R5-R6
f_[7] <- -R2-R5+R9
f_[8] <- +R8-R9
f_[9] <- +R3-R7
f_[10] <- +R7-R8
f_[11] <- +R4-R9-R10-R11-R12-R13-R14
f_[12] <- -R10
f_[13] <- +R10-R11
f_[14] <- +R11-R12
f_[15] <- +R12-R13
f_[16] <- +R13
f_[17] <- -R7+R8
f_[18] <- -R14
f_[19] <- +R14
names(f_) <- c("gaba", "gabaGABABR", "gabaGABABRGi", "gabaGABABRGibg", "GABABR", "GABABRGi", "Gi", "GiaGDP", "GiaGTP", "GiaGTPRGS", "Gibg", "GIRK", "GIRKGibg", "GIRKGibg2", "GIRKGibg3", "GIRKGibg4", "RGS", "VGCC", "VGCCGibg")
return(list(f_))
}
# ODE Jacobian: df(t,y;p)/dy
GABAB_signaling_jac <- function(t, state, parameters) {
## constants
## parameter values
kf_R0 = parameters[1];
kf_R1 = parameters[2];
kr_R1 = parameters[3];
kf_R2 = parameters[4];
kr_R2 = parameters[5];
kf_R3 = parameters[6];
kf_R4 = parameters[7];
kf_R5 = parameters[8];
kr_R5 = parameters[9];
kf_R6 = parameters[10];
kr_R6 = parameters[11];
kf_R7 = parameters[12];
kr_R7 = parameters[13];
kf_R8 = parameters[14];
kf_R9 = parameters[15];
kf_R10 = parameters[16];
kr_R10 = parameters[17];
kf_R11 = parameters[18];
kr_R11 = parameters[19];
kf_R12 = parameters[20];
kr_R12 = parameters[21];
kf_R13 = parameters[22];
kr_R13 = parameters[23];
kf_R14 = parameters[24];
kr_R14 = parameters[25];
k_gaba = parameters[26];
## state variables
gaba = state[1];
gabaGABABR = state[2];
gabaGABABRGi = state[3];
gabaGABABRGibg = state[4];
GABABR = state[5];
GABABRGi = state[6];
Gi = state[7];
GiaGDP = state[8];
GiaGTP = state[9];
GiaGTPRGS = state[10];
Gibg = state[11];
GIRK = state[12];
GIRKGibg = state[13];
GIRKGibg2 = state[14];
GIRKGibg3 = state[15];
GIRKGibg4 = state[16];
RGS = state[17];
VGCC = state[18];
VGCCGibg = state[19];
## expressions
R0 = k_gaba - kf_R0*gaba;
R1 = kf_R1*gaba*GABABR-kr_R1*gabaGABABR;
R2 = kf_R2*gabaGABABR*Gi-kr_R2*gabaGABABRGi;
R3 = kf_R3*gabaGABABRGi;
R4 = kf_R4*gabaGABABRGibg;
R5 = kf_R5*GABABR*Gi-kr_R5*GABABRGi;
R6 = kf_R6*gaba*GABABRGi-kr_R6*gabaGABABRGi;
R7 = kf_R7*GiaGTP*RGS-kr_R7*GiaGTPRGS;
R8 = kf_R8*GiaGTPRGS;
R9 = kf_R9*GiaGDP*Gibg;
R10 = kf_R10*GIRK*Gibg-kr_R10*GIRKGibg;
R11 = kf_R11*GIRKGibg*Gibg-kr_R11*GIRKGibg2;
R12 = kf_R12*GIRKGibg2*Gibg-kr_R12*GIRKGibg3;
R13 = kf_R13*GIRKGibg3*Gibg-kr_R13*GIRKGibg4;
R14 = kf_R14*VGCC*Gibg-kr_R14*VGCCGibg;
jac_ <- matrix(0.0,19,19)
jac_[1,1] <- -(kf_R0+kf_R1*GABABR+kf_R6*GABABRGi)
jac_[2,1] <- kf_R1*GABABR
jac_[3,1] <- kf_R6*GABABRGi
jac_[5,1] <- -kf_R1*GABABR
jac_[6,1] <- -kf_R6*GABABRGi
jac_[1,2] <- kr_R1
jac_[2,2] <- -(kr_R1+kf_R2*Gi)
jac_[3,2] <- kf_R2*Gi
jac_[5,2] <- kr_R1
jac_[7,2] <- -kf_R2*Gi
jac_[1,3] <- kr_R6
jac_[2,3] <- kr_R2
jac_[3,3] <- -(kr_R6+kr_R2+kf_R3)
jac_[4,3] <- kf_R3
jac_[6,3] <- kr_R6
jac_[7,3] <- kr_R2
jac_[9,3] <- kf_R3
jac_[2,4] <- kf_R4
jac_[4,4] <- -kf_R4
jac_[11,4] <- kf_R4
jac_[1,5] <- -kf_R1*gaba
jac_[2,5] <- kf_R1*gaba
jac_[5,5] <- -(kf_R1*gaba+kf_R5*Gi)
jac_[6,5] <- kf_R5*Gi
jac_[7,5] <- -kf_R5*Gi
jac_[1,6] <- -kf_R6*gaba
jac_[3,6] <- kf_R6*gaba
jac_[5,6] <- kr_R5
jac_[6,6] <- -(kr_R5+kf_R6*gaba)
jac_[7,6] <- kr_R5
jac_[2,7] <- -kf_R2*gabaGABABR
jac_[3,7] <- kf_R2*gabaGABABR
jac_[5,7] <- -kf_R5*GABABR
jac_[6,7] <- kf_R5*GABABR
jac_[7,7] <- -(kf_R2*gabaGABABR+kf_R5*GABABR)
jac_[7,8] <- kf_R9*Gibg
jac_[8,8] <- -kf_R9*Gibg
jac_[11,8] <- -kf_R9*Gibg
jac_[9,9] <- -kf_R7*RGS
jac_[10,9] <- kf_R7*RGS
jac_[17,9] <- -kf_R7*RGS
jac_[8,10] <- kf_R8
jac_[9,10] <- kr_R7
jac_[10,10] <- -(kr_R7+kf_R8)
jac_[17,10] <- kr_R7+kf_R8
jac_[7,11] <- kf_R9*GiaGDP
jac_[8,11] <- -kf_R9*GiaGDP
jac_[11,11] <- -(kf_R9*GiaGDP+kf_R10*GIRK+kf_R11*GIRKGibg+kf_R12*GIRKGibg2+kf_R13*GIRKGibg3+kf_R14*VGCC)
jac_[12,11] <- -kf_R10*GIRK
jac_[13,11] <- kf_R10*GIRK-kf_R11*GIRKGibg
jac_[14,11] <- kf_R11*GIRKGibg-kf_R12*GIRKGibg2
jac_[15,11] <- kf_R12*GIRKGibg2-kf_R13*GIRKGibg3
jac_[16,11] <- kf_R13*GIRKGibg3
jac_[18,11] <- -kf_R14*VGCC
jac_[19,11] <- kf_R14*VGCC
jac_[11,12] <- -kf_R10*Gibg
jac_[12,12] <- -kf_R10*Gibg
jac_[13,12] <- kf_R10*Gibg
jac_[11,13] <- kr_R10-kf_R11*Gibg
jac_[12,13] <- kr_R10
jac_[13,13] <- -(kr_R10+kf_R11*Gibg)
jac_[14,13] <- kf_R11*Gibg
jac_[11,14] <- kr_R11-kf_R12*Gibg
jac_[13,14] <- kr_R11
jac_[14,14] <- -(kr_R11+kf_R12*Gibg)
jac_[15,14] <- kf_R12*Gibg
jac_[11,15] <- kr_R12-kf_R13*Gibg
jac_[14,15] <- kr_R12
jac_[15,15] <- -(kr_R12+kf_R13*Gibg)
jac_[16,15] <- kf_R13*Gibg
jac_[11,16] <- kr_R13
jac_[15,16] <- kr_R13
jac_[16,16] <- -kr_R13
jac_[9,17] <- -kf_R7*GiaGTP
jac_[10,17] <- kf_R7*GiaGTP
jac_[17,17] <- -kf_R7*GiaGTP
jac_[11,18] <- -kf_R14*Gibg
jac_[18,18] <- -kf_R14*Gibg
jac_[19,18] <- kf_R14*Gibg
jac_[11,19] <- kr_R14
jac_[18,19] <- kr_R14
jac_[19,19] <- -kr_R14
## no names for return value
return(jac_)
}
# ODE parameter Jacobian: df(t,y;p)/dp
GABAB_signaling_jacp <- function(t, state, parameters) {
## constants
## parameter values
kf_R0 = parameters[1];
kf_R1 = parameters[2];
kr_R1 = parameters[3];
kf_R2 = parameters[4];
kr_R2 = parameters[5];
kf_R3 = parameters[6];
kf_R4 = parameters[7];
kf_R5 = parameters[8];
kr_R5 = parameters[9];
kf_R6 = parameters[10];
kr_R6 = parameters[11];
kf_R7 = parameters[12];
kr_R7 = parameters[13];
kf_R8 = parameters[14];
kf_R9 = parameters[15];
kf_R10 = parameters[16];
kr_R10 = parameters[17];
kf_R11 = parameters[18];
kr_R11 = parameters[19];
kf_R12 = parameters[20];
kr_R12 = parameters[21];
kf_R13 = parameters[22];
kr_R13 = parameters[23];
kf_R14 = parameters[24];
kr_R14 = parameters[25];
k_gaba = parameters[26];
## state variables
gaba = state[1];
gabaGABABR = state[2];
gabaGABABRGi = state[3];
gabaGABABRGibg = state[4];
GABABR = state[5];
GABABRGi = state[6];
Gi = state[7];
GiaGDP = state[8];
GiaGTP = state[9];
GiaGTPRGS = state[10];
Gibg = state[11];
GIRK = state[12];
GIRKGibg = state[13];
GIRKGibg2 = state[14];
GIRKGibg3 = state[15];
GIRKGibg4 = state[16];
RGS = state[17];
VGCC = state[18];
VGCCGibg = state[19];
## expressions
R0 = k_gaba - kf_R0*gaba;
R1 = kf_R1*gaba*GABABR-kr_R1*gabaGABABR;
R2 = kf_R2*gabaGABABR*Gi-kr_R2*gabaGABABRGi;
R3 = kf_R3*gabaGABABRGi;
R4 = kf_R4*gabaGABABRGibg;
R5 = kf_R5*GABABR*Gi-kr_R5*GABABRGi;
R6 = kf_R6*gaba*GABABRGi-kr_R6*gabaGABABRGi;
R7 = kf_R7*GiaGTP*RGS-kr_R7*GiaGTPRGS;
R8 = kf_R8*GiaGTPRGS;
R9 = kf_R9*GiaGDP*Gibg;
R10 = kf_R10*GIRK*Gibg-kr_R10*GIRKGibg;
R11 = kf_R11*GIRKGibg*Gibg-kr_R11*GIRKGibg2;
R12 = kf_R12*GIRKGibg2*Gibg-kr_R12*GIRKGibg3;
R13 = kf_R13*GIRKGibg3*Gibg-kr_R13*GIRKGibg4;
R14 = kf_R14*VGCC*Gibg-kr_R14*VGCCGibg;
jacp_ <- matrix(0.0,19,26)
jacp_[1,1] <- -gaba
jacp_[1,2] <- -gaba*GABABR
jacp_[2,2] <- gaba*GABABR
jacp_[5,2] <- -gaba*GABABR
jacp_[1,3] <- gabaGABABR
jacp_[2,3] <- -gabaGABABR
jacp_[5,3] <- gabaGABABR
jacp_[2,4] <- -gabaGABABR*Gi
jacp_[3,4] <- gabaGABABR*Gi
jacp_[7,4] <- -gabaGABABR*Gi
jacp_[2,5] <- gabaGABABRGi
jacp_[3,5] <- -gabaGABABRGi
jacp_[7,5] <- gabaGABABRGi
jacp_[3,6] <- -gabaGABABRGi
jacp_[4,6] <- gabaGABABRGi
jacp_[9,6] <- gabaGABABRGi
jacp_[2,7] <- gabaGABABRGibg
jacp_[4,7] <- -gabaGABABRGibg
jacp_[11,7] <- gabaGABABRGibg
jacp_[5,8] <- -GABABR*Gi
jacp_[6,8] <- GABABR*Gi
jacp_[7,8] <- -GABABR*Gi
jacp_[5,9] <- GABABRGi
jacp_[6,9] <- -GABABRGi
jacp_[7,9] <- GABABRGi
jacp_[1,10] <- -gaba*GABABRGi
jacp_[3,10] <- gaba*GABABRGi
jacp_[6,10] <- -gaba*GABABRGi
jacp_[1,11] <- gabaGABABRGi
jacp_[3,11] <- -gabaGABABRGi
jacp_[6,11] <- gabaGABABRGi
jacp_[9,12] <- -GiaGTP*RGS
jacp_[10,12] <- GiaGTP*RGS
jacp_[17,12] <- -GiaGTP*RGS
jacp_[9,13] <- GiaGTPRGS
jacp_[10,13] <- -GiaGTPRGS
jacp_[17,13] <- GiaGTPRGS
jacp_[8,14] <- GiaGTPRGS
jacp_[10,14] <- -GiaGTPRGS
jacp_[17,14] <- GiaGTPRGS
jacp_[7,15] <- GiaGDP*Gibg
jacp_[8,15] <- -GiaGDP*Gibg
jacp_[11,15] <- -GiaGDP*Gibg
jacp_[11,16] <- -GIRK*Gibg
jacp_[12,16] <- -GIRK*Gibg
jacp_[13,16] <- GIRK*Gibg
jacp_[11,17] <- GIRKGibg
jacp_[12,17] <- GIRKGibg
jacp_[13,17] <- -GIRKGibg
jacp_[11,18] <- -GIRKGibg*Gibg
jacp_[13,18] <- -GIRKGibg*Gibg
jacp_[14,18] <- GIRKGibg*Gibg
jacp_[11,19] <- GIRKGibg2
jacp_[13,19] <- GIRKGibg2
jacp_[14,19] <- -GIRKGibg2
jacp_[11,20] <- -GIRKGibg2*Gibg
jacp_[14,20] <- -GIRKGibg2*Gibg
jacp_[15,20] <- GIRKGibg2*Gibg
jacp_[11,21] <- GIRKGibg3
jacp_[14,21] <- GIRKGibg3
jacp_[15,21] <- -GIRKGibg3
jacp_[11,22] <- -GIRKGibg3*Gibg
jacp_[15,22] <- -GIRKGibg3*Gibg
jacp_[16,22] <- GIRKGibg3*Gibg
jacp_[11,23] <- GIRKGibg4
jacp_[15,23] <- GIRKGibg4
jacp_[16,23] <- -GIRKGibg4
jacp_[11,24] <- -VGCC*Gibg
jacp_[18,24] <- -VGCC*Gibg
jacp_[19,24] <- VGCC*Gibg
jacp_[11,25] <- VGCCGibg
jacp_[18,25] <- VGCCGibg
jacp_[19,25] <- -VGCCGibg
jacp_[1,26] <- 1
## no names for return value
return(jacp_)
}
# Output Function (Observables)
GABAB_signaling_func <- function(t, state, parameters) {
## constants
## parameter values
kf_R0 = parameters[1];
kf_R1 = parameters[2];
kr_R1 = parameters[3];
kf_R2 = parameters[4];
kr_R2 = parameters[5];
kf_R3 = parameters[6];
kf_R4 = parameters[7];
kf_R5 = parameters[8];
kr_R5 = parameters[9];
kf_R6 = parameters[10];
kr_R6 = parameters[11];
kf_R7 = parameters[12];
kr_R7 = parameters[13];
kf_R8 = parameters[14];
kf_R9 = parameters[15];
kf_R10 = parameters[16];
kr_R10 = parameters[17];
kf_R11 = parameters[18];
kr_R11 = parameters[19];
kf_R12 = parameters[20];
kr_R12 = parameters[21];
kf_R13 = parameters[22];
kr_R13 = parameters[23];
kf_R14 = parameters[24];
kr_R14 = parameters[25];
k_gaba = parameters[26];
## state variables
gaba = state[1];
gabaGABABR = state[2];
gabaGABABRGi = state[3];
gabaGABABRGibg = state[4];
GABABR = state[5];
GABABRGi = state[6];
Gi = state[7];
GiaGDP = state[8];
GiaGTP = state[9];
GiaGTPRGS = state[10];
Gibg = state[11];
GIRK = state[12];
GIRKGibg = state[13];
GIRKGibg2 = state[14];
GIRKGibg3 = state[15];
GIRKGibg4 = state[16];
RGS = state[17];
VGCC = state[18];
VGCCGibg = state[19];
## expressions
R0 = k_gaba - kf_R0*gaba;
R1 = kf_R1*gaba*GABABR-kr_R1*gabaGABABR;
R2 = kf_R2*gabaGABABR*Gi-kr_R2*gabaGABABRGi;
R3 = kf_R3*gabaGABABRGi;
R4 = kf_R4*gabaGABABRGibg;
R5 = kf_R5*GABABR*Gi-kr_R5*GABABRGi;
R6 = kf_R6*gaba*GABABRGi-kr_R6*gabaGABABRGi;
R7 = kf_R7*GiaGTP*RGS-kr_R7*GiaGTPRGS;
R8 = kf_R8*GiaGTPRGS;
R9 = kf_R9*GiaGDP*Gibg;
R10 = kf_R10*GIRK*Gibg-kr_R10*GIRKGibg;
R11 = kf_R11*GIRKGibg*Gibg-kr_R11*GIRKGibg2;
R12 = kf_R12*GIRKGibg2*Gibg-kr_R12*GIRKGibg3;
R13 = kf_R13*GIRKGibg3*Gibg-kr_R13*GIRKGibg4;
R14 = kf_R14*VGCC*Gibg-kr_R14*VGCCGibg;
func_ <- numeric(3)
func_[1] <- GIRKGibg4
func_[2] <- VGCCGibg
func_[3] <- gaba
names(func_) <- c("oGIRKGibg4", "oVGCCGibg", "ogaba")
return(func_)
}
GABAB_signaling_default <- function(t) {
## constants
parameters <- numeric(26)
parameters[1] <- 0.0535168989262411
parameters[2] <- 0.00285888595867763
parameters[3] <- 0.314477455454539
parameters[4] <- 0.0251449235770638
parameters[5] <- 0.0419082059617731
parameters[6] <- 0.0209541029808865
parameters[7] <- 0.167632823847092
parameters[8] <- 0.0125724617885319
parameters[9] <- 0.0209541029808865
parameters[10] <- 0.00285888595867763
parameters[11] <- 0.314477455454539
parameters[12] <- 2e-06
parameters[13] <- 0.002
parameters[14] <- 0.03
parameters[15] <- 0.00125
parameters[16] <- 0.000154358980053402
parameters[17] <- 0.0669308830676273
parameters[18] <- 0.000154358980053402
parameters[19] <- 0.0669308830676273
parameters[20] <- 0.000154358980053402
parameters[21] <- 0.0669308830676273
parameters[22] <- 0.000154358980053402
parameters[23] <- 0.0669308830676273
parameters[24] <- 0.000154358980053402
parameters[25] <- 0.0669308830676273
names(parameters) <- c("kf_R0", "kf_R1", "kr_R1", "kf_R2", "kr_R2", "kf_R3", "kf_R4", "kf_R5", "kr_R5", "kf_R6", "kr_R6", "kf_R7", "kr_R7", "kf_R8", "kf_R9", "kf_R10", "kr_R10", "kf_R11", "kr_R11", "kf_R12", "kr_R12", "kf_R13", "kr_R13", "kf_R14", "kr_R14", "k_gaba")
return(parameters)
}
GABAB_signaling_init <- function(t, parameters) {
## constants
## parameter values
kf_R0 = parameters[1];
kf_R1 = parameters[2];
kr_R1 = parameters[3];
kf_R2 = parameters[4];
kr_R2 = parameters[5];
kf_R3 = parameters[6];
kf_R4 = parameters[7];
kf_R5 = parameters[8];
kr_R5 = parameters[9];
kf_R6 = parameters[10];
kr_R6 = parameters[11];
kf_R7 = parameters[12];
kr_R7 = parameters[13];
kf_R8 = parameters[14];
kf_R9 = parameters[15];
kf_R10 = parameters[16];
kr_R10 = parameters[17];
kf_R11 = parameters[18];
kr_R11 = parameters[19];
kf_R12 = parameters[20];
kr_R12 = parameters[21];
kf_R13 = parameters[22];
kr_R13 = parameters[23];
kf_R14 = parameters[24];
kr_R14 = parameters[25];
k_gaba = parameters[26];
state <- numeric(19)
state[5] <- 400
state[7] <- 2600
state[17] <- 1000
names(state) <- c("gaba", "gabaGABABR", "gabaGABABRGi", "gabaGABABRGibg", "GABABR", "GABABRGi", "Gi", "GiaGDP", "GiaGTP", "GiaGTPRGS", "Gibg", "GIRK", "GIRKGibg", "GIRKGibg2", "GIRKGibg3", "GIRKGibg4", "RGS", "VGCC", "VGCCGibg")
return(state)
}