Non-linear least square
Biochemistry is an eminently quantitative science. The study of the relationship between variables and the determination of the parameters governing such relation, are part of the work of our discipline. Generally, the variables of interest are not linearly related to each other. Although exist linearizing transformations, these are not without risks since they can introduce significant biases.
The existence of programming languages such as R, allows us to easily and reliably address the non-linear fit of biochemical models. The function dir.MM() from the renz package carry out the non-linear least square fitting of kinetic data to the Michaelis-Menten equation.
\[\begin{equation}
\tag{1}
v = V_{max} \frac{[S]}{K_m + [S]}
\end{equation}\]
Regression analysis is a set of statistical processes for estimating the relationships between a dependent variable (in our case initial rate) and one or more independent variables (in our case the substrate concentration). The method least squares computes the unique line (or hyperplane) that minimizes the sum of squared differences between the true data and that line (see Figure).
\[\begin{equation}
\tag{2}
SSQ = \sum_{i=1}^n[v_i - V_{max} \frac{[S_i]}{K_m + [S]_i}]^2
\end{equation}\]
To find the kinetic parameters (\(K_m\) and \(V_{max}\)) that minimize the SSQ, we will have to solve the system of equations:
\[\begin{equation}
\tag{3}
\frac{\partial SSQ}{\partial V_{max}} = 0 = -2\sum_{i=1}^n[v_i - V_{max} \frac{[S]_i}{K_m + [S]_i}] [S]_i
\end{equation}\]
\[\begin{equation}
\tag{4}
\frac{\partial SSQ}{\partial K_m} = 0 = -2\sum_{i=1}^n[v_i - V_{max} \frac{[S]_i}{K_m + [S]_i}] \frac{V_{max}[S]_i}{[S]_i + K_m}
\end{equation}\]
To achieve such purpose, the function dir.MM() invokes numerical method well implemented into R.
Loading some real kinetic data
We start by loading some kinetic data obtained by students during their undergraduate laboratory training. Using \(\beta\)-galactosidase as an enzyme model, the students assess the effect of the substrate o-nitrophenyl-\(\beta\)-D-galactopynaroside (ONPG) on the initial rate (doi: 10.1002/bmb.21522). The data obtained by eight different groups of students can be loaded just typing:
data <- ONPG
library(knitr)
kable(data)
| 0.05 |
2.26 |
1.29 |
0.004 |
0.004 |
0.004 |
0.003 |
1.77 |
2.98 |
| 0.10 |
5.48 |
3.33 |
0.008 |
0.007 |
0.007 |
0.006 |
5.20 |
5.20 |
| 0.25 |
13.40 |
11.80 |
0.020 |
0.020 |
0.016 |
0.017 |
15.04 |
14.38 |
| 0.50 |
24.70 |
22.80 |
0.035 |
0.035 |
0.032 |
0.031 |
28.31 |
30.30 |
| 1.00 |
40.90 |
35.20 |
0.060 |
0.056 |
0.050 |
0.048 |
50.98 |
48.99 |
| 2.50 |
62.30 |
39.90 |
0.110 |
0.104 |
0.090 |
0.101 |
75.42 |
86.25 |
| 5.00 |
94.30 |
73.50 |
0.138 |
0.138 |
0.115 |
0.121 |
112.68 |
112.57 |
| 8.00 |
105.00 |
12.90 |
0.154 |
0.150 |
0.119 |
0.139 |
126.06 |
136.24 |
| 20.00 |
133.00 |
112.00 |
0.179 |
0.179 |
0.142 |
0.152 |
154.93 |
169.97 |
| 30.00 |
144.00 |
120.00 |
0.200 |
0.200 |
0.166 |
0.181 |
168.75 |
177.71 |
The first column gives the ONPG concentrations in mM, and the remaining 8 columns correspond to the initial rates. Note that while groups 1, 2, 7 and 8 decided to express their rates as \(\mu\)M/min, the remaining groups opted by mM/min. This information can be confirmed by checking the attributes of data:
attributes(data)
#> $names
#> [1] "ONPG" "v1" "v2" "v3" "v4" "v5" "v6" "v7" "v8"
#>
#> $row.names
#> [1] 1 2 3 4 5 6 7 8 9 10
#>
#> $`[ONPG]`
#> [1] "mM"
#>
#> $`v3, v4, v5, v6`
#> [1] "mM/min"
#>
#> $class
#> [1] "data.frame"
#>
#> $`v1, v2, v7, v8`
#> [1] "uM/min"
Thus, before continuing we are going to express all the rates using the same units: \(\mu\)M/min:
data[ , 4:7] <- 1000 * data[ , 4:7]
First thing first: scatter plot
I strongly insist to my students that when we have to analyze data, the first thing we must do is a scatter diagram, since this will give us a first impression about our data and will guide us on how to proceed with the analysis. To lead by example, we will carry out such diagrams.
The first four groups:
oldmar <- par()$mar
oldmfrow <- par()$mfrow
par(mfrow = c(2, 2))
par(mar = c(4, 4,1,1))
for (i in 2:5){
plot(data$ONPG, data[, i],
ty = 'p', ylab = 'v (uM/min)', xlab = '[ONPG] (mM)')
}
par(mar = oldmar)
par(mfrow = oldmfrow)
The next four groups:
oldmar <- par()$mar
oldmfrow <- par()$mfrow
par(mfrow = c(2, 2))
par(mar = c(4, 4,1,1))
for (i in 6:9){
plot(data$ONPG, data[, i],
ty = 'p', ylab = 'v (uM/min)', xlab = '[ONPG] (mM)')
}
par(mar = oldmar)
par(mfrow = oldmfrow)
In general, the data does not provide us with any surprises. That is, the relationship between the dependent variable (initial rate) and the independent variable ([ONPG]) is what we expect: hyperbolic curve. An exception is the rate obtained by group 2 when [ONPG] = 8 mM, which is clearly an “outlier”. No problem! We will remove that point from further analysis to prevent it from introducing artifacts.
Direct Michaelis-Menten fitting
Using the data from group 7 to illustrate the use of the dir.MM() function:
dir.MM(data[ , c(1,8)], unit_v = "mM/min")
#> $parameters
#> Km Vm
#> 3.115 181.182
#>
#> $data
#> S v fitted_v
#> 1 0.05 1.77 2.862275
#> 2 0.10 5.20 5.635521
#> 3 0.25 15.04 13.460773
#> 4 0.50 28.31 25.059751
#> 5 1.00 50.98 44.029648
#> 6 2.50 75.42 80.668744
#> 7 5.00 112.68 111.634011
#> 8 8.00 126.06 130.405398
#> 9 20.00 154.93 156.765737
#> 10 30.00 168.75 164.138910
We propose to the reader, as an exercise, to compare these results with those obtained when we use data from the remaining groups.