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/")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 same
## 1 K 1 1 15 FALSE cv% NA
## SIGMA's:
## # A tibble: 0 x 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 same
## 1 K 1 1 15 FALSE cv% NA
## SIGMA's:
## # A tibble: 0 x 0
## No variance-covariance matrix
##
## Compartments:
## A_CENTRAL (CMT=1)
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 same
## 1 K 1 1 15 FALSE cv% NA
## SIGMA's:
## name index index2 value fix type
## 1 PROP 1 1 20 FALSE cv%
## No variance-covariance matrix
##
## Compartments:
## A_CENTRAL (CMT=1)
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")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 same
## 1 K 1 1 15 FALSE cv% NA
## 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 same
## 1 K 1 1 15 FALSE cv% NA
## SIGMA's:
## name index index2 value fix type
## 1 PROP 1 1 20 FALSE cv%
## No variance-covariance matrix
##
## Compartments:
## A_CENTRAL (CMT=1)