# # Apoptosis model # # Library of R ODE solvers library("deSolve") # # ODE routine setwd("c:/R/bme_ode/chap3") source("apoptosis_1.R") # # Multiple cases for(ncase in 1:1){ # # Model parameters if(ncase==1){ a_hif=1.52; a_o2=1.8; a_p53=0.05; a_3=0.9 ; a_4=0.2 ; a_5=0.001; a_7 =0.7 ; a_8=0.06; a_9 =0.1 ; a_10=0.7; a_11=0.2; a_12=0.1 ; a_13=0.1 ; a_14=0.05;} # # Initial condition nout=101; ncall=0 yini=c(1,0,0,0,0,0) # # t interval times=seq(from=0,to=100,by=100/(nout-1)) # # ODE integration out=ode(y=yini,times=times,func=apoptosis_1,parms=NULL) # # ODE numerical solution cat(sprintf("\n\n ncase = %2d",ncase)) for(it in 1:nout){ if((it==1)|(it==26)|(it==51)|(it==101)){ cat(sprintf("\n\n t = %4.1f",out[it,1])) cat(sprintf("\n y_hif = %8.4f y_o2 = %8.4f",out[it,2],out[it,3])) cat(sprintf("\n y_p300 = %8.4f y_p53 = %8.4f",out[it,4],out[it,5])) cat(sprintf("\n y_casp = %8.4f y_kp = %8.4f",out[it,6],out[it,7])) } } # # Calls to apoptosis_1 cat(sprintf("\n\n ncall = %5d\n\n",ncall)) # # Six plots for y_hif,y_o2,y_p300,y_p53,y_casp,y_kp vs t par(mfrow=c(3,2)) # # y_hif plot(out[,1],out[,2],xlab="t",ylab="y_hif(t)",xlim=c(0,100), type="l",lwd=2,main="Hypoxia inducible factor") # # y_o2 plot(out[,1],out[,3],xlab="t",ylab="y_o2(t)",xlim=c(0,100), type="l",lwd=2,main="Oxygen level") # # y_p300 plot(out[,1],out[,4],xlab="t",ylab="y_p300(t)",xlim=c(0,100), type="l",lwd=2,main="Co-activator") # # y_p53 plot(out[,1],out[,5],xlab="t",ylab="y_p53(t)",xlim=c(0,100), type="l",lwd=2,main="Tumor suppressor gene") # # y_casp plot(out[,1],out[,6],xlab="t",ylab="y_casp(t)",xlim=c(0,100), type="l",lwd=2,main="Caspases") # # y_k plot(out[,1],out[,7],xlab="t",ylab="y_kp(t)",xlim=c(0,100), type="l",lwd=2,main="Potassium ions") # # Next case }