MCSim under R (Delay differentials example)

0.1 Prerequisites

The preliminary setting was similiar with the previous post - MCSim under R (Unix). The exercise was sourced from the example files in GNU MCSim 6.0.0. All testing was performed in Lubuntu v17.10 with R v3.4.0. The model file perc_delay.model can be also found Here. In the model file, the CalcDelay(Q_liv, tau) was used to defind the delay variable (liver) after the time (tau).

0.1.1 load desolve package

library(deSolve)

1 Compiling Model

1.1 Generating code for linking with R desolve package.

  1. Create C model file perc_delay.c and R parameter initialization file perc_delay_inits.R
system("mod -R perc_delay.model perc_delay.c") 
  1. Compile the perc_delay.c to create perc_delay.o and perc_delay_inits.so (shared object) for dynamic loading
system("R CMD SHLIB perc_delay.c") 
  1. Load perc_delay_inits.so
dyn.load(paste0("perc_delay",.Platform$dynlib.ext)) 

2 Define Parameters and Variables

Due to the model parameters and variables can be generated in perc_delay_inits.R, we can load initParms and initStats functions in perc_delay_inits.R file.

source("perc_delay_inits.R") 

Here, we want to redefine some input parameters. We can use the following method to replace the default value. In this case, we use the delayed time tau = 100 min, inhaled concentration = 1 ppma and exposure duration = 240 min (4 hr) as our setting.

newParms <- c(InhMag = 1, Period=1e10, Exposure=240, tau=100)
parms <- initParms(newParms=newParms)

For the output variables, we select all variables as default.

# Define the target state variables
Y <- initStates(parms=parms)

3 Create the Exposure Matrix

The exposue matrix include two part: Inhalation and Ingestion. We need to generate the list of exposure by combine these two exposure matricies.

3.1 Build inhalation exposure

We firstly set the periodic inhalation exposure scenario by using following step,

mag <- parms["InhMag"] # Set input
period <- 1e10
inittime <- 0 # Exposure start from 0
exposuretime <- 120 # Exposure end at 2 hr

# The output time points
times <- c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 239.99, 240, 250, 260, 270, 280, 290, 300, 310, 320, 330, 339.9, 340, 350) 

mintime <- min(times)
maxtime <- max(times)

# Define exposure scenarios: exposure and non-exposure
nperiods <- ceiling(maxtime / period) + 1 

col1 <- rep(inittime, 2 * nperiods)
I <- 1:length(col1)
J <- which(I %% 2 == 1)
col1[J]  <- col1[J] + (ceiling((J / 2)) - 1) * period
J <- which(I %% 2 == 0)
col1[J] <- col1[J-1] + exposuretime
PerExpo1 <- cbind(col1, rep(c(mag,0), length(col1)/2))

# The matrix of inhalation exposure
PerExpo1
##           col1  
## InhMag 0.0e+00 1
##        1.2e+02 0
## InhMag 1.0e+10 1
##        1.0e+10 0

Here, we can see the setting expoure is start at 0 and 1e10 then stop at 240 and 1e10, respectively.

3.2 Plot exposure scenario

The exposure scenario can be visualized by following method,

plot(PerExpo1[1:2,], type = "s", xlim=c(0,maxtime), main = "C_inh", xlab = "Time (min)")

4 Build Ingestion dose Matrix

4.1 Setting 1: without constant exposure

Here we generate the matrix for ingestion dose,

mag1 <- 0 # the backgroud exposure dose = 0
IngDose <- rbind(c(min(times), mag1), c(max(times) + 1, 0))
IngDose
##      [,1] [,2]
## [1,]    0    0
## [2,]  351    0

As above, we can find that there is no given exposure dose from start to the end (0 to 351 min). We can generate the list of exposure by combine the PerExpo and IngDose. After that, we can further apply the solver for delay differential equations dede from desolve package.

# Add code to implement in the inputs built into the model.
Forcings1 <- list(PerExpo1, IngDose)
out1 <- dede(Y, times, func="derivs", parms=parms, dllname="perc_delay",
            initfunc="initmod", nout=length(Outputs), outnames=Outputs,
            initforc="initforc",
            fcontrol=list(method="constant", rule=1, f=0),
            forcings=Forcings1)

4.2 Setting 2: with constant exposure

Here, we assume the constant exposure concentration = 0.2 ppm.

PerExpo2 <- PerExpo1
PerExpo2[2,2] <- 0.2 # 0.2 ppm as constant exposure 
Forcings2 <- list(PerExpo2, IngDose)
out2 <- dede(Y, times, func="derivs", parms=parms, dllname="perc_delay",
            initfunc="initmod", nout=length(Outputs), outnames=Outputs,
            initforc="initforc",
            fcontrol=list(method="constant"),
            forcings=Forcings2)

We can do the inspection to check the output result by total quantity as,

stop_point <- which(out1[,"time"] == 240)

Then, check the exposure dose without constant dose. The output should approximate to 0.5 (ppm).

out1[,"Q_tot"][stop_point] / (240 * 5.6 * 0.006778) 
## [1] 0.4999836

Check the exposure dose with constant dose. The output should approximate to 0.6 (ppm).

out2[,"Q_tot"][stop_point] / (240 * 5.6 * 0.006778)
## [1] 0.5999803

5 Visualization

Plot and compare the simulation result

par(mfrow = c(3,5), mar = c(4,2,3,1))
for(i in 2:ncol(out1)){
  plot(out2[,1], out2[,i], type="l", xlab = "Time", ylab="" , 
       col = 2, main = names(out1[3,i])) 
  lines(out1[,1], out1[,i], type = "l", col = 3, lty = 2)
}

The green dashed line is periodic exposure without constant level. The red solid line is periodic exposure with constant level The simulation of delay output can obtain in old_Q_liv, Q_met and Pct_metabolized.

Related

comments powered by Disqus