In a previous post I used R to solve **O**rdinary **D**ifferential **E**quations (ODE). Since I'm still more familiar with Matlab to work with ODE, I decided to go one step further and learn how to use R to estimate parameters in ODE.

In this short tutorial, I will use the ODE presented here which quantifies the salt concentration in a well-stirred tank:

\[ \begin{align} &\frac{dx}{dt}=a\frac{150-x(t)}{200}\\[10pt] \end{align} \]

With \[ x(0) = 20 \]

The analytic solution is:

\[ \begin{align} &x(t)=150-130e^{-at/200} \end{align} \]

**Step 1**

Generate data using the analytic solution and add some random noise. Note that the ODE uses one parameter that has been fixed at \( a = 0.75 \).

```
t = seq(1, 800, by = 10)
a = 0.75 ## Fixed parameter used to simulate data.
## Simulate data and add noise.
conc.observed = 150 - 130 * exp(-a * t/200) + rnorm(length(t), sd = 5)
## Plot
plot(t, conc.observed, pch = 21, bg = "gray", ylim = c(20, 180), xlab = "Time",
ylab = "Salt (kg)", axes = F)
box(bty = "l")
axis(1, tck = -0.02)
axis(2, tck = -0.02, las = 2)
legend("bottomright", legend = c("Observed data"), pch = 21, pt.bg = "gray",
bty = "n")
```

**Step 2**

Write a function that will be used to solve the ODE. While we there, solve it with fixed parameter \( a = 0.75 \).

```
library(deSolve)
salttank = function(t, state, parameters) {
with(as.list(c(state, parameters)), {
# rate of change
dX = a * ((150 - X)/200)
# return the rate of change
list(c(dX))
})
}
## Solve the ODE. Not necessary at this point.
## Define the initial value for the state variable
state = c(X = 20)
## Time range.
times = seq(1, 800, by = 10)
## Parameters
parameters = c(a = 0.75)
conc.modeled = ode(y = state, times = times, func = salttank, parms = parameters,
method = "rk4")
```

Plot the results.

```
plot(t, conc.observed, pch = 21, bg = "gray", ylim = c(20, 180), xlab = "Time",
ylab = "Salt (kg)", axes = F)
box(bty = "l")
axis(1, tck = -0.02)
axis(2, tck = -0.02, las = 2)
lines(conc.modeled[, "time"], conc.modeled[, "X"], col = "red")
legend("bottomright", legend = c("Observed data", "True solution (a = 0.75)"),
lty = c(NA, 1), pch = c(21, NA), col = c("black", "red"), pt.bg = "gray",
bty = "n")
```

**Step 3**

This is where the magic happens. To fit parameters, I will use **nls.lm()** from the **minpack.lm** package. In its simplest form, the function uses *par* which is a list of guessed parameters and *fn* a function used to minimize the sum square of the residuals using the Levenberg-Marquardt algorithm.

Now we have to write a function (here I named it *ssq*) that will use ODE parameters as input (in this case only \( a \)) and produces a residuals vector as output.

```
ssq = function(params) {
## Range and initial condition as before.
times = seq(1, 800, by = 10)
state = c(X = 20)
## Resolve the ODE.
out = ode(y = state, times = times, func = salttank, parms = params, method = "rk4")
## modeled - observed
ssq = out[, "X"] - conc.observed
}
```

**Step 4**

Finally, we can estimate the parameter \( a \).

```
library(minpack.lm)
## Start with a guess a = 1
params.guessed = c(a = 1)
params.fitted = nls.lm(par = params.guessed, fn = ssq)
summary(params.fitted)
```

```
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 0.7438 0.0113 66 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.93 on 79 degrees of freedom
## Number of iterations to termination: 5
## Reason for termination: Relative error in the sum of squares is at most `ftol'.
```

We see from the summary of *params.fitted* that the estimated parameter is \( a = 0.7438 \) which is obviously very close to \( a=0.75 \) used to produce the observed data.

**Step 5**

Finally, simulate the data using the estimated parameter and plot the results.

```
## Simulate using fitted parameters
params = coef(params.fitted)
conc.modeled3 = ode(y = state, times = times, func = salttank, parms = params,
method = "rk4")
```

```
plot(t, conc.observed, pch = 21, bg = "gray", ylim = c(20, 180), xlab = "Time",
ylab = "Salt (kg)", axes = F)
box(bty = "l")
axis(1, tck = -0.02)
axis(2, tck = -0.02, las = 2)
lines(conc.modeled[, "time"], conc.modeled[, "X"], col = "red")
lines(conc.modeled3[, "time"], conc.modeled3[, "X"], col = "blue")
legend("bottomright", legend = c("Observed data", "True solution (a = 0.75)",
paste("Simulated with fitted parameters (a = ", sprintf("%2.2f", params),
")", sep = "")), lty = c(NA, 1, 1), pch = c(21, NA, NA), col = c("black",
"red", "blue"), pt.bg = "gray", bty = "n")
```

Super,

ReplyDeleteJe pensais que tu aurais fait cette modÃ©lisation pour la croissance d'une population de levure de biÃ¨re! Peut Ãªtre pour une prochaine.

Sinon quand est-ce que tu t'attaque Ã PARAFAC?

Bravo pour le travail de ton blog!

Thomas

Merci Thomas!

ReplyDeleteL'idÃ©e de la biÃ¨re serait trÃ¨s intÃ©ressante. D'ailleurs, cela ferait un excellent sujet pour un cours Ã la maÃ®trise avec les Ã©tudiants. Suivre la cinÃ©tique des levures dans une expÃ©rience de fermentation.

Bon courage,

Phil

Many thanks for making the truthful effort to explain this. I feel very strong about it and would like to read more. If you

ReplyDeletecan, as you find out more in depth knowledge, would you mind posting more posts similar to this one with more information.

Qassim University

.,

Thank you for your comment. I have some posts coming soon with more detailed examples.

DeleteThank you for the example. I am also more familiar with matlab, maple and mathematica. Do you know how to fit also the initial condition ? i.e state = c(X = 20), intead of 20 a new parameter X0 ? Thanks.

ReplyDelete