update SI

Author: Deming-Yang

code/2 Driver Sr turnover.R

@@ -145,7 +145,21 @@ MCMC.CI.c$CI_high
 log(2)/MCMC.CI.c$CI_low
 log(2)/MCMC.CI.c$CI_high
 
+MAP.flux.ratio <- map_estimate(post.misha.pc2p3$BUGSoutput$sims.list$flux.ratio)
+MAP.flux.ratio[1]
+MCMC.CI.flux.ratio <- hdi(post.misha.pc2p3$BUGSoutput$sims.list$flux.ratio,0.89)
+MCMC.CI.flux.ratio$CI_low
+MCMC.CI.flux.ratio$CI_high
+
+MAP.pool.ratio <- map_estimate(post.misha.pc2p3$BUGSoutput$sims.list$pool.ratio)
+MAP.pool.ratio[1]
+MCMC.CI.pool.ratio <- hdi(post.misha.pc2p3$BUGSoutput$sims.list$pool.ratio,0.89)
+MCMC.CI.pool.ratio$CI_low
+MCMC.CI.pool.ratio$CI_high
+
 #get posterior of intake ratio
 post.misha.pc2p3.Rin.89<- MCMC.CI.bound(post.misha.pc2p3$BUGSoutput$sims.list$Rin, 0.89)
 #calculate the mean of intake after switch
 intake.af <- mean(post.misha.pc2p3.Rin.89[[1]][90:750])
+
+

code/5 Driver inversion case study mammoth.R

@@ -173,3 +173,17 @@ MCMC.CI.c.m$CI_low
 MCMC.CI.c.m$CI_high
 log(2)/MCMC.CI.c.m$CI_low
 log(2)/MCMC.CI.c.m$CI_high
+
+flux.ratio.m <- post.misha.invmamm.i$BUGSoutput$sims.list$a.m/post.misha.invmamm.i$BUGSoutput$sims.list$b.m
+MAP.flux.ratio.m <- map_estimate(flux.ratio.m)
+MAP.flux.ratio.m[1]
+MCMC.CI.flux.ratio.m <- hdi(flux.ratio.m,0.89)
+MCMC.CI.flux.ratio.m$CI_low
+MCMC.CI.flux.ratio.m$CI_high
+
+pool.ratio.m <- post.misha.invmamm.i$BUGSoutput$sims.list$c.m/post.misha.invmamm.i$BUGSoutput$sims.list$b.m
+MAP.pool.ratio.m <- map_estimate(pool.ratio.m)
+MAP.pool.ratio.m[1]
+MCMC.CI.pool.ratio.m <- hdi(pool.ratio.m,0.89)
+MCMC.CI.pool.ratio.m$CI_low
+MCMC.CI.pool.ratio.m$CI_high

code/6 Driver forward model.R

@@ -224,8 +224,7 @@ res60 <- forw.m(t = length(syn.input.60), input = syn.input.60, a = a, b = b, c
 a <- MAP.a[1]
 b <- MAP.b[1]
 c <- MAP.c[1]
-c/b
-
+c/b #0.2909254
 a/b
 #to change pool ratio, change c; to change flux ratio, change a
 # flux.ratio <- a/b
@@ -255,7 +254,7 @@ reac.prog.tk <- reac.prog[101:length(reac.prog)]
 lm.res.misha.2 <- lm(log(reac.prog.tk[200:800]) ~ c(200:800))
 summary(lm.res.misha.2)
 lm.res.misha.2$coefficients[[1]] #intercept
-lm.res.misha.2$coefficients[[2]] #slope
+lm.res.misha.2$coefficients[[2]] #slope lambda2
 
 exp(lm.res.misha.2$coefficients[[1]]) #f2
 -1/lm.res.misha.2$coefficients[[2]]*log(2) #t1/2 = 330 days
@@ -270,7 +269,7 @@ plot(c(1:150),reac.prog.res1,ylab="residual: ln(1-F)",xlab="Days")
 lm.res.misha.1 <- lm(reac.prog.res1 ~ c(1:150))
 summary(lm.res.misha.1)
 lm.res.misha.1$coefficients[[1]] #intercept
-lm.res.misha.1$coefficients[[2]] #slope
+lm.res.misha.1$coefficients[[2]] #slope lambda1
 
 exp(lm.res.misha.1$coefficients[[1]]) #f1 = 0.48
 -1/lm.res.misha.1$coefficients[[2]]*log(2) #t1/2 = 20.5 days
@@ -281,13 +280,17 @@ exp(lm.res.misha.1$coefficients[[1]]) #f1 = 0.48
 # f1/f2 ~= flux ratio
 exp(lm.res.misha.1$coefficients[[1]])/exp(lm.res.misha.2$coefficients[[1]])
 
+# lambda1/(lambda2*f1*f2) ~= pool ratio
+lm.res.misha.2$coefficients[[2]]/lm.res.misha.1$coefficients[[2]]/exp(lm.res.misha.1$coefficients[[1]])/exp(lm.res.misha.2$coefficients[[1]])
+
 #################use a different ratio###################
 ####begin forward model####
 #loading MAP estimates from posterior of the calibration
 a <- 0.04
 b <- MAP.b[1]
 c <- MAP.c[1]
 a/b
+c/b
 #to change pool ratio, change c; to change flux ratio, change a
 # flux.ratio <- a/b
 # pool.ratio <- c/b
@@ -317,7 +320,7 @@ lm.res.frl.2 <- lm(log(reac.prog.frl.tk[200:800]) ~ c(200:800))
 plot(1:length(reac.prog.frl.tk), log(reac.prog.frl.tk),ylab="ln(1-F)",xlab="Days",main="Slow turnover pool")
 summary(lm.res.frl.2)
 lm.res.frl.2$coefficients[[1]] #intercept
-lm.res.frl.2$coefficients[[2]] #slope
+lm.res.frl.2$coefficients[[2]] #slope lambda2
 
 exp(lm.res.frl.2$coefficients[[1]]) #f2 = 0.2885
 -1/lm.res.frl.2$coefficients[[2]]*log(2) #t1/2 = 271 days
@@ -331,7 +334,7 @@ plot(c(1:150),reac.prog.frl.res1,ylab="residual: ln(1-F)",xlab="Days")
 lm.res.frl.1 <- lm(reac.prog.frl.res1 ~ c(1:150))
 summary(lm.res.frl.1)
 lm.res.frl.1$coefficients[[1]] #intercept
-lm.res.frl.1$coefficients[[2]] #slope
+lm.res.frl.1$coefficients[[2]] #slope lambda1
 exp(lm.res.frl.1$coefficients[[1]]) #f1 = 0.711
 -1/lm.res.frl.1$coefficients[[2]]*log(2) #t1/2 = 12.3 days
 
@@ -340,3 +343,58 @@ exp(lm.res.frl.1$coefficients[[1]]) #f1 = 0.711
 
 # f1/f2 ~= flux ratio
 exp(lm.res.frl.1$coefficients[[1]]) /exp(lm.res.frl.2$coefficients[[1]])  
+
+# lambda1/(lambda2*f1*f2) ~= pool ratio
+lm.res.frl.2$coefficients[[2]]/lm.res.frl.1$coefficients[[2]]/exp(lm.res.frl.1$coefficients[[1]])/exp(lm.res.frl.2$coefficients[[1]])
+##############use a different c #################
+a <- MAP.a[1]
+b <- MAP.b[1]
+c <- 0.001 #MAP.c = 0.0041
+c/b
+
+#2 number of days in the simulation
+t <- 900
+
+#3 generate input series with fixed duration#
+input.misha <- initiate.switch(t, n.switch=1, day.switch=100, a=0.706, gap=0.005, duration=360)
+
+#4 generate serum and bone series based on input series and turnover parameters#
+res <- forw.m(t = 900, input = input.misha, a = a, b = b, c = c, R1.int = NULL, R2.int = NULL)
+res.prs <- res[[1]]
+
+reac.prog.prs <- (0.711-res.prs)/(0.711-0.706)
+
+reac.prog.prs.tk <- reac.prog.prs[101:length(reac.prog.prs)]
+
+#after the inflection point: [200:800]
+
+lm.res.prs.2 <- lm(log(reac.prog.prs.tk[200:800]) ~ c(200:800))
+plot(1:length(reac.prog.prs.tk), log(reac.prog.prs.tk),ylab="ln(1-F)",xlab="Days",main="Slow turnover pool")
+summary(lm.res.prs.2)
+lm.res.prs.2$coefficients[[1]] #intercept
+lm.res.prs.2$coefficients[[2]] #slope lambda2
+
+exp(lm.res.prs.2$coefficients[[1]]) #f2 = 0.4700503
+-1/lm.res.prs.2$coefficients[[2]]*log(2) #t1/2 = 1284.937 days
+
+#estimate of parameter c:
+-lm.res.prs.2$coefficients[[2]]/(1-exp(lm.res.prs.2$coefficients[[1]]))
+
+#fit for pool 1
+reac.prog.prs.res1 <- log(reac.prog.prs.tk[1:150]-exp(1:150 * lm.res.prs.2$coefficients[[2]] +lm.res.prs.2$coefficients[[1]])) #first residual
+plot(c(1:150),reac.prog.prs.res1,ylab="residual: ln(1-F)",xlab="Days")
+lm.res.prs.1 <- lm(reac.prog.prs.res1 ~ c(1:150))
+summary(lm.res.prs.1)
+lm.res.prs.1$coefficients[[1]] #intercept
+lm.res.prs.1$coefficients[[2]] #slope lambda1
+exp(lm.res.prs.1$coefficients[[1]]) #f1 = 0.5343916
+-1/lm.res.prs.1$coefficients[[2]]*log(2) #t1/2 = 21.53086 days
+
+#estimate of parameter a:
+-lm.res.prs.1$coefficients[[2]]*exp(lm.res.prs.1$coefficients[[1]]) 
+
+# fin/f2 ~= flux ratio
+exp(lm.res.prs.1$coefficients[[1]]) /exp(lm.res.prs.2$coefficients[[1]])  
+
+# lambda1/(lambda2*f1*f2) ~= pool ratio
+lm.res.prs.2$coefficients[[2]]/lm.res.prs.1$coefficients[[2]]/exp(lm.res.prs.1$coefficients[[1]])/exp(lm.res.prs.2$coefficients[[1]])