求解R中的微分方程组

我在R中有一个简单的通量模型.它归结为两个微分方程,模拟模型中的两个状态变量,我们称之为A和B.它们被计算为四分量通量flux1-flux4,5个参数p1-p5和第六个参数的简单差分方程of_interest,其可以取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))
  })

我想编写一个函数来取的导数dAdt相对于of_interest,派生方程设置为0,然后重新排列并求解的值of_interest.这将是of_interest最大化函数的参数的值dAdt.

到目前为止,我已经能够在稳定状态下解决模型,跨越可能的值of_interest来证明应该有一个最大值.

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])

然后绘图:

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

我怎样才能在分析上解决R of_interest中最大化的价值dAdt？如果无法获得分析解决方案,我怎么知道,如何以数字方式解决这个问题呢？

更新:我认为这个问题可以通过R中的deSolve包解决,链接在这里,但是我在使用我的特定示例实现它时遇到了麻烦.