R gsl solver interface
This R package solves a series of initial value problems given as an ordinary differential equation, an initial state and a numerical parameter vector. The C code calls gsl odeiv2 module functions to solve the problem or problems.
The parameter vector can be replaced by an n×m matrix of several parameterisations (each column is a distinct parameter vector). The solver wil be called m times.
The result is returned to R as a 3-dimensional array.
This project was funded by the EU as part of the Human Brain Project and supported by EBRAINS infrastructure.
Using the remotes package:
remotes::install_github("a-kramer/rgsl")
Ensure that the GNU Scientific Library (gsl) is installed in your system and pkg-config can find it:
$ pkg-config --libs gsl
The above command should print something similar to:
-lgsl -lgslcblas -lm
In this package we consider the initial value problem:
ẋ = f(x, t; p)x(t0) = x0
where p is a parameter vector.
We think that it is helpful to distinguish between constants (parameters that never change), unknown parameters, and input parameters.
| Type | Symbol | Description |
|---|---|---|
| constants | NONE | never change, are built into the model |
| unknown | k | may be inferred from data |
| known | u | these distinguish different experiments |
The overall vector p contains both the unknown k and the known
u, e.g. as p=c(k,u) (aka a direct sum). But any ordering is possible (here c means
concatenate, like in R).
To determine the parameters k, the model f needs to be simulated in the provided experimental scenarios and the parameters are found through optimisation, or sampling (or perhaps another, similar method).
Typically, these methods need to solve the model many times. This is especially true for Bayesian methods (sampling). So, solving a batch of similar problems in one function call is desirable. In the next two Sections we discuss two common Scenarios.
Let us assume that we have obtained a data set D, that includes measured observables form different experiments. Each experiment is distinguished from the others by the values in the input to the model.
The input may well be dynamic, and if it is, the model needs to contain the possible input dynamics within the model files. But, let us assume here that we can switch between input dynamics by setting specific input parameter values: u₁, u₂, u₃, …. We shall evaluate the model for a given suggested parameter vector k. Then we can construct a p for each input: p₁=c(k,u₁). Example R code:
# two experiments, "Stim" and "Control"
u<-load_inputs_from_file("my_file_of_inputs.dat")
k<-c(1,0,2,3) # let's say this makes sense
p<-matrix(c(k,u$Stim,k,u$Control),ncol=2)
t=seq(0,10,length.out=12)
y0=c(1,0,0) # initial state vector
library(rgsl)
y<-r_gsl_odeiv2("NameOfModel",t,y0,p)
# the below plot will correspont to parameters c(k,u$Stim)
plot(t,y[1,,1])
# and this will be the ratio between Stim and Control:
lines(t,y[1,,1]/y[1,,2])
The result is a three dimensional array y, where y[i,j,l]
will be the state variable i at time j for parameter vector l.
Similarly, if an MCMC sample K for a given model already exists, the user may want to simulate a prediction based on this sample. To do this the whole sample may be provided to the C solver used here and the result will be a batch of trajectories, one for each sampled parameter vector. The sample matrix K is a set of column-vectors, each a valid parameterisation of the model. We want to predict the system’s behaviour for a given input of interest: u.
We construct a parameter matrix:
K<-load_sample_from_file("sample_file.hdf5")
u<-c(1,2,3)
n<-dim(K)
P<-rbind(K,matrix(u,nrow=3,ncol=n(2)))
So, the input is now repeated for each column in K. Then we call the
r_gsl_odeiv2() interface function to simulate the scenario u for
each parameter vector. The result is a matrix y[,,l] where l can
be any number in the range 1:n(2).
This package requires the model to be provided as a shared library:
${MODEL}.so. The function names within the file need to correspond to the model name:
${MODEL}_vf();
${MODEL}_jac();
See the example file HarmonicOscillator to inspect an auto-generated right-hand-side and Jacobian function compatible with the solvers from the gsl.
The example script HarmonicOscillator.R
contains a demo() function.
These C code files (for the gsl module odeiv2) can be automatically obtained using from VFGEN or written manually.
This project integrates well with SBtabVFGEN as this package prints vfgen files from source files written using the SBtab format. SBtab is a format in systems biology and is suited for writing/drafting biological models. The workflow could be:
User written: generated Simulation
+-----------+ +-----------+ +------------+ +----------+
| | | | | (CVODE) | | |
| SBtab (M) +--+-->+ VFGEN +----->+ ODE code +--+-->+ rgsl |
| | | | | | +jacobian | | | |
+-----------+ | +-----------+ +------------+ | +----------+
| |
+----+--------------+ +----------------+
| | | |
| sbtab_to_vfgen() | | r_gsl_odeiv2() |
| | | |
+-------------------+ +----------------+