(1) Introduction to MCSimMod

Dustin Kapraun

2025-08-19

The MCSimMod package can be used to define ordinary differential equation (ODE) models and solve initial value problems (IVPs) for such models. To define an ODE model, one first creates a “model specification” as a text string or a text file. (More details about the structure of the model specification text are provided in a separate vignette: “Model Specification”.) The model specification text can then be used to create a model object (i.e., an instance of the MCSimMod Model class) via the MCSimMod function createModel().

Consider, for example, a simple ODE model given by \[\begin{align} \frac{\textrm{d}}{\textrm{d}t}y(t) &= m, \\ y(0) &= y_0. \end{align}\] We can create a model object corresponding to this ODE model using a text string that both provides the ODE for the state variable, \(y\), and sets the (default) values of the model parameters to \(y_0 = 2\) and \(m = 0.5\).

First, we load the MCSimMod package as follows.

library(MCSimMod)

Then, we define the model string as follows.

mod_string <- "
States = {y};
y0 = 2;
m = 0.5;
Initialize {
    y = y0;
}
Dynamics {
    dt(y) = m;
}
End.
"

Next, we create a Model object called model using the following command.

model <- createModel(mString = mod_string)

Once the model object is created, we can “load” the model (so that it’s ready for use in a given R session) as follows.

model$loadModel()

Then, we can perform a simulation that provides results for the desired output times (\(t = 0, 0.1, 0.2, \ldots, 20.0\)) using the following commands.

times <- seq(from = 0, to = 20, by = 0.1)
out <- model$runModel(times)

The final command shown above performs the model simulation and stores the simulation results in a “matrix” data structure called out. The first five rows of this data structure are shown below. Note that the independent variable, which is \(t\) in the case of the linear model we’ve created here, is always labeled “time” in the output data structure.

time y
0.0 2.00
0.1 2.05
0.2 2.10
0.3 2.15
0.4 2.20

We can examine the parameter values and initial conditions that were used for this simulation with the following commands.

model$parms
#>  y0   m 
#> 2.0 0.5
model$Y0
#> y 
#> 2

Finally, we can create a visual representation of the simulation results. For example, we can plot the value of \(y\) vs. time (\(t\)) using the following command.

# Plot simulation results.
plot(out[, "time"], out[, "y"],
  type = "l", lty = 1, lwd = 2, xlab = "Time",
  ylab = "y"
)