(2) Using Model Specification Files

A Mathematical Model for Exponential Growth or Decay

A quantity undergoes exponential growth or decay when the rate of change of the quantity is proportional to the amount currently present. So, if \(A(t)\) is the amount at time \(t\) and \(r\) is the exponential rate of change (with units of one over time), then \[\begin{equation} \frac{\textrm{d}}{\textrm{d}t}A(t) = r \cdot A(t). \end{equation}\] When \(r\) is positive, the process is known as exponential growth, and when \(r\) is negative, the process is known as exponential decay.

In order to solve an initial value problem involving an exponential growth or decay differential equation, one needs to specify the exponential rate \(r\) and an initial amount \(A_0\) such that \(A(0) = A_0\).

MCSim Model Specification

We will use the GNU MCSim model specification language to implement the exponential model. Recall that this model has a single state variable, \(A\), and two parameters, \(r\) and \(A_0\), the latter of which is an initial condition for the state variable. All of these components of the exponential model can be included in a text file that we refer to as an “MCSim model specification file.” The complete MCSim model specification file for the exponential model, exponential.model, can be found in the extdata subdirectory of the MCSimMod package. Also, a full description of rules and syntax for MCSim model specification files can be found in Section 5 of the GNU MCSim User’s Manual.

The first element of the MCSim model specification file is a listing of the state variables in the ordinary differential equation (ODE) model. We use the text symbol A to represent the state variable \(A\) in the exponential model.

# STATE VARIABLES for the model (for which ODEs are provided).
States = {
    A,        # Amount of substance.
};

Note that the # character indicates the start of a comment. That is, any text following the # character on any line will be ignored when translating the model specification file into machine language. Note also that the state variables should be provided in a comma-delimited list that begins with { and ends with }; (with a semicolon).

The model specification file can also include a listing of the “output variables” for the model. These are variables for which values can be calculated as analytic functions of state variables, “input variables” (which will be described next), and/or parameters. For the exponential model, we will not include any output variables, so we will use a blank list. (This element of the model specification file could be omitted.)

# OUTPUT VARIABLES for the model (which can be obtained at any point in time
# as analytic functions of state variables, input variables, and parameters).
Outputs = {};

Another optional element of the model specification file is a listing of the “input variables” for the model. These variables are independent of other variables and can vary in time. For the exponential model, we will not include any input variables. (This element of the model specification file could be omitted.)

# INPUT VARIABLES for the model (which are independent of other variables, and
# which may vary in time).
Inputs = {};

The next element of the model specification file allows one to name the parameters and to provide default values for those parameters. (Note that the parameter values to be used for specific model simulations can be changed without editing the model specification file.) Recall that the parameters for the exponential model are \(A_0\) and \(r\), for which we use the text symbols A0 and r, respectively.

# PARAMETERS for the model (which are independent of time).
A0 = 0;          # Initial amount.
r = 0;           # Exponential rate of change (time^-1).

The “Initialize” section is the next element of the model specification file. In this section, the values of any static parameters that need to be calculated (e.g., based on the values of other parameters) can be determined and initial values of the state variables can be provided. Note that this section begins with Initialize { and ends with } (with no semicolon).

# MODEL INITIALIZATION section.
Initialize {
    # Assign an initial value for each state variable.
    A = A0;
}

Finally, we have the “Dynamics” section. In this section of the model specification file, the state equations (ODEs) for all state variables should be provided. This section begins with Dynamics { and ends with } (with no semicolon). For each state variable in the model, the state equation should be provided using dt followed by the text symbol for the state variable in parentheses, the = symbol, and then a mathematical expression that represents the time rate of change of the state variable.

# DYNAMICS section.
Dynamics {
    # Time rate of change (ODE) for each state variable.
    dt(A) = r * A;
}

Building the Model

Next, we will use the MCSimMod package to define the ODE model and solve an initial value problem. Solving an initial value problem is synonymous with “performing a simulation” with the model.

First, we load the MCSimMod package as follows.

library(MCSimMod)

Using the following commands, we create an exponential model object (i.e., an instance of the Model class) using the model specification file exponential.model that is included in the MCSimMod package.

# Get the full name of the package directory that contains the example MCSim
# model specification file.
mod_path <- file.path(system.file(package = "MCSimMod"), "extdata")

# Create a model object using the example MCSim model specification file
# "exponential.model" included in the MCSimMod package.
exp_mod_name <- file.path(mod_path, "exponential")
exp_mod <- createModel(exp_mod_name)

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. If necessary (e.g., when the model has not previously been “built”), this step will create output files (with names ending in “.c”, “.o”, “_inits.R”, and “.dll” or “.so”) in a temporary directory. Some of these model files (the ones with names ending in “_inits.R” and “.dll” or “.so”) are used to perform simulations with the model.

# Load the model.
exp_mod$loadModel()
#> C compilation complete. Full details are available in the file C:\Users\dkapraun\AppData\Local\Temp\Rtmp4YuNTn\compiler_output.txt.
#> Hash created and saved in the file C:\Users\dkapraun\AppData\Local\Temp\Rtmp4YuNTn\exponential_model.md5.

Solving an Initial Value Problem (i.e., Performing a Model Simulation)

Next, we can change the parameter values from their default values (which are given in the model specification file) to the values we wish to use for our simulation (i.e., \(A_0 = 100\) and \(r = -0.5\)).

# Change the values of the model parameters from their default values: set the
# initial amount (A0) to 100 and the exponential rate (r) to -0.5.
exp_mod$updateParms(c(A0 = 100, r = -0.5))

# Update the initial value(s) of the state variable(s) based on the updated
# parameter value(s).
exp_mod$updateY0()

Note that executing the command exp_mod$updateParms(c(A0=100, r=-0.5)) updated the parameter values (replacing the default values that were provided in the model specification file) and that executing the command exp_mod$updateY0() updated the initial value of the state variable A using the updated value of the parameter A0. (The class methods updateParms() and updateY0() implement any logic provided in the “Initialize” section of the model specification file.)

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

# Define output times for simulation.
times <- seq(from = 0, to = 20, by = 0.1)

# Run simulation.
out <- exp_mod$runModel(times)

Examining the Results

The final command shown above, out <- exp_mod$runModel(times), performs a model simulation and stores the simulation results in a “matrix” data structure called out. There is one row for each output time, and one column for each state variable (and each output variable when such variables are included in the model). The first five rows of this data structure are shown below. Note that the independent variable, which is \(t\) in the case of the exponential model, is always labeled “time” in the output data structure.

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

exp_mod$parms
#>    A0     r 
#> 100.0  -0.5
exp_mod$Y0
#>   A 
#> 100

Finally, we can create a visual representation of the simulation results. For example, we can plot the amount vs. time using the following command.

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