## Example in modeling second order kinetics, by
## David Nicolaides.
## On simulated data.
##############################
## load TIMP
##############################
library("TIMP")
##############################
## SIMULATE DATA
##############################
## set up the Example problem, a la in-situ UV-Vis spectroscopy of a simple
## reaction.
## A + 2B -> C + D, 2C -> E
cstart <- c(A = 1.0, B = 0.8, C = 0.0, D = 0.0, E = 0.0)
times <- c(seq(0,2, length=21), seq(3,10, length=8))
k <- c(kA = 0.5, k2C = 1)
## stoichiometry matrices
rsmatrix <- c(1,2,0,0,0,0,0,2,0,0)
smatrix <- c(-1,-2,1,1,0,0,0,-2,0,1)
concentrations <- calcD(k, times, cstart, rsmatrix, smatrix)
wavelengths <- seq(500, 700, by=2)
spectra <- matrix(nrow = length(wavelengths), ncol = length(cstart))
location <- c(550, 575, 625, 650, 675)
delta <- c(10, 10, 10, 10, 10)
spectra[, 1] <- exp( - log(2) *
(2 * (wavelengths - location[1])/delta[1])^2)
spectra[, 2] <- exp( - log(2) *
(2 * (wavelengths - location[2])/delta[2])^2)
spectra[, 3] <- exp( - log(2) *
(2 * (wavelengths - location[3])/delta[3])^2)
spectra[, 4] <- exp( - log(2) *
(2 * (wavelengths - location[4])/delta[4])^2)
spectra[, 5] <- exp( - log(2) *
(2 * (wavelengths - location[5])/delta[5])^2)
sigma <- .001
Psi_q <- concentrations %*% t(spectra) + sigma *
rnorm(dim(concentrations)[1] * dim(spectra)[1])
## store the simulated data in an object of class "dat"
kinetic_data <- dat(psi.df=Psi_q , x = times, nt = length(times),
x2 = wavelengths, nl = length(wavelengths))
##############################
## DEFINE MODEL
##############################
## starting values
kstart <- c(kA = 1, k2C = 0.5)
## model definition for 2nd order kinetics
kinetic_model <- initModel(mod_type = "kin", seqmod = FALSE,
kinpar = kstart,
numericalintegration = TRUE,
initialvals = cstart,
reactantstoichiometrymatrix = rsmatrix,
stoichiometrymatrix = smatrix )
##############################
## FIT INITIAL MODEL
## adding constraints to non-negativity of the
## spectra via the opt option nnls=TRUE
##############################
kinetic_fit <- fitModel(data=list(kinetic_data),
modspec = list(kinetic_model),
opt = kinopt(nnls = TRUE, iter=80,
selectedtraces = seq(1,kinetic_data@nl,by=2)))
## look at estimated parameters
parEst(kinetic_fit)
## various results
## concentrations
conRes <- getX(kinetic_fit)
matplot(times, conRes, type="b", col=1,pch=21, bg=1:5, xlab="time (sec)",
ylab="concentrations", main="Concentrations (2nd order kinetics)")
## spectra
specRes <- getCLP(kinetic_fit)
matplot(wavelengths, specRes, type="b", col=1,pch=21, bg=1:5,
xlab="wavelength (nm)",
ylab="amplitude", main="Spectra")
## see help(getResults) for how to get more results information from
## kinetic_fit
##############################
## CLEANUP GENERATED FILES
##############################
# This removes the files that were generated in this example
# (do not run this code if you wish to inspect the output)
file_list_cleanup = c(Sys.glob("*paramEst.txt"), Sys.glob("*.ps"), Sys.glob("Rplots*.pdf"))
# Iterate over the files and delete them if they exist
for (f in file_list_cleanup) {
if (file.exists(f)) {
unlink(f)
}
}
Run the code above in your browser using DataLab