#Fitting a Michaelis-Menten model
rm(list = ls())
##################################
library(drc) ## for fitting Michaelis Menten model
## Loading required package: MASS
##
## 'drc' has been loaded.
## Please cite R and 'drc' if used for a publication,
## for references type 'citation()' and 'citation('drc')'.
##
## Attaching package: 'drc'
## The following objects are masked from 'package:stats':
##
## gaussian, getInitial
library(ggplot2) ## for drawing
mm <- structure(list(S = c(3.6, 1.8, 0.9, 0.45, 0.225, 0.1125, 3.6,
1.8, 0.9, 0.45, 0.225, 0.1125, 3.6, 1.8, 0.9, 0.45, 0.225, 0.1125,
0), v = c(0.004407692, 0.004192308, 0.003553846, 0.002576923,
0.001661538, 0.001064286, 0.004835714, 0.004671429, 0.0039, 0.002857143,
0.00175, 0.001057143, 0.004907143, 0.004521429, 0.00375, 0.002764286,
0.001857143, 0.001121429, 0)),
.Names = c("S", "v"), class = "data.frame", row.names = c(NA, -19L))
mm[order(mm$S),]
## S v
## 19 0.0000 0.000000000
## 6 0.1125 0.001064286
## 12 0.1125 0.001057143
## 18 0.1125 0.001121429
## 5 0.2250 0.001661538
## 11 0.2250 0.001750000
## 17 0.2250 0.001857143
## 4 0.4500 0.002576923
## 10 0.4500 0.002857143
## 16 0.4500 0.002764286
## 3 0.9000 0.003553846
## 9 0.9000 0.003900000
## 15 0.9000 0.003750000
## 2 1.8000 0.004192308
## 8 1.8000 0.004671429
## 14 1.8000 0.004521429
## 1 3.6000 0.004407692
## 7 3.6000 0.004835714
## 13 3.6000 0.004907143
###################################method_1
#We will use package drm to fit the model
model.drm <- drm (v ~ S, data = mm, fct = MM.2())
summary(model.drm)
##
## Model fitted: Michaelis-Menten (2 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## d:(Intercept) 0.00542523 0.00011758 46.141 < 2.2e-16 ***
## e:(Intercept) 0.43980191 0.03051108 14.415 5.814e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.0001727146 (17 degrees of freedom)
mml <- data.frame(S = seq(0, max(mm$S), length.out = 100))
mml$v <- predict(model.drm, newdata = mml)
head(mml)
## S v
## 1 0.00000000 0.0000000000
## 2 0.03636364 0.0004143122
## 3 0.07272727 0.0007698341
## 4 0.10909091 0.0010782499
## 5 0.14545455 0.0013483403
## 6 0.18181818 0.0015868311
#Using the below code, we can visualize the result.
ggplot(mm, aes(x = S, y = v)) +
theme_bw() +
xlab("Concentration [mM]") +
ylab("Speed [dE/min]") +
ggtitle("Michaelis-Menten kinetics") +
geom_point(alpha = 0.5) +
geom_line(data = mml, aes(x = S, y = v), colour = "red")
########################################method_2
#Using nls()
#We need to provide the starting values Vm and K by ourself.
#Vm is the asymptote and K is half way between 0 and the asymptote.
model.nls <- nls(v ~ Vm * S/(K+S), data = mm,
start = list(K = max(mm$v)/2, Vm = max(mm$v)))
summary(model.nls)
##
## Formula: v ~ Vm * S/(K + S)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## K 0.4398016 0.0311612 14.11 8.1e-11 ***
## Vm 0.0054252 0.0001193 45.47 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0001727 on 17 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 4.666e-07
mml_2 <- data.frame(S = seq(0, max(mm$S), length.out = 100))
mml_2$v <- predict(model.nls, newdata = mml_2)
head(mml_2)
## S v
## 1 0.00000000 0.0000000000
## 2 0.03636364 0.0004143124
## 3 0.07272727 0.0007698344
## 4 0.10909091 0.0010782504
## 5 0.14545455 0.0013483408
## 6 0.18181818 0.0015868316
##############
ggplot(mm, aes(x = S, y = v)) +
theme_bw() +
xlab("Concentration [mM]") +
ylab("Speed [dE/min]") +
ggtitle("Michaelis-Menten kinetics") +
geom_point(alpha = 0.5) +
geom_line(data = mml_2, aes(x = S, y = v),colour = "green")
#To get V_{max} and K_M use coef() in R to get the coefficients of the model.
coef(model.nls) #
## K Vm
## 0.439801576 0.005425232
coef(model.drm)
## d:(Intercept) e:(Intercept)
## 0.005425233 0.439801912
#############################plot3
ggplot(mm, aes(x = S, y = v)) +
geom_point(col = "yellowgreen") +
geom_smooth(method = nls,col = "skyblue",
method.args = list(formula = y ~ Vmax * x / (Km + x),
start = list(Km = coef(model.nls)[1], Vmax = coef(model.nls)[2])), se = F) +
xlab("Concentration [mM]") +
ylab("Speed [dE/min]") +
ggtitle("Michaelis-Menten kinetics") +
theme(legend.position = "",
panel.background = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, size=0),
panel.grid = element_blank(),
axis.text.x = element_text(size= 12, color = "darkblue",
family = "sans",hjust = 0.5),
axis.text.y = element_text(size= 12, color = "darkblue",
family = "sans",vjust = 0.5,hjust = 0.5),
axis.title = element_text(size=12, color = "darkred",family = "sans"),
axis.ticks = element_line(size= 0.5),
axis.ticks.length = unit(3, "pt"))
## `geom_smooth()` using formula 'y ~ x'
#output p1
ggsave(paste0(Sys.Date(),"-plot21.tiff"),
plot = last_plot(), device = NULL,
scale = 1, width = 10, height = 10, units ="cm",dpi = 300,
limitsize = TRUE,compression = "lzw")
## `geom_smooth()` using formula 'y ~ x'
#for mean se plot see https://rpubs.com/TX-YXL/656451
##############ref and acknowledgement
#https://rpubs.com/RomanL/6752