Describing and simulating the ODEs representing a chemical network, and plotting the simulated output.

ChemReactionNet.R

# Install and load ODE solver package
install.packages("deSolve")
library(deSolve)

# Named vector of parameter values
parameters=c(k1=0.15,k2=0.075,V=1.0)

# Named vector of initial conditions
state=c(A=4,B=2,AB=0)

# Function describing ODEs
RxnNet=function(t,state,parameters){
	with(as.list(c(state,parameters)),{
		dA    <- -(k1/V)*A*B
		dB    <- -(k1/V)*A*B -(k2/V)*B
		dAB   <- (k1/V)*A*B
    		return(list(c(dA, dB, dAB)))
	})
}

# Report concentrations at this vector of times
times=seq(0,20,by=0.01)

# Simulate from the model, given initial conditions and parameter values
out=ode(y=state,times=times,func=RxnNet,parms=parameters)

# Plot output
plot(out[,1],out[,2],type="l",ylim=c(0,max(out[,2:length(colnames(out))])),xlab="Time (s)",ylab="Concentration (mmol/L)",lwd=3,cex.lab=1.5)
points(out[,1],out[,3],type="l",col="blue",lwd=3)
points(out[,1],out[,4],type="l",col="red",lwd=3)
legend("topright",colnames(out)[2:length(colnames(out))],col=c("black","blue","red"),lwd=3)