# Simulating basic enzyme kinetics - no assumptions
#
#  Craig Martin    Chem 728   2021


# Wrap everything in a user function (not necessary, but easier to tweak below)

# the following sets defaults for some parameters (that then become optional)
EnzKin <- function(t0=0, tf, tincr=(tf-t0)/50, k1p, km1p, kcp, krp=0, E0, S0, ES0=0, P0=0) {

	# create a time range
	time <- seq(t0, tf, by = tincr)

	# parameters:
	parameters <- c(k1=k1p, km1=km1p, kc=kcp, kr=krp)

	# initial condition:
	state <- c(E=E0, S=S0, ES=ES0, P=P0)

	# R function to calculate the value of the derivatives at each time value
	# Use the names of the variables as defined in the vectors above
	# define in same order as specified above in state
	multiKin <- function(t, state, parameters){
	  with(as.list(c(state, parameters)), {
		dE <-  -k1*E*S + km1*ES + kc*ES - kr*E*P
		dS <-  -k1*E*S + km1*ES
		dES <- k1*E*S - km1*ES - kc*ES + kr*E*P    
		dP <-  kc*ES - kr*E*P
		
		# return derivatives in same relative order as specified above in state
		return(list(c(dE, dS, dES, dP)))
	  })
	}  ## end of function multiKin

	## Integration with 'ode' - ordinary differential equations
	out1 <- ode(y = state, times = time, func = multiKin, parms = parameters)

	# get results as a dataframe
	df <- as.data.frame(out1)
	
	Km <- (km1p+kcp)/k1p   # not used right now

	# return a list with things we might want to access directly in plotting, etc
	outl <- list("output"=out1, "t"=time, "E"=df$E, "S"=df$S, "ES"=df$ES, "P"=df$P )

	return(outl)
 }

# You can see from above, how you might output a "y-value" that could feed into nls

# call the above function (you can omit optional parameters)
 out <- EnzKin(tf=15,  tincr=0.01, k1p=0.1, km1p=0.1, kcp=0.2, E0=0.2, S0=20.0)


# plot the results (this will auto-scale to the product range; over-ride if desired)
plot( out$t, out$P, col="green", type="l", xlab="Time (s)", ylab="Concentration")
lines(out$t, out$E, col="purple")
lines(out$t, out$ES, col="red")
lines(out$t, out$S, col="blue")

legend(200, 3, c("P", "E", "ES", "S"), col = c("green", "purple", "red", "blue"), lty=c("solid","solid","solid","solid"))