---

title: "dMod - Dynamic Modeling and Parameter Estimation in R"

output:

github_document:

toc: false

---

```{r setup, include=FALSE}

knitr::opts_chunk$set(

fig.width = 8,

fig.height = 3,

out.width = "100%",

dpi = 150

)

```

<!-- badges: start -->

[](https://github.com/JetiLab/dMod/actions/workflows/R-CMD-check.yaml)

<!-- badges: end -->

The dMod package is a framework that provides functions to generate ODEs of reaction networks, parameter transformations, observation functions, residual functions, etc. The framework follows the paradigm that derivative information should be used for optimization whenever possible. Therefore, all major functions produce and can handle expressions for symbolic derivatives.

## System requirements

dMod uses the package [cOde](https://github.com/dkaschek/cOde) to set up ODE models as compiled C code (deSolve) or [CppODE](https://github.com/simonbeyer1/CppODE) to autogenerate C++ code (Boost.Odeint). This means that **C and C++ compilers** are required on the system. On Linux, the compilers are installed by default. Windows users need to install [RTools](https://cran.r-project.org/bin/windows/Rtools/).

For **parallelization**, dMod uses `mclapply()` on Linux and macOS. On Windows, parallelization is implemented via the `foreach` package using `%dopar%`.

## Installation from GitHub

To **install dMod from the GitHub repository**, it is convenient to use **RStudio**.

1. Create a new project via\

**File → New Project → Version Control → Git**.

2. Use the repository URL\

```

jetilab/dMod

```

and create the project.

3. Install all required package dependencies (including GitHub dependencies specified via `Remotes`) by running:

``` r

remotes::install_deps(dependencies = TRUE)

```

4. Open the **Build** tab in RStudio and click **Build and Reload** (or **Install**).

Once these steps are completed, it should be possible to run the following example.

------------------------------------------------------------------------

If **PEtab support** is required, **libSBML** is needed in addition. Installation and usage instructions can be found in the wiki under [Support for PEtab](https://github.com/dkaschek/dMod/wiki/Support-for-PEtab).

------------------------------------------------------------------------

## Simple example: enzyme kinetics

### Load required packages

```{r libraries, message=FALSE}

library(dMod)

library(ggplot2)

```

### Generate an ODE model of enzyme kinetics with enzyme degradation

```{r model}

# Reactions

f <- NULL

f <- addReaction(f,

from = "Enz + Sub",

to = "Compl",

rate = "k1*Enz*Sub",

description = "production of complex")

f <- addReaction(f,

from = "Compl",

to = "Enz + Sub",

rate = "k2*Compl",

description = "decay of complex")

f <- addReaction(f,

from = "Compl",

to = "Enz + Prod",

rate = "k3*Compl",

description = "production of product")

f <- addReaction(f,

from = "Enz",

to = "" ,

rate = "k4*Enz",

description = "enzyme degradation")

# ODE model

model <- odemodel(f, modelname = "enzymeKinetics")

# Prediction function

x <- Xs(model)

```

### Define observables and generate observation function `g`

```{r observables}

observables <- eqnvec(

product = "Prod",

substrate = "(Sub + Compl)",

enzyme = "(Enz + Compl)"

)

# Generate observation functions

g <- Y(observables, x, compile = TRUE, modelname = "obsfn", attach.input = FALSE)

```

### Define parameter transformation for two experimental conditions

```{r trafo}

# Get all parameters

innerpars <- getParameters(g*x)

# Identity transformation

trafo <- repar("x~x", x = innerpars)

# Set some initial value parameters

trafo <- repar("x~0", x = c("Compl", "Prod"), trafo)

# Explicit log-transform of all parameters

trafo <- repar("x~exp(x)", x = innerpars, trafo)

## Split transformation into two

trafo1 <- trafo2<- trafo

# Set the degradation rate in the first condition to 0

trafo1["k4"] <- "0"

# Generate parameter transformation functions

p <- NULL

p <- p + P(trafo1, condition = "noDegradation")

p <- p + P(trafo2, condition = "withDegradation")

```

### Initialize parameters and make prediction

```{r prediction}

# Initialize with randomly chosen parameters

set.seed(1)

outerpars <- getParameters(p)

pouter <- structure(rnorm(length(outerpars), -2, .5), names = outerpars)

times <- 0:100

plot((g*x*p)(times, pouter))

```

### Define data to be fitted by the model

```{r data}

data <- datalist(

noDegradation = data.frame(

name = c("product", "product", "product", "substrate", "substrate", "substrate"),

time = c(0, 25, 100, 0, 25, 100),

value = c(0.0025, 0.2012, 0.3080, 0.3372, 0.1662, 0.0166),

sigma = 0.02),

withDegradation = data.frame(

name = c("product", "product", "product", "substrate", "substrate", "substrate", "enzyme", "enzyme", "enzyme"),

time = c(0, 25, 100, 0, 25, 100, 0, 25, 100),

value = c(-0.0301, 0.1512, 0.2403, 0.3013, 0.1635, 0.0411, 0.4701, 0.2001, 0.0383),

sigma = 0.02)

)

timesD <- sort(unique(unlist(lapply(data, function(d) d$time))))

# Compare data to prediction

plot(data) + geom_line()

plot((g*x*p)(times, pouter), data)

```

### Define an objective function to be minimized and run minimization by `trust()`

```{r trust}

# Define prior values for parameters

prior <- structure(rep(0, length(pouter)), names = names(pouter))

# Set up objective function

obj <- normL2(data, g*x*p) + constraintL2(mu = prior, sigma = 10)

# Optimize the objective function

myfit <- trust(obj, pouter, rinit = 1, rmax = 10)

plot((g*x*p)(times, myfit$argument), data)

```

### Compute the profile likelihood to analyze parameter identifiability

```{r profiles, fig.height = 5}

# Compute the profile likelihood around the optimum

bestfit <- myfit$argument

profiles <- profile(obj, bestfit, names(bestfit), limits = c(-10, 10), cores = 4)

# Take a look at each parameter

plotProfile(profiles)

```

```{r, echo=FALSE}

# Cleaning step

dlls <- list.files(pattern = "\\.(so|dll)$")

for (d in dlls) try(dyn.unload(d), silent = TRUE)

files <- list.files(pattern = "^(enzymeKinetics|obsfn|expl_parfn_noDegradation|expl_parfn_withDegradation).*(c|cpp|o|dll|so)$")

unlink(files)

```