Introduction
The algorithms contained in the seeds package are implementation of the dynamic elastic net approach. There are two algorithms implemented that aim to calculate missing interactions between states of an ode in order to explain and minimize a systematic error within the model-equation. The hidden inputs should be an hint as to which states are missing interactions.
Both algorithm makes heavy use of the deSolve package. The use of Rtools is recommended to significantly shorten the computing time of the algorithm. It is not possible to use the BDEN algorithm without Rtools.
The DEN (dynamic elastic net) Algorithm
The first Algorithm is a greedy method for generating a sparse solution for the hidden inputs. The algorithm is based on Pontryagin’s Maximum Principle and implements an iterative gradient descent method which aims at minimizing the Hamilton function of the system.
Rtools is suggested, because of the immense computational cost of the step size estimation of the gradient step.
The odeModel Class
The class aims to be used to describe the experiment you want to analyze. It has the following slots:
- func: A function containing the ode-equations of the model. For syntax look at the given examples of the deSolve package.
- times: time steps at which the model should be evaluated
- parms: the parameters of the model
- input: matrix containing the inputs with the time points
- measFunc: function that converts the output of the ode solution
- y: initial (state) values of the ODE system, has to be a vector
- meas: matrix with the (experimental) measurements of the system
- sd: optional standard deviations of the measurements, is used by the algorithms as weights in the cost function
- custom: customized link function (is internally used)
- nnStates: bit vector that indicates if states should be observed by the root function
- nnTollerance: tolerance at which a function is seen as zero (optional)
- resetValue: value a state should be set to by an event (optional)
The class also provides methods to simply evaluate and display the nominal solution of an created odeModel object with the nominalSol and plot() function.
Example
The following example uses the UVB protein g model and some simulated measurement data. An synthetic error has been added to the states \(x_3\) and \(x_{11}\).
The states of the model have to be written with the letter \(x\) following by a number representing the state. If an input is needed for the model the input takes the letter \(u\). This restrictions are made to ensure that the automatic differentiation will work as intended.
library(seeds)
# usb network example
uvbParameter = c( ks1=0.23,
ks2=4.0526,
kdr1=0.1,
kdr2=0.2118,
k1=0.0043,
k2=161.62,
ka1=0.0372,
ka2=0.0611,
ka3=4.7207,
kd1=94.3524,
kd2=50.6973,
kd3=0.5508,
ks3=0.4397,
kdr3=1.246,
uv=1,
ka4=10.1285,
kd4=1.1999,
n1=3,
n2=2,
n3=3.5,
kdr3a=0.9735,
kdr3b=0.406,
ksr=0.7537,
fhy3_s=5)
# initial state
x0 = c(0.2,10,2,0,0,20,0,0,0,4.2,0.25,20,0)
# model ode
uvbModel <- function(t,x,parameters) {
with (as.list(parameters),{
dx1 = ((-2) * ((ka1 * (x[1]^2) * (x[4]^2)) - (kd1 * x[5])) +
(-2) * ((ka2 * (x[1]^2) * x[2]) - (kd2 * x[3])) +
((ks1 *((1) + (uv * n3 * (x[11] + fhy3_s)))) -
(kdr1 * ((1) + (n1 * uv)) * x[1])))
dx2 = ((-1) * ((ka2*(x[1]^2) * x[2]) - (kd2 * x[3])) +
(-1) * ((ka4 * x[2] * x[12]) - (kd4 * x[13])))
dx3 = (((ka2 * (x[1]^2) * x[2]) - (kd2* x[3])))
dx4 = ((-2) * (k1*(x[4]^2)) + (2) * (k2 * x[6]) +
(-2) * ((ka1 * (x[1]^2)* (x[4]^2)) - (kd1 * x[5])) +
(-1)* (ka3 * x[4] *x[7]))
dx5 = (((ka1 * (x[1]^2) * (x[4]^2)) -(kd1 * x[5])))
dx6 = ((-1) * (k2 * x[6]) + (k1 * (x[4]^2)) +(kd3 * (x[8]^2)))
dx7 = ((-1) * (ka3 * x[4] * x[7]) + ((ks2 * ((1) + (uv * x[5]))) -
(kdr2 * x[7])) + (2) * (kd3 * (x[8]^2)))
dx8 = ((-2) * (kd3 * x[8]^2) + (ka3 * x[4] * x[7]))
dx9 = 0
dx10 = 0
dx11 = (((ks3 * ((1) + (n2 * uv))) -(kdr3 * (((x[3] / (kdr3a + x[3])) +
(x[13] / (kdr3b + x[13]))) -(x[5] / (ksr + x[5]))) * x[11])))
dx12 = ((-1) * (ka4 * x[2] * x[12]) + (kd4 * x[13]))
dx13 =((ka4 * x[2] * x[12]) - (kd4 * x[13]))
list(c(dx1,dx2,dx3,dx4,dx5,dx6,dx7,dx8,dx9,dx10,dx11,dx12,dx13))
})
}
# measurement function
uvbMeasure <- function(x) {
y1 = 2*x[,5] + x[,4] + x[,8]
y2 = 2*x[,5] + 2* x[,3] + x[,1]
y3 = x[,6]
y4 = x[,11]
y5 = x[,4]
return(cbind(y1,y2,y3,y4,y5))
}
# measurement data
y <- uvbData[,1:6]
# timesteps
t <- uvbData$t
# standard deviation of the measurements
sd <- uvbData[,7:11]
uvbModel <- odeModel(func = uvbModel, parms = uvbParameter, times = t,
measFunc = uvbMeasure, y = x0, meas = y, sd = sd)
res <- DEN(odeModel = uvbModel, alphaStep = 500, alpha2 = 0.0001, epsilon = 0.2, plotEstimates = FALSE)The results are returned as list of resultsSeeds objects, which store the estimated states, measurements and the hidden inputs of a model. A standard plot function is given to display the results. Functions for annotating the plots are also included in the package. The use of the BDEN algorithm displays not only the calculated states and inputs, but also confidence intervals for the corresponding states of the system.
Plotting
The standard plot gives you three figures:
- The estimated states (blue without hidden inputs, red with hidden inputs)
- The estimated hidden inputs
- The estimated measurements. (blue without hidden inputs, red with hidden inputs)
All calculated values can be output by methods of the class odeModel, or extracted directly by the corresponding slots. The structure and methods of the class can be found in the corresponding documentation.
The plotting function will give you three plots.
Estimated States for the model
Estimated hidden inputs
Measurements resulting from the estimated hidden inputs