---
title: "Model"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r cran, eval=FALSE}
install.packages("PKPDsim")
```

# Introduction

`new_ode_model` is the function that creates a new ODE model that can be used in the `sim()` command. It defines the ODE system and sets some attributes for the model. The model can be specified in three different ways:

- `model`: a string that references a model from the library included in `PKPDsim`. Examples in the current library are e.g. `pk_1cmt_oral`, `pk_2cmt_iv`. To show the available models, run `new_ode_model()` without any arguments.
- `code`: using code specifying derivatives for ODE specified in *pseudo-R* code
- `file`: similar to `code`, but reads the code from a file

# Model from library

For example, a 1-compartment oral PK model can be obtained using:

```{r load-lib, echo=FALSE}
library("PKPDsim")
```

```{r new-model}
pk1 <- new_ode_model(model = "pk_1cmt_oral")
```

Run the `new_ode_model()` function without arguments to see the currently available models:

```{r available-models, error=TRUE}
new_ode_model()
```

# Custom model from code

The custom model needs to be specified as a string or text block:

```{r custom-model}
pk1 <- new_ode_model(code = "
  dAdt[1] = -KA * A[1]
  dAdt[2] = +KA * A[1] -(CL/V) * A[2]
")
```

The input code should adhere to the follow rules:

- the derivatives in the ODE system are defined using `dAdt`omputate
- array indices for the derivatives and compartments are indicated with `[ ]`. 
Compartments indices can start at either 0 or 1. If the latter, all indices will
be reduced by 1 in the translation to C++ code.
- equations are defined using `=`
- power functions need to be written out as `pow(base,exp)`.
- force any numbers to be interpreted as `real`, to avoid them being interpreted
as `integer`. So if an equation involves the real number `3`, it is usually safer
to write this as `3.0` in code.

The input code is translated into a C++ function. You can check that the model compiled correctly by typing the model name on the R command line, which prints the model information:

```{r print-custom-model}
pk1
```

If you're interested, you can also output the actual C++ function that is compiled by specifying the `cpp_show_code=TRUE` argument to the `new_ode_model()` function.

# More custom model options

You can introduce new variables in your code, but you will have to define them using `declare_variables` argument too:

```{r custom-model-variables}
pk1 <- new_ode_model(code = "
  KEL = CL/V
  dAdt[1] = -KA * A[1]
  dAdt[2] = +KA * A[1] -KEL * A[2]
", declare_variables = c("KEL"))
```

Also, when you want to use covariates in your ODE system (more info on how to define covariates is in the *Covariates* vignette), you will have to define them, both in the code and in the function call:

```{r custom-model-covariates}
pk1 <- new_ode_model(code = "
  CLi = WT/70
  KEL = CLi/V
  dAdt[1] = -KA * A[1]
  dAdt[2] = +KA * A[1] -(CL*(WT/70)/V) * A[2]
", declare_variables = c("KEL", "CLi"), covariates = c("WT"))
```

One exception to the input code syntax is the definition of power functions. `PKPDsim` does not translate those from the *pseudo-R* code to valid C++ syntax automatically. C/C++ does not use the `^` to indicate power functions, but uses the `pow(value, base)` function instead, so for example an allometric PK model should be written as:

```{r custom-model-power-functions}
pk1 <- new_ode_model(code = "
  CLi = CL * pow((WT/70), 0.75)
  dAdt[1] = -KA * A[1]
  dAdt[2] = +KA * A[1] -(CLi/V) * A[2]
", declare_variables = c("CLi"))
```

## Dosing / bioavailability

The default dosing compartment and bioavailability can be specified using the `dose` argument. By default, the dose will go into compartment `1`, with a bioavailability of `1`. The `bioav` element in the list can be either a number or a character string referring a parameter.

```{r bioav}
pk1 <- new_ode_model(code = "
    dAdt[1] = -KA * A[1]
    dAdt[2] = +KA * A[1] -(CL/V) * A[2]
  ",
  dose = list(cmt = 1, bioav = "F1"),
  parameters = list(KA = 1, CL = 5, V = 50, F1 = 0.7)
)
```

Bioavailability can also be used for dosing based on mg/kg, since that is not supported in  `new_regimen()`. The way to implement this is by scaling the dose by the "weight" covariate using the bioavailability:

```{r scale}
mod <- new_ode_model(code = "
    dAdt[1] = -(CL/V)*A[1];
  ",
  dose = list(cmt = 1, bioav = "WT"),
  obs = list(cmt = 1, scale = "V"),
  covariates = list("WT" = new_covariate(value = 70))
)
```

## Observations

The observation compartment can be set by specifying a list to the `obs` argument, with either the elements `cmt` and `scale`, or `variable`.

```{r obs}
pk1 <- new_ode_model(code = "
    dAdt[1] = -KA * A[1]
    dAdt[2] = +KA * A[1] -(CL/V) * A[2]
  ", 
  obs = list(cmt = 2, scale = "V")
)
```

The `scale` can be either a parameter or a number, the `cmt` can only be a number.

*Note that the variables specified inside the differential equation block are not available as scaling parameters. E.g. for allometry you will have to redefine the scaled volume as follows:*

```{r scale-2}
pk1 <- new_ode_model(code = "
    Vi = V * (WT/70)
    dAdt[1] = -KA * A[1]
    dAdt[2] = +KA * A[1] -(CL/Vi) * A[2]
  ", 
  obs = list(cmt = 2, scale = "V * (WT/70)")
)
```

Or define the observation using a variable:

```{r variable}
pk1 <- new_ode_model(code = "
    dAdt[1] = -KA * A[1]
    dAdt[2] = +KA * A[1] -(CL/V) * A[2]
    CONC = A[2]
  ", 
  obs = list(variable = "CONC"),
  declare_variables = "CONC"
)
```


# Custom model from file

Using the `file=` argument, the model code is read from the specified files. This is just a convenience function, i.e. it allows you to separate models from R code more easily.

```{r model-from-file, eval=FALSE}
pk1 <- new_ode_model(
  file = "pk_1cmt_oral_nonlin_v1.txt",
  declare_variables = c("KEL", "CLi")
)
```