Adapt your model • campsismod

Prerequisite

For this vignette, please load the campsismod package and load the minimalist model that we have created in the previous vignette.

library(campsismod)
model <- read.campsis("resources/minimalist_model/")

Define the concentration

A concentration is defined by the amount of drug present in the central compartment divided by the apparent volume of distribution.
Let’s now define the volume as a fixed constant of 100 in the model. This can be achieved as follows:

model <- model %>% add(Equation("V", "100"))
model

## [MAIN]
## K=THETA_K*exp(ETA_K) # Elimination constant
## V=100
## 
## [ODE]
## d/dt(A_CENTRAL)=-K*A_CENTRAL
## 
## [INIT]
## A_CENTRAL=1000 # Initial value
## 
## 
## THETA's:
##   name index value   fix
## 1    K     1  0.06 FALSE
## OMEGA's:
##   name index index2 value   fix type
## 1    K     1      1    15 FALSE  cv%
## SIGMA's:
## # A tibble: 0 × 0
## No variance-covariance matrix
## 
## Compartments:
## A_CENTRAL (CMT=1)

By default, this new equation has been appended to the model parameters, in the MAIN code record.

Now, we would like to define the plasma concentration. This can be done by adding an equation into the ODE code record. To do that, we’ll use the argument pos to define the position of this equation. It can be inserted, for instance, right after the ordinary differential equation:

model <- model %>% add(Equation("CP", "A_CENTRAL/V"), pos=Position(Ode("A_CENTRAL")))
model

## [MAIN]
## K=THETA_K*exp(ETA_K) # Elimination constant
## V=100
## 
## [ODE]
## d/dt(A_CENTRAL)=-K*A_CENTRAL
## CP=A_CENTRAL/V
## 
## [INIT]
## A_CENTRAL=1000 # Initial value
## 
## 
## THETA's:
##   name index value   fix
## 1    K     1  0.06 FALSE
## OMEGA's:
##   name index index2 value   fix type
## 1    K     1      1    15 FALSE  cv%
## SIGMA's:
## # A tibble: 0 × 0
## No variance-covariance matrix
## 
## Compartments:
## A_CENTRAL (CMT=1)

Define an error model

Say we want to add a proportional error model on the concentration with a coefficient of variation of 20%.
We therefore need to add a new parameter SIGMA into the model:

model <- model %>% add(Sigma("PROP", value=20, type="cv%"))

Then, we need to add a new ERROR code record with the appropriate equations:

error <- ErrorRecord()
error <- error %>% add(Equation("OBS_CP", "CP*(1 + EPS_PROP)"))
model <- model %>% add(error)
model

## [MAIN]
## K=THETA_K*exp(ETA_K) # Elimination constant
## V=100
## 
## [ODE]
## d/dt(A_CENTRAL)=-K*A_CENTRAL
## CP=A_CENTRAL/V
## 
## [INIT]
## A_CENTRAL=1000 # Initial value
## 
## [ERROR]
## OBS_CP=CP*(1 + EPS_PROP)
## 
## 
## THETA's:
##   name index value   fix
## 1    K     1  0.06 FALSE
## OMEGA's:
##   name index index2 value   fix type
## 1    K     1      1    15 FALSE  cv%
## SIGMA's:
##   name index index2 value   fix type
## 1 PROP     1      1    20 FALSE  cv%
## No variance-covariance matrix
## 
## Compartments:
## A_CENTRAL (CMT=1)

Simulate the observed concentration

Let’s now simulate a few individuals and show OBS_CP.

library(campsis)
dataset <- Dataset(3) %>% add(Observations(seq(0,24,by=3)))
results <- model %>% simulate(dataset=dataset, seed=0)
spaghettiPlot(results, "OBS_CP")

A couple of useful functions in action

We can check the existence of an equation (or any other type of model statement), by calling the function contains.

model %>% contains(Equation("CP"))

## [1] TRUE

model %>% contains(Ode("A_CENTRAL"))

## [1] TRUE

In the same way, we can retrieve any model statement easily using the function find:

model %>% find(Equation("CP"))

## CP=A_CENTRAL/V

model %>% find(Ode("A_CENTRAL"))

## d/dt(A_CENTRAL)=-K*A_CENTRAL

For instance, right-hand side formula of equation CP can be retrieved as follows:

(model %>% find(Equation("CP")))@rhs

## [1] "A_CENTRAL/V"

Any model statement may be replaced using the function replace:

model %>% replace(Equation("V", "50")) # Previous value of 100 is overridden

## [MAIN]
## K=THETA_K*exp(ETA_K) # Elimination constant
## V=50
## 
## [ODE]
## d/dt(A_CENTRAL)=-K*A_CENTRAL
## CP=A_CENTRAL/V
## 
## [INIT]
## A_CENTRAL=1000 # Initial value
## 
## [ERROR]
## OBS_CP=CP*(1 + EPS_PROP)
## 
## 
## THETA's:
##   name index value   fix
## 1    K     1  0.06 FALSE
## OMEGA's:
##   name index index2 value   fix type
## 1    K     1      1    15 FALSE  cv%
## SIGMA's:
##   name index index2 value   fix type
## 1 PROP     1      1    20 FALSE  cv%
## No variance-covariance matrix
## 
## Compartments:
## A_CENTRAL (CMT=1)

Finally, model statements can also be deleted forever (making the model broken in the following case):

model %>% delete(Equation("V"))

## [MAIN]
## K=THETA_K*exp(ETA_K) # Elimination constant
## 
## [ODE]
## d/dt(A_CENTRAL)=-K*A_CENTRAL
## CP=A_CENTRAL/V
## 
## [INIT]
## A_CENTRAL=1000 # Initial value
## 
## [ERROR]
## OBS_CP=CP*(1 + EPS_PROP)
## 
## 
## THETA's:
##   name index value   fix
## 1    K     1  0.06 FALSE
## OMEGA's:
##   name index index2 value   fix type
## 1    K     1      1    15 FALSE  cv%
## SIGMA's:
##   name index index2 value   fix type
## 1 PROP     1      1    20 FALSE  cv%
## No variance-covariance matrix
## 
## Compartments:
## A_CENTRAL (CMT=1)