rodedifferential-equationsrstandesolve

Error in eval(expr, envir, enclos) : object 'B' not found

I have been stuck with this code for almost two months, and any help would be highly appreciated.

I would like to integrate three differential equations with the deSolve package in R. Here is my code.

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)

And these are the mathematical equations that i would like to explore

However I have not succeeded yet since I take this error:

Error in eval(expr, envir, enclos) : object 'B' not found

Do you have any idea what I am doing wrong?

Update

After trying some time I have found the answer to the question. I ll post a github link in a bit with the solution and the graphs

Solution

  • the code runs perfectly like this

    library(deSolve);library(ggplot2);library(reshape2) # load packages deSolve is for solving ODE

##############################################################################################################
#
# model parameters
#
##############################################################################################################

a = c(0.25, 0.225, 0.2, 0.175, 0.15) #alga growth rate for each one of the five types
epsilon = 10^-3                      #mutation rate

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)           #mutation-transition matrix for the host the sum of each column should be one


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)            #mutation-transition matrix for the virus the sum of each column should be one


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,                   #infection matrix where 1 means infection and 0 resistance
          byrow=TRUE)


# here I put all the parameters into a vector 
parameters <- c(D = 0.1,
                epsilon = 10^-3,
                a = c(0.25, 0.225, 0.2, 0.175, 0.15),
                S0 = 10, # inflow rate 
                c = 2.3 * 10 ^-5, # algae conversion coefficiency 
                H = 1, # algae conversion efficiency 
                phi = 7.5 * 10^-8,  #  virus adsorption rate 
                beta = 200, # virus burst size 
                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
#
##############################################################################################################

time <- seq(0,300, by = 0.01)

##############################################################################################################
#
# initial values
#
##############################################################################################################

yini= c(S=30, B=c(100,0,0,0,0), P=c(1,0,0,0))

##############################################################################################################
#
# function
#
##############################################################################################################


nutrients <- function(t, yini, parameters) {
  with(as.list(c(yini, parameters)), {
    g = ((a*S) / (H + S))
    B = yini[paste("B",1:5, sep = "")]   #### like this you convert B from a vector to an array and this is what the program needs 
    P = yini[paste("P",1:4, sep = "")]
    
    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)))  
  })
}

##############################################################################################################
#
# solve the model
#
##############################################################################################################

out <- ode(y = yini, times = time, func = nutrients, parms = parameters, method="lsoda",atol=10^-15,rtol=10^-15)