MCSim under R (Unix)

0.1 Prerequisites

  • R
  • R Studio (optional)
  • MCSim

This exercise was sourced from the file MCSim under R in GNU MCSim. The testing operation system is Lubuntu v17.10 with MCSim v5.6.6 and R v3.4.0 installed. \(\mu\)

Be sure you’re in the right place. Check .Platform$OS.type == "unix" first! Windows user can read this instruction.

Also, we can use this method,

cat(system("lsb_release -a", intern = T), sep = '\n')

1 Solve ODE in pure R

library(deSolve)

model <- function(t, Y, parameters) {
  with (as.list(parameters), {
    dy1 = -k1*Y[1] + k2*Y[2]*Y[3];
    dy3 = k3*Y[2]*Y[2];
    dy2 = -dy1 - dy3;
    list(c(dy1, dy2, dy3));
    })
}

jac <- function (t, Y, parameters) {
  with (as.list(parameters), {
       PD[1,1] <- -k1;
       PD[1,2] <- k2*Y[3];
       PD[1,3] <- k2*Y[2];
       PD[2,1] <- k1;
       PD[2,3] <- -PD[1,3];
       PD[3,2] <- k3*Y[2];
       PD[2,2] <- -PD[1,2] - PD[3,2];

       return(PD) 
    }) 
}


parms <- c(k1 = 0.04, k2 = 1e4, k3=3e7) 

Y <- c(y1 = 1.0, y2 = 0.0, y3 = 0.0)

times  <- c(0, 0.4*10^(0:11))

PD <- matrix(nrow = 3, ncol = 3, data = 0)

out.1 <- ode(Y, times, model, parms = parms, jacfunc = jac)

plot(out.1,type="b",log="x")

2 Compile the MCSim model file and solve by deSolve

2.1 Prepare model file

The testing model simple.model file looks like this:

#------------------------------------------------------------------------------
# the deSolve example model (simple.model)
#------------------------------------------------------------------------------

States = {y0, y1, y2};

Outputs = {yout};

# Parameters
k1 = 1;
k2 = 1;
k3 = 1;

Dynamics {

  dt(y0) = -k1 * y0 + k2 * y1 * y2; 
  dt(y2) =  k3 * y1 * y1; 
  dt(y1) = -dt(y0) - dt(y2);

  yout = y0 + y1 + y2;

} # End of Dynamics

Events {
  y0 = 1;
}

Roots { # here we need inlining, otherwise 'gout' will not be understood 
  Inline ( gout[0] = y[0] - 0.5; );
}

End.

The following step can compile the simple.model file and create the executable program

Make sure that the simple.model is put in the working directory!

mName <- "simple.model"

# Create .c file and use -R to generate "simple.model_inits.R" initialization file
system (paste("mod -R ", mName, " ", mName, ".c", sep = "")) 
system (paste("R CMD SHLIB ", mName, ".c", sep = "")) # create .o and .so files
dyn.load(paste(mName, .Platform$dynlib.ext, sep=""))
source ("simple.model_inits.R")

parms_init <- initParms () # define parameter values
Y_init <- initStates (parms) # define initial state values

# Check the default parameters, initial values
parms; Y
##    k1    k2    k3 
## 4e-02 1e+04 3e+07
## y1 y2 y3 
##  1  0  0
# list output variable, not a nice name will will call it "Sum" below
Outputs 
## [1] "yout"
# set output times, time 0 must be given
times <- c(0, 0.4*10^(0:11)) # This setting is same as section 1

out.2.1 <- ode(Y_init, times, func = "derivs", parms = parms_init, 
               jacfunc = "jac", dllname = "simple.model", 
               initfunc = "initmod", nout = 1, outnames = "Sum")

plot(out.2.1, log="x",type="b")

We can’t find any change in the dynamic process. It is due to the all initial value equal to 0.

2.2 Customize the parameter and initial values

# Use the samee parameters and initial values in section 1
out.2.2 <- ode(Y, times, func = "derivs", parms = parms,
               dllname = mName, initfunc = "initmod",
               nout = 1, outnames = "Sum")

# Here is the output data frame
head(out.2.2)
##       time        y1           y2         y3 Sum
## [1,] 0e+00 1.0000000 0.000000e+00 0.00000000   1
## [2,] 4e-01 0.9851723 3.386399e-05 0.01479384   1
## [3,] 4e+00 0.9055179 2.240437e-05 0.09445974   1
## [4,] 4e+01 0.7158298 9.185501e-06 0.28416101   1
## [5,] 4e+02 0.4505174 3.222887e-06 0.54947938   1
## [6,] 4e+03 0.1831994 8.942298e-07 0.81679972   1
# Plotting
plot(out.2.2, log="x",type="b")

# Use "which" to select specific variables
plot(out.2.2, log="x",type="b", which = c("y1", "y2", "y3"))

3 Compile the MCSim model file and solve by MCSim

3.1 Prepare in file

In this part, we need to put simple.in file in the working directory. The in file define the initial and parameter values, respectively. Also,we need to set the output time steps and intervals. The source code looks like this:

#------------------------------------------------------------------------------
# the deSolve example model run with GNU MCSim (simple.in)
#------------------------------------------------------------------------------

Integrate (Lsodes, 1E-7, 1E-7, 1);

Simulation {

  y0 = 1;

  # Parameters
  k1 = 0.04;
  k2 = 1e4;
  k3 = 3e7;

  Print (y0, 0, 4e-1, 4e+0, 4e+1, 4e+2, 4e+3, 4e+4, 4e+5, 4e+6, 4e+7, 4e+8, 
         4e+9, 4e+10);
  Print (y1, 0, 4e-1, 4e+0, 4e+1, 4e+2, 4e+3, 4e+4, 4e+5, 4e+6, 4e+7, 4e+8, 
         4e+9, 4e+10);
  Print (y2, 0, 4e-1, 4e+0, 4e+1, 4e+2, 4e+3, 4e+4, 4e+5, 4e+6, 4e+7, 4e+8, 
         4e+9, 4e+10);

  Print (yout, 0, 4e-1, 4e+0, 4e+1, 4e+2, 4e+3, 4e+4, 4e+5, 4e+6, 4e+7, 4e+8, 
         4e+9, 4e+10);

}

End.

The following step compile the simple.model file and create the executable program named mcsim.simple

Name <- "simple"
mName <- "simple.model"

# compile the "simple.model" to "mcsim.simple"
system(paste("makemcsim ", mName, sep = ""))

# run the simulation input file
# put its name here
inName <- "simple.in"

# run!
system(paste("./mcsim.", Name, " ", inName, sep = ""))

# the output file is called out by default, load it
out.3 <- read.delim("sim.out", skip = 2)

head (out.3)
##    Time       y0          y1        y2 yout
## 1 0e+00 1.000000 0.00000e+00 0.0000000    1
## 2 4e-01 0.985172 3.38640e-05 0.0147937    1
## 3 4e+00 0.905519 2.24048e-05 0.0944587    1
## 4 4e+01 0.715826 9.18549e-06 0.2841650    1
## 5 4e+02 0.450518 3.22290e-06 0.5494780    1
## 6 4e+03 0.183202 8.94234e-07 0.8167970    1
# look

lo = layout(t(matrix(1:4, ncol = 2)), respect = TRUE)
plot(out.3[,1], out.3[,2], log="x",type="b", xlab = "time", ylab = "", main = "y1")
plot(out.3[,1], out.3[,3], log="x",type="b", xlab = "time", ylab = "", main = "y2")
plot(out.3[,1], out.3[,4], log="x",type="b", xlab = "time", ylab = "", main = "y3")

Related

comments powered by Disqus