Dynamical Systems in Biology

Exercise 1: Non-linear regression

The following table presents data for the hemoglobin saturation dependin on $O_2$ pressure:

mydata <- data.frame(
  O2 = c(3.08, 4.61, 6.77, 10.15, 12.31, 15.38, 18.77, 22.77, 25.85, 30.15, 36.00, 45.23, 51.69, 61.85, 75.38, 87.08, 110.5),
  Y = c(2.21, 3.59, 6.08, 10.50, 14.09, 19.34, 28.45, 40.33, 50.00, 60.50, 69.89, 80.11, 83.98, 88.95, 93.37, 95.86, 98.07)
)

knitr::kable(mydata, align = c("c", "c"))

Fit these data to a curve of the form

$$\begin{align} Y = 100\cdot\frac{\sum_{j = 0}^4 j\alpha_j [O_2]^j}{4\sum_{j = 0}^4 \alpha_j [O_2]^j}. \end{align}$$ Instruction: Suppose we have data points $\{x_i, y_i\}$, $i = 1, \dots, n$, that we wish to fit to some function $y = f(x)$, and that the function $f$ depends on parameters $\{\alpha_i\}$, $i = 1, \dots, m$. A fit of the data is achieved when the parameters are picked such that the function $$\begin{align} F = \sum_{i = 1}^n (f(x_i) - y_i)^2 \end{align}$$ is minimized. Fo find this fit, start with reasonable estimates for the parameters and then allow them to change dynamically (as a function of a time-like variable $t$) according to $$\begin{align} \forall j:\quad \frac{d\alpha_j}{dt} = - \sum_{i = 1}^n \big(f(x_i)-y_i\big)\frac{\partial f(x_i)}{\partial \alpha_j}.\qquad\textrm{(1)} \end{align}$$ With this choice, $$\begin{align} \frac{dF}{dt} = \sum_{j = 1}^{m}\sum_{i = 1}^n 2\big(f(x_i)-y_i\big)\frac{\partial f(x_i)}{\partial \alpha_j}\frac{d\alpha_j}{dt}\leq 0, \end{align}$$ so that $F$ is a decreasing function of $t$. A fit is found when numerical integration of $\frac{d\vec\alpha}{dt}$ reaches a steady-state solution.

Implementation:

• Define a function model(x, pars) that computes the model prediction $Y$ from oxygen levels $[O_2]$, denoted by x and values $\alpha_j$, denoted by pars. Compute derivatives $\frac{\partial Y_i}{\partial \alpha_j}$ within this function and return the matrix of derivatives together with the model prediction vector, e.g. as an attribute or in a list. Advise: Write a small function to check whether the analytically computed derivates are correct.

## Set-up model with derivatives
model <- function(x, pars) {
  j <- 0:4
  predictor <- outer(x, j, function(a, b) a^b)
  jpredictor <- outer(x, j, function(a, b) b*a^b)
  
  value <- 100*as.vector((jpredictor %*% pars)/(4*predictor %*% pars))
  gradient <- 100*(jpredictor/as.vector(4*predictor%*%pars)- 
    4*predictor*as.vector(predictor %*% (j * pars))/as.vector(4*predictor %*% pars)^2)
  
  attr(value, "gradient") <- gradient
  
  return(value)
  
}

## Check if derivatives are implemented correctly
pars <- rnorm(5)
x <- rnorm(3)
grad.analytical <- attr(model(x, pars), "gradient")
grad.numerical <- numDeriv::jacobian(function(pars) model(x, pars), pars)

precision <- max(abs(grad.analytical - grad.numerical))
rel.precision <- precision/mean(abs(grad.numerical))

cat("precision:", precision, ", rel. precision:", rel.precision, "\n")

## precision: 1.031212e-07 , rel. precision: 3.695935e-11

• Define a function func(times, states, parameters, mydata) for the dynamic system, eq. (1), to be used with ode(). Integrate the system and show the integration result on log-time scale. Check, if the integration time was sufficient to reach the steady-state.

## Set-up dynamic model for integration
func <- function(times, states, parameters, mydata) {
  prediction <- model(mydata$O2, states)
  dalpha <- -(prediction - mydata$Y) %*% attr(prediction, "gradient")
  list(dalpha)
}

alpha <- exp(rnorm(5, -1))
times <- c(0, exp(seq(-22, 0, .1)))
out <- ode(alpha, times, func, mydata = mydata)


estimate <- tail(out, 1)[,-1]
cat("estimated parameters:", estimate, "\n")

## estimated parameters: 1.365444 0.03490661 0.000829054 2.69911e-05 3.30215e-06

library(tidyr)
out <- gather(as.data.frame(out), parameter, value, -time)
ggplot(out, aes(x = time, y = value)) + facet_wrap(~parameter, scales = "free") + 
  geom_line() + scale_x_log10(labels = scales::trans_format("log10", scales::math_format(10^.x)))

## Warning: Transformation introduced infinite values in continuous x-axis

• Generate a plot showing the data points, the model prediction prior to optimization and after optimization.

x <- seq(min(mydata$O2), max(mydata$O2), 1)
prediction <- data.frame(O2 = x, Y = model(x, estimate))
prediction.ini <- data.frame(O2 = x, Y = model(x, alpha))

ggplot(mydata, aes(x = O2, y = Y)) + geom_point() + 
  geom_line(data = prediction) +
  geom_line(data = prediction.ini, lty = 2, color = "red")

• What is the maximal degree $j_{\max}$ needed to fit the data accurately? What is the measurement uncertainty of the data?

Cathedral exercise:

Episcopal churches have two towers. Freiburg is a diocesan town. Why does the cathedral have only one tower?