Problems in making a simulation with "deSolve" for a biological system

Question

I have a hypothesis about a biological system I'm working with, and I need to use the program "deSolve" to formulate differential equations and parameters, in order to simulate the system. But I don't know what's wrong in my code..

This is the system I'm working with:

This is what the simulation graph SHOULD look like:

This is what my graph currently looks like:

This is the code I have currently:

#1. Define Initial conditions
#2. Define time-steps
#3. Define differential-equations for all states
#4.Simulate the model

#----------------------------------------------------------------------------
library(deSolve)

#Define initial conditions
states <- c(R=1, 
            Rp=0,
            RS=1,
            RSp=0,
            S=1)

#Define time steps
times <- c(seq(0, 20, 0.1))

#Define ordinary differential equations for the model
model1 <- function(time, states, parameters) {
  with (as.list(c(states, parameters)), {
    dR = -(R*k1*S) + (Rp*k2) + (Rp*kfeed*RSp)      
    dRp = (R*k1*S) - (Rp*k2) - (Rp*kfeed*RSp)
    dRS = -(RS*k5*Rp) + (RSp*k4)
    dRSp = (RS*k5*Rp) - (RSp*k4) - (Rp*kfeed*RSp)
    dS = 0;
    return(list(c(dR, dRp, dRS, dRSp, dS)))
  })
}

#Simulate the model
simModel1 <- function(parameters) {
  return(as.data.frame(ode(y = states, times = times, func = model1, parms = parameters)))
}

pstart1 <- c(k1=1, 
             k2=0.001,
             kfeed=100,
             k4=0.01,
             k5=0.01)

sim1 <- simModel1(pstart1)

plot(sim1$time, sim1$Rp, type = "l", xlab = "time", ylab = "Rp", col = 'violet')

Does anyone know what might be wrong in my code and how do I fix it so my graph looks like it should (the second image)? :/ — Lina2, Feb 12 '23 at 17:11

tpetzoldt · Answer 1 · 2023-02-12T22:37:22.817

The figure is somewhat out of context, so I am not sure if I understand its syntax correctly. After some creative interpretation of the provided code, it may look something like the code below. As you wrote "something like this", the parameters may differ, so I just extended the simulation time to 300.

library(deSolve)

states <- c(R=1, Rp=0, RS=1, RSp=0, S=1)
times <- seq(0, 300, 0.1)

model1 <- function(time, states, parameters) {
  with (as.list(c(states, parameters)), {
    dR <- -(R*k1*S) + (Rp*k2) + (Rp*kfeed*RSp)      
    dRp <- (R*k1*S) - (Rp*k2) - (Rp*kfeed*RSp) - (RS*k5*Rp)
    dRS <- -(RS*k5*Rp) + (RSp*k4)
    dRSp <- (RS*k5*Rp) - (RSp*k4) - (Rp*kfeed*RSp)
    dS <- -(R*S*k1)
    return(list(c(dR, dRp, dRS, dRSp, dS)))
  })
}

parameters <- c(k1=1, k2=0.001, kfeed=100, k4=0.01, k5=0.01)

sim <- ode(y = states, times = times, func = model1, parms = parameters)

plot(sim)

ode simulation

Problems in making a simulation with "deSolve" for a biological system

1 Answers1