#####################################
# Simple model of reaction kinetics
#####################################

# Scheme of reaction:
# E + S <-> C -> P + E
#
# E: free enzyme
# S: free substrate
# C: enzyme-substrate complex
# P: product


###################################
# FUNCTION DEFINITIONS
###################################

###
# reaction <- function(t, x, parms)
# Use: calculates the derivatives for the reaction kinetics model
# Input:
#      t: time (not used here, because there is no explicit time dependence)
#      x: vector of the current values of S and E
#      parms: vector of parameters
# Output:
#      der: vector of derivatives

# To use the lsoda function, the model function has to be a function of t (time),
# x (the vector of the variables) and parms (for parameters of the model).

reaction <- function(t, x, parms){
with(as.list(c(x,parms)),{
dS <- - k1*E*S + k2*(etot-E)
dE <- - k1*E*S + (k2+k3)*(etot-E)
der <- c(dS, dE)  # order should be the same as in x!
list(der)  # output must be returned

}) # end of 'with'
}  # end of function definition

###########################
# MAIN PROGRAM
###########################

### LOAD LIBRARIES
#load R library for ordinary differential equation solvers
library(odesolve)

# set the initial values
inits <- c(E=10, S=20, C=0, P=0)

# calculate derived parameters
etot = inits["E"]+inits["C"]
stot = sum(inits) - inits["E"]

dt    <- seq(0,50,0.1)  # set the time points for evaluation

### SIMULATE THE MODEL

## Use lsoda to solve the differential equations numerically. The syntax should be
## lsoda(initial values, time points, model function, parameters)

simulation <- as.data.frame(lsoda(inits[c("S","E")], dt, reaction, parms=c(k1=0.25,k2=0.2,k3=0.1))) # note that etot is defined globally 

### PLOT THE OUTPUT

attach(simulation) # this command allows you to refer to the columns of the data frame directly.

# calculate C and P
C <- etot - E
P <- stot-etot-S+E

pdf("reaction.pdf")
par(lwd=2)  # set default line width
plot(dt, S, type="l", col="blue", ylim=c(0,max(c(S,C,P,E))), xlab="time", ylab="concentrations")
lines(dt, E, type="l", col="red")
lines(dt, C, type="l", col="green")
lines(dt, P, type="l", col="brown")

# Add a legend to the graph
legend(max(dt)*0.7,max(inits)*0.8, legend=c("S","E","C","P"), col=c("blue", "red", "green","brown"), lty=1)
dev.off()

detach(simulation)

write.table(simulation,file="reaction.data")
copy <- read.table("reaction.data",header=TRUE)