Team:TU-Munich/Modeling/Enzyme

From 2013.igem.org

Revision as of 08:16, 25 October 2013 by ChristopherW (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

Introduction

To estimate the efficiency of our remediation rafts in detoxification of contaminated water, we saw the need to assess the activity of our target effectors. As in vivo experiments were still to be carried out we, started to analyse data we obtained in in vitro experiments.

We analysed the data for

**Table 1:** Experimental Conditions
Experiment
[#Laccase (Substrate Dependence) Laccase (Substrate Dependence)]
Erethromycin esterase B (Substrate Dependence)
Laccase (pH Dependence)
Erethromycin esterase B (pH Dependence)
Laccase (NaCl Dependence)
Erethromycin esterase B (NaCl Dependence)

Laccase (Substrate Dependence)

Data

To assess the activity of <partinfo>BBa_K1159002</partinfo> we created an experimental setup with a dilution series for the concentration of the substrate ABTS with 16 different concentrations with a 2:1 dilution ratio resulting in concentrations ranging from 2.429 mM to 74.1 nM in triplicates and added 40.8 nM of purified enzyme. The concentration of the product of the enzymatic reaction was measured via absorption over 5 minutes with time steps of 10 seconds. As absorption measurements for very low substrate concentrations were almost undistinguishable from noise, we decided only to fit the data for the 8 highest concentrations. The data for the 8 highest concentration including standard deviation is shown in Figure 1.

Figure 1: Experimental data with error bars plotted at 1-σ interval.

Model

Calibration Curve

For quantitative modelling it is essential to relate the measured absorption with the concentration of the product of the enzymatic reaction. After initial investigation of the data we realised that the ratio between the absorption in the steady states was not scaling linearly with the dilution ratios. This led us to the conclusion that there was a nonlinear relationship between the concentration of the product of the enzymatic reaction and the absorption. Hence we contacted the wetlab and they measured a calibration curve for us such that we could carry out rigorous quantitative modeling.

Figure 2: Calibration data indicated by error bars displaying the 1-σ interval in blue and the fitted logarithmic calibration line in red.

To use this curve later for modelling we fitted a parametric logarithmic curve in accordance with the [http://en.wikipedia.org/wiki/Beer–Lambert_law Beer-Lambert law]

by minimizing the corresponding log-likelihood assuming a normally distributed error. The obtained data and the corresponding fitted logarithmic curve are shown in figure 2, the resulting best fitting parameters are shown in table 1.

**Table 1:** Calibration fitting
Parameter	Best Fit
k₁	1.651 a.u.
k₂	3891 1/M

Model Specification

For the Laccase kinetics we used the standard model for enzyme kinetics in reduced form.

We assumed a normally distributed measurement error which yielded us the following negative log likelihood

where is the mean of the experimental measurements at time point t_i for the j-th experiment in the dilution series, σ_ij denotes the corresponding standard deviation and y_j is the numerical simulation of system Θ₁ for parameters θ = [k_f, k_r, k_cat] with the calibration curve applied to the output P at time t_i.

Fitting

To obtain the best fitting parameters we used a multistart optimization in logarithmic parameter space with latin hypercube sampling as initialisation. Multistart optimization is employed to be sure that the resulting best fit is not a local optimum but a global one and logarithmisation of the parameter space is well-known to yield better results http://www.plosone.org/article/info:doi/10.1371/journal.pone.0074335 Raue et al., 2013. The parameter boundaries for the hypercube are shown in table 2.

**Table 2:** Hypercube parameters
Parameter	Lower Bound	Upper Bound
k_f	1E+03 1/s	1E+07 1/s
k_r	1E-05 1/s	1E+06 1/s
k_cat	1E+00 1/s	1E+08 1/s

We used 20 initialisations for the multistart. Figure 4 shows the final value for the optimization routine for the different initial values, sorted by the resulting final value. As one can see there are 11 initialisations that yield approximately the same final result, which suggests that this might be a global optimum. Figure 3 displays the obtained parameters of the best fits colored according to their likelihood, where black is more likely and white is less likely.

Figure 3: Fitted parameter values

Figure 4: Final Negative log-likelihood of different optimization initialisations

These results suggest that parameter 3 (k_cat) is very well identifiable, as all optimization runs yielded approximately the same parameter value. In contrast, the estimation of the two other parameters show strong uncertainties as the different initialisations yielded a range of different parameter values. Using the obtained parameters we can also compute K_M which is used when Michaelis-Menten kinetics are assumed. The formula for this is

The obtained parameters for the best fit are shown in Table 3:

**Table 3:** Best fit parameters
Parameter	Best fit
k_f	4.196 E+06 1/s
k_r	1.399 E+03 1/s
k_cat	61.892 1/s
K_M	33.447 μM

The simulation of the resulting parameter compared against the experimental data is shown in Figure 5.

Figure 5: Fitted data

These results suggest that we have obtained a reliable estimate for some of the parameter values and that the model seems sufficient to explain the experimental data. As there seems to be a certain uncertainty with respect to some of the parameters, we decided to carry out further uncertainty analysis to asses the reliability of the parameter estimates.

Uncertainty Analysis

The uncertainties in measurement can be directly transferred into uncertainties in the parameter estimation with the help of the likelihood function. As the computation of the likelihood includes solving a non-trivial differential equation it is not possible to write down an analytical expression for it that could be used to estimate the parameter distributions. Hence we will have to resort to sampling from the parameter distribution. We did this using the Metropolis-Hastings algorithm implemented in the adaptive MCMC toolbox DRAM http://helios.fmi.fi/~lainema/mcmc/ Haario et al. (2006).We used a burn-in of 1,000 samples and then generated 400,000 samples from the distribution. Using these samples we can also generate samples for the parameter K_M by computing K_M for every sample via the same formula that was used above. For visualization 500 samples out of these 200,000 are shown in Figure 6. Using these samples we used standard kernel density estimation to estimate the marginal densities of single parameters. The resulting approximate marginal distributions are shown in Figure 7.

Figure 6: Thinned samples and their respective parameter values

Figure 7: Marginal densities of the individual parameters obtained by using kernel density estimation on the MCMC kernels. The best fit parameters are indicated by red dotted lines.

In Figure 6 and 7 one can see that all parameters except k_r have localized parameter ranges of high probability, which means that we can expect our optimization should yield good results. Figure 7 shows the marginal densities of the parameters, obtained via kernel density estimation using MCMC samples, with the best parameter fits, obtained via optimization, indicated by red lines. As one can see, all fitted parameter except for k_f are located in regions of high probability density.

We can explain this phenomenom by looking at the correlation between parameters. The correlation matrix of the parameter samples obtained via MCMC is presented in Table 4:

**Table 4:** Correlation Matrix
	k_f	k_r	k_cat
k_f	1	0.698	-0.054
k_r	0.698	1	0.027
k_cat	-0.054	0.027	1

Investigating the values presented in Table 4 one can see that there is a high correlation between parameters k_f and k_cat. The joint distribution of the samples in their values for k_f and k_r is shown in Figure 8. The likelihood of the individual samples is indicated by color. The best fit parameters are again displayed via a red dotted line. One can clearly see that there is a strong nonlinear correlation between the two parameters. Moreover it is evident, that despite having a low marginal likelihood for the obtained value of k_f, in the joint distribution the best fit lies in a strongly localised area and we can accept the results with high certainty.

Figure 8: Correlation between k_f and k_r

Erythromycin Esterase B (Substrate Dependence)

Data

To assess the activity of <partinfo>BBa_K1159000</partinfo> we created an experimental setup with a dilution series for the concentration of the substrate p-NPB with 16 different concentrations with 2:1 dilution ratio resulting in concentrations ranging from 1.195 mM to 36.469 nM in triplicates and added 5.804 μM of purified enzyme. The concentration of the product of the enzymatic reaction was measured via absorption for the duration of 2 hours with a time steps of 30 seconds. As measurements for very low substrate concentrations were almost undistinguishable from noise, we decided only to fit the data for the 6 highest concentrations.

Figure 9: Experimental Data with error bars plotted at 1-σ interval.