Dynamical Systems in Biology

Tutorial 10

Exercise 1: The profile-likelihood method

The profile likelihood method is used to assess confidence intervals of parameters in non-linear models. A model of the reaction flux of an enzymatic reaction associated with the substrate concentration \([S]\) is given by \[ \begin{align*} \dot v = \frac{10^{\rm Vmax}\cdot [S]}{10^{\rm Km} + [S]}. \end{align*} \] Here, Km and Vmax denote the Michaelis-Menten constant and maximal flux on the \(\log_{10}\) scale, respectively. Let us assume that the reaction flux \(\dot v\) can be measured as a function of the substrate concentration and that the parameters Vmax and Km are to be determined from the experimental data.

1. Simulate the dose-response curve \(\dot v([S])\) for parameters Km = 0 and Vmax = 1 and concentrations of \([S]\) between 0 and 2. Plot the result.

2. For concentrations \([S]\in\{.1, .5, 1, 2\}\) simulate the model prediction and pick normally distributed random numbers around the prediction with a standard deviation of 7%. These values represent the simulated data with measurement uncertainty. Add the data points with error bars to the previous plot.

3. The basis of the profile-likelihood method is the log-likelihood function and a maximum-likelihood fit serving as the point of departure for the profile likelihood analysis. Define the negative log-likelihood as a function objfn <- function(Km, Vmax) of the parameters. Make explicit use of the previously defined data.frame with simulated data and its uncertainty values. You can either analytically derive the expression for the negative log-likelihood from the Gauß function or use the dnorm() function. Finally, use mle2() from the bbmle package to minimize the negative log-likelihood and obtain parameter estimates.

4. Implement an algorithm for the profile-likelihood method. In principle, the algorithm should include:

• a while-loop to compute the deviance between the original LL-value and the new LL-value after refitting for the stepwise increased parameter value.
• a for-loop to include both directions (sign) for the profile computation starting from the optimum
• a for-loop over the parameters whichPar

The algorithm should have the general form:

For whichPar in parameters
  
  For sign in (-1, 1)
    
    Set mypars to parameters from original fit
    Set deviance to 0
  
    While deviance < 4 and parameter values in [-5, 5]
  
      fit the model with value of whichPar fixed to value(whichPar) + sign*stepsize
      overwrite mypars by fit result
      compute deviance = -2 ( logLik(fit_0) - logLik(fit_new) )
  
    End
  
  End    
      
End

Generate a data.frame with columns deviance (the difference of twice the negative log-likelihood between the “refit” and the original fit), fixed (the value of the parameter held fixed during max-likelihood estimation), Km, Vmax (the fit result per iteration of the while-loop) and parameter (whichPar, the name of the parameter for which the profile was computed). Show the results in (1) a plot showing the fixed parameter value vs. the deviance and (2) a plot showing Km vs Vmax for each of the profiles.

5. For the above choice of substrate concentrations we find that the profiles of both parameters exceed the deviance threshold associated with a 95% confidence level. Repeat steps 1-4 for simulated data points with \([S]\in(0, .1]\) or with \([S]\in[5, 10]\). Interpret the results obtained for the likelihood-profiles and profile paths.

Cathedral exercise:

How does the cathedral relate to the symposium of Plato? Study the Apocalypse of Saint John and the writings of Proklos, Plotinus and the Abbot Suger of St. Denis.