Solving the system of differential equations in R

Question

I have a simple flow model in R. It comes down to two differential equations that model two state variables in the model, we will call them A and B They are calculated as simple difference equations of the four-component flux1-flux4 , 5 parameters p1-p5 and 6th parameter of_interest , which can take values ​​between 0-1.

 parameters<- c(p1=0.028, p2=0.3, p3=0.5, p4=0.0002, p5=0.001, of_interest=0.1) state <- c(A=28, B=1.4) model<-function(t,state,parameters){ with(as.list(c(state,parameters)),{ #fluxes flux1 = (1-of_interest) * p1*(B / (p2 + B))*p3 flux2 = p4* A #microbial death flux3 = of_interest * p1*(B / (p2 + B))*p3 flux4 = p5* B #differential equations of component fluxes dAdt<- flux1 - flux2 dBdt<- flux3 - flux4 list(c(dAdt,dBdt)) }) 

I would like to write a function to take the derivative of dAdt with respect to of_interest , set the derivative equation to 0, and then rearrange and solve for the value of_interest . This will be the value of_interest , which maximizes the dAdt function.

So far, I have managed to solve the model in a stable state according to the possible values of_interest , to demonstrate that there should be a maximum.

 require(rootSolve) range<- seq(0,1,by=0.01) for(i in range){ of_interest=i parameters<- c(p1=0.028, p2=0.3, p3=0.5, p4=0.0002, p5=0.001, of_interest=of_interest) state <- c(A=28, B=1.4) ST<- stode(y=y,func=model,parms=parameters,pos=T) out<- c(out,ST$y[1]) 

Then build:

 plot(out~range, pch=16,col='purple') lines(smooth.spline(out~range,spar=0.35), lwd=3,lty=1) 

How can I analytically solve the value of_interest that maximizes dAdt in R? If an analytical solution is impossible, how can I know, and how can I solve it numerically?

Update: I think this problem can be solved with the deSolve package in R linked here , however I am having problems implementing it using my example.

Answer 1

Your equation in B(t) simply split, because you can separate B(t) from which you can get

 B(t) = C * exp{-p5 * t} * (p2 + B(t)) ^ {of_interest * p1 * p3} 

This is an implicit solution for B(t) , which we will solve in different ways.

You can solve for C given initial value of B I suppose t = 0 initially? In this case

 C = B_0 / (p2 + B_0) ^ {of_interest * p1 * p3} 

It also gives a slightly nicer expression for A(t) :

 dA(t) / dt = B_0 / (p2 + B_0) * p1 * p3 * (1 - of_interest) * exp{-p5 * t} * ((p2 + B(t) / (p2 + B_0)) ^ {of_interest * p1 * p3 - 1} - p4 * A(t) 

This can be solved by integrating the factor (= exp{p4 * t} ), by numerically integrating the term with B(t) . We define the lower limit of the integral as 0 so that we never need to evaluate B outside the range [0, t] , which means that the integrating constant is just A_0 and therefore:

 A(t) = (A_0 + integral_0^t { f(tau; parameters) d tau}) * exp{-p4 * t} 

The main point of B(t) governs everything in this system - the approach will be: decide for the behavior of B(t) , then use this to find out what happens to A(t) , and then maximize it.

First, the "external" parameters; we also need nleqslv to get B :

 library(nleqslv) t_min <- 0 t_max <- 10000 t_N <- 10 #we'll only solve the behavior of A & B over t_rng t_rng <- seq(t_min, t_max, length.out = t_N) #I'm calling of_interest ttheta ttheta_min <- 0 ttheta_max <- 1 ttheta_N <- 5 tthetas <- seq(ttheta_min, ttheta_max, length.out = ttheta_N) B_0 <- 1.4 A_0 <- 28 #No sense storing this as a vector when we'll only ever use it as a list parameters <- list(p1 = 0.028, p2 = 0.3, p3 = 0.5, p4 = 0.0002, p5 = 0.001) 

Hence the main plan:

• Given the values ​​of the parameters (in particular, ttheta ), solve for BB over t_rng by solving nonlinear equations
• Given BB and parameter values, solve for AA over t_rng by numerical integration
• Provided by AA and your expression for dAdt, plug and maximize.

derivatives <sapply (tthetas, function (th) {#append current ttheta params <- c (parameters, ttheta = th)

 #declare a function we'll use to solve for B (see above) b_slv <- function(b, t) with(params, b - B_0 * ((p2 + b)/(p2 + B_0)) ^ (ttheta * p1 * p3) * exp(-p5 * t)) #solving point-wise (this is pretty fast) # **See below for a note** BB <- sapply(t_rng, function(t) nleqslv(B_0, function(b) b_slv(b, t))$x) #this is f(tau; params) that I mentioned above; # we have to do linear interpolation since the # numerical integrator isn't constrained to the grid. # **See below for note** a_int <- function(t){ #approximate t to the grid (t_rng) # (assumes B is monotonic, which seems to be true) # (also, if t ends up negative, just assign t_rng[1]) t_n <- max(1L, which.max(t_rng - t >= 0) - 1L) idx <- t_n:(t_n+1) ts <- t_rng[idx] #distance-weighted average of the local B values B_app <- sum((-1) ^ (0:1) * (t - ts) / diff(ts) * BB[idx]) #finally, f(tau; params) with(params, (1 - ttheta) * p1 * p3 * B_0 / (p2 + B_0) * ((p2 + B_app)/(p2 + B_0)) ^ (ttheta * p1 * p3 - 1) * exp((p4 - p5) * t)) } #a_int only works on scalars; the numeric integrator # requires a version that works on vectors a_int_v <- function(t) sapply(t, a_int) AA <- exp(-params$p4 * t_rng) * sapply(t_rng, function(tt) #I found the subdivisions constraint binding in some cases # at the default value; no trouble at 1000. A_0 + integrate(a_int_v, 0, tt, subdivisions = 1000L)$value) #using the explicit version of dAdt given as flux1 - flux2 max(with(params, (1 - ttheta) * p1 * p3 * BB / (p2 + BB) - p4 * AA))}) Finally, simply run `tthetas[which.max(derivs)]` to get the maximizer. 

Note:

This code is not optimized for efficiency. There are several places where there are some potential accelerations:

• it might be faster to run a recursive solution algorithm, because it will converge faster with better initial guesses - using the previous value instead of the initial value is definitely better
• It would be faster to simply use Riemann sums for integration; the compromise is accurate, but there should be a penalty if you have a fairly dense grid. One beauty of Riemann is that you don’t have to interpolate at all, but numerically it is a simple linear algebra. I ran this with t_N == ttheta_N == 1000L , and it started in a few minutes.
• Perhaps you can simply vectorize a_int , and not just sapply on it, which accompanies acceleration, a more direct appeal to BLAS.
• Download other little things. Pre-compute ttheta * p1 * p3 , since it is reused so much, etc.

I didn’t worry about including any of these things because, frankly, you better translate this into a faster language - Julia is my favorite favorite favorite, but of course, R speaks well with C ++, C, Fortran etc.