r - eval (expr、envir、enclos) のエラー: オブジェクト 'B' が見つかりません

Question

私はこのコードでほぼ 2 か月間立ち往生しており、どんな助けも大歓迎です。

RのdeSolveパッケージで3つの微分方程式を統合したいと思います.これが私のコードです

library(deSolve)
library(ggplot2)


### Parameters 

D = 0.1  
S0= 6 
c = 2.3 * 10 ^-5 
a = c (0.25, 0.225, 0.2, 0.175, 0.15) # algae maximum growth rate 
H = 1 # algae conversion efficiency 
phi = 7.5 * 10^-8  
beta = 100  
epsilon = 10^-3

M_B =  matrix(c(1-epsilon, epsilon/2, 0,0,0,epsilon, (1-epsilon),    (epsilon/2), 0, 0, 0, epsilon/2, (1-epsilon), epsilon/2, 0, 0, 0 , epsilon/2, (1-epsilon), epsilon,0,0,0, epsilon/2, 1-epsilon),
          nrow=5,
          ncol=5,
          byrow=TRUE)

 M_P =  matrix(c(1-epsilon, epsilon/2,0,0,epsilon, (1-epsilon),(epsilon/2), 0, 0, epsilon/2, (1-epsilon), epsilon, 0,0,  epsilon/2, (1-epsilon)),
          nrow=4,
          ncol=4,
          byrow=TRUE)



 A= matrix(c(1,1,1,1,0,1,1,1,0,0,1,1,0,0,0,1,0,0,0,0), 
      nrow=5, 
      ncol=4,
      byrow=TRUE)

  ## time sequence
 time <- seq(0,1000, by = 1)

# parameters: a named vector
parameters <- c(D = 0.1,
            c = 2.3, 
            H = 1,
            a = c (0.25, 0.225, 0.2, 0.175, 0.15),
            S0= 30, 
            c = 2.3 * 10 ^-5,  
            H = 1, 
            phi = 7.5 * 10^-8,  
            beta = 100,
            epsilon = 10^-3,
            M_B =  matrix(c(1-epsilon, epsilon/2, 0,0,0,epsilon, (1-epsilon), (epsilon/2), 0, 0, 0, epsilon/2, (1-epsilon), epsilon/2, 0, 0, 0 , epsilon/2, (1-epsilon), epsilon,0,0,0, epsilon/2, 1-epsilon),
                          nrow=5,
                          ncol=5,
                          byrow=TRUE),
            M_P =  matrix(c(1-epsilon, epsilon/2,0,0,epsilon, (1-epsilon),(epsilon/2), 0, 0, epsilon/2, (1-epsilon), epsilon, 0,0,  epsilon/2, (1-epsilon)),
                          nrow=4,
                          ncol=4,
                          byrow=TRUE),
            A= matrix(c(1,1,1,1,0,1,1,1,0,0,1,1,0,0,0,1,0,0,0,0), 
                        nrow=5, 
                        ncol=4,
                        byrow=TRUE)) 


    nutrients <- function(t, state, parameters){
 with(as.list(c(state, parameters)),{
 g= a*S / (H + S)
 dS= D*(S0 - S) - c*sum(g,B)
 dB = M_B %*% (g * B) - (phi * (A %*% P)) * B - D*B
 dP= (M_P * beta) %*% (phi*(t(A)%*%B)*P) - (phi*(t(A)%*%B)*P) - D*P   
 return(list(c(dS,dB,dP)))
 })
 }

  out <- ode(y = c(S=30, B=c(10000,0,0,0,0), P=c(100,0,0,0)), times = time,    func = nutrients, parms = parameters)

ただし、このエラーが発生したため、まだ成功していません。

eval (expr、envir、enclos) のエラー: オブジェクト 'B' が見つかりません

私が間違っていることを知っていますか？

アップデート

---

しばらく試した後、質問に対する答えを見つけました。ソリューションとグラフを含む github リンクを後で投稿します。