######################################## # Hormones and Homeostasis: Keeping Blood Glucose Concentrations Stable # Wicked Problems, 2018 # Authors: B. Sanft and A. Walter # Solutions to Exercises ######################################## # Load libraries library(deSolve) library(manipulate) library(fBasics) library(FME) ######################################## ######################################## # Exercise 3(b) # Sketch f(G) versus G (uptake term) k=4e-4; r1=9.8e-3; G=seq(0,200,length=1000); plot(G,k*G+r1,type="l",ylab="uptake coefficient") ######################################## # Exercise 3(c) # Sketch insulin secretion term versus G v=250; Gm=200; G=seq(0,400,length=1000); manipulate({ plot(G,v*G^n/(Gm^n+G^n),type="l",ylab="secretion term",ylim=c(0,250))}, n=slider(1,6,step=1)) ######################################## ######################################## # Exercise 4(b) # Plot solution to differential equation model to find homeostatic concentration gluc <- function(t,y,parameters) { with(as.list(c(y,parameters)), { G_ing=0 G_uptake = (k*G+r1)*I I_secretion = v*G^6/(Gm^6+G^6) I_degradation = r2*I dG = G_ing + G_rel - G_uptake dI = I_secretion - I_degradation list(c(dG, dI)) })} t = seq(0,100,length=100) parameters = c(G_rel=1,k=4e-4,r1=9.8e-3,Gm=200,v=250,r2=0.2) yinit = c(G=100,I=20) out = ode(y = yinit, times = t, func = gluc, parms=parameters) plot(out[,1],out[,2],type="l",xlab="time (minutes)",ylab="ECF concentration",ylim=c(0,150)) lines(out[,1],out[,3],col="blue") legend("topright",c("glucose","insulin"),lty=c(1,1),col=c("black","blue")) tail(out) # show last six rows of output ######################################## # Exercise 4(c) # Simulate model using a step function for G_ing(t) gluc <- function(t,y,parameters) { with(as.list(c(y,parameters)), { if (t<30) G_ing=3 else G_ing=0 G_uptake = (k*G+r1)*I I_secretion = v*G^6/(Gm^6+G^6) I_degradation = r2*I dG = G_ing + G_rel - G_uptake dI = I_secretion - I_degradation list(c(dG, dI)) })} out = ode(y = yinit, times = t, func = gluc, parms=parameters) plot(out[,1],out[,2],type="l",xlab="time (minutes)",ylab="ECF glucose concentration") ######################################## ######################################## # Exercise 5(a) # Compare step function and trangular funtion for G_ing gluc <- function(t,y,parameters) { with(as.list(c(y,parameters)), { # Constant function # if (t<30) # G_ing=(30000/140/30) # else # G_ing=0 # Triangular function if (t<15) G_ing=(2*30000/140/30)/15*t else if (t>15 && t<30) G_ing=-(2*30000/140/30)/15*(t-30) else G_ing=0 G_uptake = (k*G+r1)*I I_secretion = v*G^6/(Gm^6+G^6) I_degradation = r2*I dG = G_ing + G_rel - G_uptake dI = I_secretion - I_degradation list(c(dG, dI)) })} out = ode(y = yinit, times = t, func = gluc, parms=parameters) plot(out[,1],out[,2],type="l",ylim=c(0,220),xlab="time (minutes)",ylab="ECF concentration") lines(out[,1],out[,3],col="blue",type="l",xlab="time (minutes)",ylab="glucose concentration") legend("topright",c("glucose (mg/dl)","insulin (ng/dl)"),lty=c(1,1),col=c("black","blue")) ######################################## # Exercise 5(e) # Assume a candy bar is consumed at t=0 and again at t=20 gluc <- function(t,y,parameters) { with(as.list(c(y,parameters)), { # Triangular function - first candy bar if (t<15) G_ing1=(2*30000/140/30)/15*t else if (t>15 && t<30) G_ing1=-(2*30000/140/30)/15*(t-30) else G_ing1=0 # Triangular function - second candy bar if (t>20 && t<35) G_ing2=(2*30000/140/30)/15*(t-20) else if (t>35 && t<50) G_ing2=-(2*30000/140/30)/15*(t-50) else G_ing2=0 G_ing = G_ing1 + G_ing2; G_uptake = (k*G+r1)*I I_secretion = v*G^6/(Gm^6+G^6) I_degradation = r2*I dG = G_ing + G_rel - G_uptake dI = I_secretion - I_degradation list(c(dG, dI)) })} t = seq(0,200,length=1000) out = ode(y = yinit, times = t, func = gluc, parms=parameters) plot(out[,1],out[,2],type="l",ylim=c(0,220),xlab="time (minutes)",ylab="ECF concentration") lines(out[,1],out[,3],col="blue",type="l",xlab="time (minutes)",ylab="glucose concentration") legend("topright",c("glucose (mg/dl)","insulin (ng/dl)"),lty=c(1,1),col=c("black","blue")) ######################################## ######################################## # Exercise 6 # Complex carbohydrates gluc <- function(t,y,parameters) { with(as.list(c(y,parameters)), { if (t<120) G_ing=3 else G_ing=0 G_uptake = (k*G+r1)*I I_secretion = v*G^6/(Gm^6+G^6) I_degradation = r2*I dG = G_ing + G_rel - G_uptake dI = I_secretion - I_degradation list(c(dG, dI)) })} t = seq(0,200,length=1000) out = ode(y = yinit, times = t, func = gluc, parms=parameters) plot(out[,1],out[,2],type="l",ylim=c(0,180),xlab="time (minutes)",ylab="ECF concentration") lines(out[,1],out[,3],col="blue",type="l",xlab="time (minutes)",ylab="glucose concentration") legend("topright",c("glucose (mg/dl)","insulin (ng/dl)"),lty=c(1,1),col=c("black","blue")) ######################################## ######################################## # Exercise 7 # Fit model to data (30,000 mg glucose consumed) t=c(0,20,40,60,80,120) data=c(108,131,134,107,118,109) # Spring 2018 data # data=c(94,117,145,122,123,91) # Spring 2018 data derivs <- function(t,y,parms) { with(as.list(parms), { G = y[1]; I = y[2]; load=30000 Td=Ta; Vmax=load/Ta/140; # Divide by volume of ECF if (tTa && t<(Ta+Td)) G_ing=Vmax-Vmax/Td*(t-Ta) else G_ing=0 G_release = 1; G_uptake = (4e-4*G+9.8e-3)*I I_secretion = 250*G^6/(200^6+G^6); I_degradation = I*0.2; dG = G_ing + G_release - G_uptake dI = I_secretion - I_degradation list(c(dG, dI)) })} time=seq(0,140,length=100) g0=data[1] # Set initial glucose value to first glucose reading at t=0 load=30000; # Amount of glucos consumed #parms=c(Ta=30) #yinit = c(g0,20) # Initial values: insulin = 60 mU #out=ode(y = yinit, times = time, func = derivs, parms=parms) #plot(out[,1],out[,2],type="l",ylim=c(0,200)) #points(t,data) #lines(out[,1],out[,3],col="blue") ############################################################## # Create a slider for the unknown parameter, Ta. # Determine a reasonable guess for Ta. manipulate({ parameters=c(Ta=Ta) out=ode(y = c(g0,41.6), times = time, func = derivs, parms=parameters) plot(out[,1],out[,2],type="l",ylim=c(20,150)) lines(out[,1],out[,3],lty=2) points(t,data) }, Ta=slider(0,100)) #################################################### # Objective function (sum of squared errors) # Calculates squared residuals sse.gluc = function(parms0){ #function to calculate squared errors out=ode(y = c(g0,20), times = t, func = derivs, parms=c(Ta=parms0)) return((out[,2]-data)^2) # squared residuals } ################################################### # Optimization: Find Ta that minimizes error between model output and data parms0<-c(20) #initial guess fit = modFit(f=sse.gluc,p=parms0,lower=0,upper=100) fit$par # parameters ################################################### # Obtain model prediction plot(t,data,col='blue', pch=21,ylim=c(15,150),ylab="ECF concentration") Ta=fit$par[1] mod.pred = ode(c(g0,20),times=time,func=derivs,parms=c(Ta)) lines(time,mod.pred[,2],col="blue") #and plot as a line lines(time,mod.pred[,3],col="black") legend("topright",c("data","glucose (mg/dl)","insulin (ng/dl)"),lty=c(NA,1,1),col=c("blue","blue","black"),pch=c(1,NA,NA))