It looks like MATLAB, Octave and Python seem to be the preferred tools for scientific and engineering analysis (especially those involving physical models with differential equations). However as part of my learning R experience, I wanted to check out some of R tools for parameter fitting of models involving ordinary differential equations. R has packages deSolve for solving differential equations and FME for parameter fitting. The specific example here is taken from the computational appendix (A.6) of the book Chemical Reactor Analysis and Design Fundamentals by Rawlings and Ekerdt. In fact, all examples in this book are available in Octave and MATLAB.

This example involves two reactions in series

where we need to estimate the rate constants of the the two reactions and from data

The rate equations are captured in a function that is an input parameter to the ODE solver

# rate function
rxnrate=function(t,c,parms){
# rate constant passed through a list called parms
k1=parms$k1
k2=parms$k2
# c is the concentration of species
# derivatives dc/dt are computed below
r=rep(0,length(c))
r[1]=-k1*c["A"] #dcA/dt
r[2]=k1*c["A"]-k2*c["B"] #dcB/dt
r[3]=k2*c["B"] #dcC/dt
# the computed derivatives are returned as a list
# order of derivatives needs to be the same as the order of species in c
return(list(r))
}

Computing the predicted concentration for a given set of rate constants involves just solving the ODEs with intial conditions. This is illustrated for the parameter values (which is the actual parameters based on which data was generated).

This is wrapped into a function whose input is the parameters to be estimated and the output is the residuals. Here there are three concentrations that are fitted. In general, we may want to put different weights to these but in this example, they have the same weight. Also, the package FME has several nice features for fitting and uncertainty estimation. Here I have not used FME but directly tried to do the fit more to learn how this works.

Another way of estimating the uncertainty is using bootstrapping procedure. Here several simulated datasets are generated using the current estimated model and adding random normal noise to each data point (with mean=0 and variance = mean square error from model). Then the parameter is estimated for each simulated dataset. The set of parameters thus generated indicate the uncertainty in parameter estimates.

Next, the percentage of parameter estimates that fall within the ellipse is computed and found to be 93% (expected is 95%). In this case, the ellipsoidal approximation of parameter uncertainty seems adequate. But this might not be the case in general.

As I mentioned previously, the package FME has several functions for parameter fitting, MCMC simulation, sensitivity analysis but was not used in this post. Also, when running multiple simulations using ODE models, it might be better to use the C compiled version of reaction rate function. deSolve package provides that ability and describes it in this vignette. The data and R code is available here.

I don’t have other examples that I have done in R. But one of the reference I find useful is:Chemical Reactor Analysis and Design
It uses Octave as language but hopefully the method in this post can be used to apply to other problems.

I don’t know why this might be happening. One thought as I was looking through the equations was following. It looks like the set of reactions you are modeling is:
Inf -> A (forward rate kinf)
A = B (forward rate k12, reverse rate k21)
A -> C (forward rate ke)
So I am wondering if the expression for r(2) is the following (c[“A”] in the last term instead of c[“C”]) so that mole balance will be satisfied:
r[2]=+kinf*c[“inf”]-k12*c[“A”]+k21*c[“B”]-ke*c[“A”]
Maybe that change might fix things. If you find a solution, please post it in the comment here.

I’ve tried to change rdabbler code and run it. Unfortunately example runs with no problems, but after adjusting equations to my problem, I get a warning messages:
50: In r[2] = (1/k1) * 2.42 * ((1 + 0.26 * cP)/(1 + 0.004 * … :
number of items to replace is not a multiple of replacement length (cP is the variable from the data.frame) and an error message:
Error in chol.default(object$hessian) :
the leading minor of order 1 is not positive definite and further:
Error in solve.default(object$hessian) :
Lapack routine dgesv: system is exactly singular: U[1,1] = 0
My equations are of the form:
r=rep(0,length(c))
r[1]=(500+(4500*cE^8)/(200^8+cE^8))/(1+cP/12.2)-2.42*((1+0.26*cP)/(1+0.004*cE))*c[“pLH”]; #dRP_LH/dt
r[2]=(1/k1)*2.42*((1+0.26*cP)/(1+0.004*cE))*c[“pLH”]-k2*c[“L”] #dL/dt
return(list(r))

I’ve embedded values of some parameters, as the code didn’t like bigger number of parameters to be found that the number of equations.

Reblogged this on the math researcher and commented:
A good tutorial on how to do parameter fitting for models using differential equations. It has helped me with my research, so I am here to pass along the information.

Thanks for this great example and I have been using it quite a lot with my work. I have run into a unique problem in that your example works great for long time series but one of the new data sets I have uses a lot short time series data (n=6) that all needs to be fitted with the same function.

My low tech solution to this has been to create columns for the subsequent time series and add them to the code:

The problem with this is that I quickly run out of letters (e.g. c[“A”]) in addition to the fact that it is extremely inefficient to hard code all of these time series (over 100 in total).

My question is because the paired equations are all the same, is there a solution where I can write them once and have the function apply it to all the subsequent paired time series. I have tested this method against other analytical methods for solving differential equations and the outputs are pretty comparable.

Great example. I am also trying to use R in PK/PD modeling and simulations. This is a great starting point for me.

One question:
When dealing with data from multiple subjects there’s two types of uncertainties: the inter-individual uncertainty and residual error. When we solve a differential equation, is there any good way to add uncertainty to a specific parameter, and add residual error also?

For example: we assume drug concentration dy/dt = -k*y + ε. When we simulate for y changes over time for several subjects, how shall we let k change randomly? Is it also possible to generate a random number (ε) for each observation (this would be each time t).

Thanks for your interest. I don’t have much experience with PK/PD modeling. But it seems like some kind of bayesian modeling would help with what you are looking to do. I came across this youtube video (https://www.youtube.com/watch?v=Q_8NSvPtehc) on a package called PMXStan which does bayesian modeling for PK. But I couldn’t find the actual package in CRAN.

Hi Shankar,
A wonderful example for parameter estimation for model using R. This is my first time to use R, but this tutorial is easy to follow and understand.

I have one more question to ask: if the model initial conditions (cinit) are also considered as the unknow parameters (like k1,k2 and k3). How should I change the above code to achieve this purpose? Would I change the code in ssq=function(parms) as well as its corresponding positions, or do some change in rxnrate function (in fact, it is not necessary to change rxnrate function because cinit is not required in this function. Cinit is only required to get a solution for desolve).

I already tried many changes according to my understanding. However, I did not get the correct solutions at the moment. Hoperfully, I can get your response as soon as possible.

Thank you for your code. Should be very useful for me, I’m new to R. I have problems with running the code (the problem is with the version of R). What version did u use? I have 3.3.2.
Also, I have a problem with plotting. I get an error message:
expdf=melt(df,id.var=”time”,variable.name=”species”,value.name=”conc”)
Error in match.names(clabs, names(xi)) :
names do not match previous names
> ssqres=preddf$conc-expdf$conc
Error: object ‘expdf’ not found
and much more (I couldn’t install minpack on my version of R)

Hi There. Thanks for sharing the code. I was wondering if we have only one column of data (let say we have only data for the second variable c_b), and we want to fit the model with data only related to c_b, and we don’t care about the rest of the variables. how can we change the code just for that, although still we have 3 equations of ODE.
Thanks for your time.