MCSim under R

The post was updated on 2019-12-31 under Windows 10 x64.

Introduction

Using GNU MCSim under Windows is more complicated than the Unix-like System, such as macOS and Linux. Also, it does not provide the function to analyze and visualize the output result. Therefore, I write this post to demo how to run GNU MCSim under Windows OS and conduct the follow-up analysis with R. This exercise was sourced from the file MCSim under R in GNU MCSim.

In this post, I used Rtools (v3.5) which has some essential set of Unix/Linux tools to conduct the model compilation. Also, the GNU MCSim (v6.1.0) is used to do the following modeling and simulation works. The portable version are put in my GitHub repo (https://github.com/nanhung/MinGM)

wd <- getwd()
url <- "https://github.com/nanhung/MinGM/archive/v6.1.0.zip"
tf <- tempfile()
download.file(url, tf, mode = "wb")
unzip(zipfile = tf, exdir = wd)

If you have installed devtools you can use devtools::find_rtools() to check the avaliable of Rtools.

devtools::find_rtools()
## [1] TRUE

Also, we need to chek the Rtools is in the PATH by using

Sys.getenv("PATH")
## [1] "C:\\Program Files\\R\\R-3.6.2\\bin\\x64;C:\\Windows\\system32;C:\\Windows;C:\\Windows\\System32\\Wbem;C:\\Windows\\System32\\WindowsPowerShell\\v1.0\\;C:\\Windows\\System32\\OpenSSH\\;C:\\Program Files\\NVIDIA Corporation\\NVIDIA NvDLISR;C:\\Program Files\\Git\\cmd;C:\\Users\\nan_1\\AppData\\Local\\Microsoft\\WindowsApps;;C:\\Users\\nan_1\\AppData\\Local\\Programs\\Microsoft VS Code\\bin"

If you can not find the Rtools in the PATH. such as, c:\\Rtools\\bin and c:\\Rtools\\mingw_64\\bin. You will need to manually add it as follow,

Sys.setenv(PATH = paste("c:\\Rtools\\bin", Sys.getenv("PATH"), sep=";"))
Sys.setenv(PATH = paste("c:\\Rtools\\mingw_64\\bin", Sys.getenv("PATH"), sep=";"))

Then, use following command to check the function.

Sys.which("gcc")
##                                  gcc 
## "c:\\Rtools\\mingw_64\\bin\\gcc.exe"
Sys.which("make")
##                        make 
## "c:\\Rtools\\bin\\make.exe"

Create GNU MCSim program

The following steps is used to generate the executable program (mod.exe) that can traslate the format of model code to C.

libSBML = "" # libSBML = "-lsbml" 
# uncomment and execute if you have and want to use libSBML.

# create mod.exe in "mod" folder
system(paste("gcc -o MinGM-6.1.0/mod.exe MinGM-6.1.0/mod/*.c ", libSBML, sep = "")) 

NOTE: In this case, we don’t need to load libSBML. Therefore, check the config.h file in mod folder. We can see the line 26 in the config.h file like THIS. The #define HAVE_LIBSBML 1 had been disabled. It can prevent the problem when compile to the mod.exe file. (Thanks Yu-Sheng Lin for provding this issue; Mar. 23, 2018).

Finally, check if we successfully generate the program mod.exe.

file.exists("MinGM-6.1.0/mod.exe")
## [1] TRUE

If everything is going well, we can move on to next section.

R package deSolve solution

Solve ODE in pure R

This section is a demo of solving the ODE in R environment with deSolve package.

# Load deSolve package
library(deSolve)
# Create the ODE function based on the **deSolve** format
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) 
    }) 
}
# Define the parameter values and some conditions
parms <- c(k1 = 0.04, k2 = 1e4, k3=3e7) # parameter values
Y <- c(y1 = 1.0, y2 = 0.0, y3 = 0.0) # Initial conditions
times  <- c(0, 0.4*10^(0:11)) # Time steps
PD <- matrix(nrow = 3, ncol = 3, data = 0)
# Solve ODE 
out.1 <- ode(Y, times, model, parms = parms, jacfunc = jac)
head(out.1)
##       time        y1           y2         y3
## [1,] 0e+00 1.0000000 0.000000e+00 0.00000000
## [2,] 4e-01 0.9851723 3.386399e-05 0.01479384
## [3,] 4e+00 0.9055179 2.240437e-05 0.09445974
## [4,] 4e+01 0.7158298 9.185501e-06 0.28416101
## [5,] 4e+02 0.4505174 3.222887e-06 0.54947938
## [6,] 4e+03 0.1831994 8.942298e-07 0.81679972
# Visualization
plot(out.1,type="b",log="x")
## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from
## logarithmic plot

## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from
## logarithmic plot

## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from
## logarithmic plot

Use deSolve with GNU MCSim model file

Here, we want to use model code that was based on the GNU MCSim format. The model code look like following format,

##################
## 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;

} 

Events {
  y0 = 1;
}

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

End.

First, we need to create model script that is based on C (simple.model.c) and generate initialization file (simple.model_inits.R).

system("MinGM-6.1.0/mod.exe -R simple.model simple.model.c") 
# create simple.model.o and simple.model.dll files
system ("R CMD SHLIB simple.model.c") 

# load simple.model.dll
dyn.load(paste("simple.model", .Platform$dynlib.ext, sep="")) 

Note, the compile the GNU MCSim model, needs to be done each time you modify the model.

source ("simple.model_inits.R")

# define parameter values
parms_init <- initParms () 

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

# set output times, time 0 must be given
times <- c(0, 4e-1, 4e+0, 4e+1, 4e+2, 4e+3, 4e+4, 4e+5, 4e+6, 4e+7, 4e+8, 4e+9)
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")
## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from
## logarithmic plot

## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from
## logarithmic plot

## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from
## logarithmic plot

## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from
## logarithmic plot

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

parms_init
## k1 k2 k3 
##  1  1  1
Y_init
## y0 y1 y2 
##  0  0  0

Parameters and initial values

Since the parameters and initial state values should define before simulation, we need to use the following code to assign them.

# Define the parameters and initial values
parms <- c(k1 = 0.04, k2 = 1e4, k3=3e7) 
Y <- c(y0 = 1.0, y1 = 0.0, y2 = 0.0)

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

# Here is the output data frame
head(out.2.2)
##       time        y0           y1         y2 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")
## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from
## logarithmic plot

## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from
## logarithmic plot

## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from
## logarithmic plot

## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from
## logarithmic plot

GNU MCSim solution

Here we want to use GNU MCSim to solve ODE. In this part, we need to have the input file (simple.in) in the working directory. The input file can define the initial and parameter values, respectively. Also,we need to set the output time steps and intervals. The source code looks like this:

###############
## simple.in ##
###############

Simulation {

  # Initial value
  y0 = 1;

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

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

End.

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

mcsim <- "MinGM-6.1.0"
model <- "simple.model"
input <- "simple.in"
system(paste0(mcsim, "/mod.exe ", model, " ", model, ".c")) 

Unlike above exercice, we need to remove the -R flag that can only create c file. Then, compile the simple.model.c to mcsim.simple.model.exe.

system(paste0("gcc -O3 -I.. -I ", mcsim,"/sim -o mcsim.", model, ".exe ", 
       model, ".c ", mcsim, "/sim/*.c -lm "))
# run!
system(paste0("./mcsim.", model, ".exe ", input))
# 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 4e-01 0.9851720 3.38640e-05 0.0147937    1
## 2 4e+00 0.9055160 2.24043e-05 0.0944616    1
## 3 4e+01 0.7158290 9.18560e-06 0.2841620    1
## 4 4e+02 0.4505040 3.22272e-06 0.5494920    1
## 5 4e+03 0.1831970 8.94202e-07 0.8168020    1
## 6 4e+04 0.0389846 1.62182e-07 0.9610150    1
# look
par(mfrow=c(2,2))
plot(out.3[,1], out.3[,2], log="x",type="b", xlab = "time", ylab = "", main = "y0")
plot(out.3[,1], out.3[,3], log="x",type="b", xlab = "time", ylab = "", main = "y1")
plot(out.3[,1], out.3[,4], log="x",type="b", xlab = "time", ylab = "", main = "y2")
plot(out.3[,1], out.3[,5], log="x",type="b", xlab = "time", ylab = "", main = "yout")

In addition, creating the R function can help us conduct the simulation more easily,

# Create function
mingm <- function(mcsim, model, input){
  system(paste0(mcsim, "/mod.exe ", model, " ", model, ".c")) 
  paste0("gcc -O3 -I.. -I ", mcsim,"/sim -o mcsim.", model, ".exe ",
         model, ".c ", mcsim, "/sim/*.c -lm ")
  system(paste0("./mcsim.", model, ".exe ", input))
}
mingm(mcsim = "MinGM-6.1.0", model = "simple.model", input = "simple.in")
out.4 <- read.delim("sim.out", skip = 2)
head (out.4)
##    Time        y0          y1        y2 yout
## 1 4e-01 0.9851720 3.38640e-05 0.0147937    1
## 2 4e+00 0.9055160 2.24043e-05 0.0944616    1
## 3 4e+01 0.7158290 9.18560e-06 0.2841620    1
## 4 4e+02 0.4505040 3.22272e-06 0.5494920    1
## 5 4e+03 0.1831970 8.94202e-07 0.8168020    1
## 6 4e+04 0.0389846 1.62182e-07 0.9610150    1

Session info

sessionInfo()
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] deSolve_1.25
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.3        compiler_3.6.2    prettyunits_1.0.2 remotes_2.1.0    
##  [5] tools_3.6.2       testthat_2.3.1    digest_0.6.23     pkgbuild_1.0.6   
##  [9] pkgload_1.0.2     evaluate_0.14     memoise_1.1.0     rlang_0.4.2      
## [13] cli_2.0.0         rstudioapi_0.10   yaml_2.2.0        blogdown_0.17    
## [17] xfun_0.11         withr_2.1.2       stringr_1.4.0     knitr_1.26       
## [21] fs_1.3.1          desc_1.2.0        devtools_2.2.1    rprojroot_1.3-2  
## [25] glue_1.3.1        R6_2.4.1          processx_3.4.1    fansi_0.4.0      
## [29] rmarkdown_2.0     bookdown_0.16     sessioninfo_1.1.1 callr_3.4.0      
## [33] magrittr_1.5      backports_1.1.5   ps_1.3.0          ellipsis_0.3.0   
## [37] htmltools_0.4.0   usethis_1.5.1     assertthat_0.2.1  stringi_1.4.3    
## [41] crayon_1.3.4

Related