--- title: "Simulation of Dissolution Profiles" author: "Zhengguo XU" date: "`r format(Sys.Date())`" output: rmarkdown::html_vignette: keep_md: true toc: true toc_depth: 3 fig_width: 7 fig_height: 4.5 highlight: "tango" bibliography: ref.bib notes-after-punctuation: false link-citations: yes csl: ref.csl vignette: > %\VignetteIndexEntry{Simulation of Dissolution Profiles} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup0, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#" ) ``` ```{r setup, echo=FALSE} library(bootf2) ``` ## Introduction Dissolution profiles are simulated either by mathematical models, or by multivariate normal distribution based on function `mvrnorm()` from R package `MASS`. ### Mathematical Models Two models are implemented at the moment: *first-order* model and *Weibull* model (the default). The first-order model is expressed as $$f_t = f_\mathrm{max}\left(1 - e^{-k\cdot (t - t_\mathrm{lag})}\right),$$ and the Weibull model either as $$f_t = f_\mathrm{max}\left(1 - e^{-\left(\frac{t - t_\mathrm{lag}} {\mathrm{MDT}}\right)^\beta}\right),$$ or as $$f_t = f_\mathrm{max}\left(1 - e^{-\frac{\left(t - t_\mathrm{lag}\right)^\beta} {\alpha}}\right).$$ Obviously, the two expressions of Weibull model are mathematically equivalent with $\alpha = \mathrm{MDT}^\beta$. Parameter $f_\mathrm{max}$ is the maximum dissolution. In theory, it is 100, but in practice, it might be slightly high than 100, or much less if the dissolution is not complete. For immediate-release formulation, the lag time $t_\mathrm{lag}$ is typically zero, but not necessary so for the extended release formulation. Parameter $k$ is the first-order rate constant, and $\alpha$ and $\beta$ are scale and shape parameters, and $\mathrm{MDT}$ is the mean dissolution time. ### Simulation based on mathematical models Let $P_\mu$ be the set of *population* model parameters, i.e., $P_\mu = \{f_\mathrm{max}, k, t_\mathrm{lag}\}$ for the first-order model and $P_\mu = \{f_\mathrm{max}, t_\mathrm{lag}, \mathrm{MDT}, \beta\}$ or $P_\mu = \{f_\mathrm{max}, t_\mathrm{lag}, \alpha, \beta\}$ for Weibull model, given any time points `tp`, the dissolution profile $f_t$ will be calculated according to the equations above with parameters $P_\mu$. The generated profile is considered as *population* dissolution profile `dp`. Individual dissolution profiles are generated with the same equation, but with individual model parameters $P_i$, which is simulated according to exponential error model $$P_i = P_\mu \cdot e^{N\left(0,\, \sigma^2\right)},$$ where $N\left(0,\, \sigma^2\right)$ represents the normal distribution with mean zero and variance $\sigma^2$ ($\sigma = \mathrm{CV}/100$) ### Simulation based on multivariate normal distribution Note that the *population* dissolution profile `dp` in the previous section is generated by equation with *population* model parameters $P_\mu$. Sometime we might have a mean dissolution profile `dp`, which is considered as population profile, and we want to simulate many individual profiles such that the calculated mean profile is equal to `dp`. In such case, multivariate normal distribution approach will be used. Given any `tp`, `dp`, and the required number of individual units, `n.units`, the function will simulate individual profiles fulfil the above mentioned requirements. ### Recommendation *Depending on the variability and the target profile, the run time of simulation based on multivariate normal distribution might be very long, sometimes it might even fail to finish. So in general, the modelling approach is recommended.* ## Usage The complete list of arguments of the function is as follows: ```{r simdp-code0, eval=FALSE} sim.dp(tp, dp, dp.cv, model = c("Weibull", "first-order"), model.par = NULL, seed = NULL, n.units = 12L, product, max.disso = 100, ascending = FALSE, message = FALSE, time.unit = c("min", "h"), plot = TRUE, plot.max.unit = 36L) ``` The approach used to simulate individual dissolution profiles depends on if the target mean dissolution profile dp is provided by the user or not. - If `dp` is not provided, then it will be calculated using `tp`, `model`, and `model.par`. The parameters defined by `model.par` are considered as the *population parameter*; consequently, the calculated `dp` is considered as the *targeted population profile*. In addition, `n.units` sets of *individual model parameters* will be simulated based on exponential error model, and *individual dissolution profiles* will be generated using those individual parameters. The mean of all simulated individual profiles, `sim.mean`, included in one of the outputs data frames, `sim.summary`, will be *similar, but not identical, to `dp`*. The difference between `sim.mean` and `dp` is determined by the number of units and the variability of the model parameters. In general, the larger the number of units and the lower of the variability, the smaller the difference. - If `dp` is provided, then `n.units` of individual dissolution profiles will be simulated using multivariate normal distribution. The mean of all simulated individual profiles, `sim.mean`, will be *identical to `dp`*. In such case, it is recommended that `dp.cv`, the CV at time points `tp`, should also be provided by the user. If `dp.cv` is not provided, it will be generated automatically such that the CV is between 10% and 20% for time points up to 10 min, between 1% and 3% for time points where the dissolution is more than 95%, between 3% and 5% for time points where the dissolution is between 90% and 95%, and between 5% and 10% for the rest of time points. Whether the `dp.cv` is provided or generated automatically, the CV calculated using all individual profiles will be equal to `dp.cv`. ### Minimum required arguments that must be provided by the user If `dp` is provided by the user, logically, `tp` must also be provided, and the approach based on multivariate normal distribution is used, as explained above. If `dp` is not provided, `tp` could be missing, i.e., a simple function call such as `sim.dp()` will simulate dissolution profiles. In such case, a default `tp` will be generated depending on the `time.unit`: if the `time.unit` is `"min"`, then `tp` would be `c(5, 10, 15, 20, 30, 45, 60)`; otherwise the `tp` would be `c(1, 2, 3, 4, 6, 8, 10, 12)`. The rest arguments are either optional or required by the function but default values will be used. ### Notes on function arguments 1. Model parameter `model.par` has to be specified as *named list*, e.g., `list(fmax = 100, fmax.cv = 3, tlag = 0, tlag.cv = 0, k = 0.9, k.cv = 30)`. The parameter `fmax/k/tlag` are used to generate `dp`, and the corresponding `.cv` parts are used to generate individual model parameters. If `model.par` is missing, it will be generated randomly. 1. When multivariate normal distribution approach is used, depending on the target profile and the variability, sometimes the simulated individual profiles will *decrease* in the end. This could also happen to the real-life profiles, especially for those products that are unstable in the dissolution media. To force the profiles *always increase* with time, set option `ascending = TRUE`. However, in that case, it is possible that the function takes long time to run or even fails. 1. Parameter `product` is not really necessary for the simulation. so if missing, it will be generated automatically with 3 letters and 4 digits. It might be useful in the situations such as many simulation will be run and output are pooled together for analysis. 1. Read the function manual by `help("sim.dp")` for more details on each argument. ## Examples ### Simple case For the most basic use, the function can be run without any user provided arguments, e.g., `sim.dp()`. In such case, 12 units of individual dissolution profiles will be simulated using Weibull model with a typical sampling time points of 5, 10, 15, 20, 30, 45, and 60 min. A `seed` number will be randomly generated, if not provided by the user, and included in the output for reproducibility purpose. ```{r simple01-dat} # simulation tmp1 <- sim.dp(seed = 1234) ``` The output is a list of 5 components: 1. `sim.summary`: A *data frame* with summary statistics of all individual dissolution profiles. ```{r simple01-summary, error=TRUE} tmp1$sim.summary ``` `dp` is the population mean profile obtained with the parameters in `sim.info`, as explained in previous section. The rest columns with prefix `sim` are basic descriptive statistics calculated from all simulated individual profiles. `sim.qt05`, `sim.qt25`, ..., are 5%, 25%, .., quantiles. 2. `sim.disso`: A *data frame* with all individual dissolution profiles. ```{r simple01-disso, error=TRUE} tmp1$sim.disso ``` 3. `sim.info`: A *data frame* with information of the simulation. ```{r simple01-info, error=TRUE} tmp1$sim.info ``` 4. `model.par.ind`: A *data frame* of individual model parameters that are used to simulate the individual dissolution profiles if mathematical models are chosen for the simulation. ```{r simple01-par, error=TRUE} tmp1$model.par.ind ``` 5. `sim.plot`: A plot since the default argument `plot` is set to be `TRUE`. ```{r simple01-plot, error=TRUE} tmp1$sim.plot ``` When the number of individual units is not large, all individual profiles will be plotted, together with the population profile in green (as target profile), and mean profile in blue calculated with simulated individual profiles. When the number is small, there will be visible difference between the mean and the target profile due to random error. When the number of units increase, the difference will become smaller. The argument `plot.max.unit` control how individual profile will be represented in the plot. When the actual number of units is greater than the value of `plot.max.unit`, the individual profile will be represented as boxplots at each time points, as shown below. ```{r simple01-boxplot, error=TRUE} # default plot.max.unit = 36 sim.dp(n.units = 100)$sim.plot ``` ### Example with model parameters To obtain better controlled simulation, model parameters should be provided. CV should be specified in *percentage*. ```{r model01-dat, error=TRUE} fo.par <- list(fmax = 100, fmax.cv = 3, k = 0.1, k.cv = 20, tlag = 0, tlag.cv = 0) fo.dat <- sim.dp(model = "first-order", model.par = fo.par, seed = 123) fo.dat$sim.plot fo.dat$sim.summary ``` Similarly, for Weibull model, the model parameter should be provided like the example below. The order of the parameter does not matter, but the parameter names have to be specified exactly as the following: ```{r model01-par, eval=FALSE} # with alpha = xx and alpha.cv = yy to replace beta/beta.cv if alternative # expression of Weibull model is used. mod.par <- list(fmax = 100, fmax.cv = 3, tlag = 0, tlag.cv = 0, mdt = 20, mdt.cv = 25, beta = 2, beta.cv = 30) ``` ### Example with multivariate normal distribution ```{r mvn01-dat, error=TRUE} # target mean profile dp <- c(39, 56, 67, 74, 83, 90, 94) # CV at each time points dp.cv <- c(19, 15, 10, 8, 8, 5, 3) mvn.dat <- sim.dp(tp, dp = dp, dp.cv = dp.cv, seed = 1234) mvn.dat$sim.summary ``` Notice that the the mean and CV of the simulated individual profile ( `sim.mean` and `sim.cv`) are equal to the target profile `dp` and `dp.cv`. The plot look like this: ```{r mvn01-plot, error=TRUE} mvn.dat$sim.plot ``` Another example with missing `dp.cv`. ```{r mvn02-dat, error=TRUE} mvn.dat2 <- sim.dp(tp, dp = dp, seed = 123) mvn.dat2$sim.summary ``` Notice that the `dp.cv` were generated automatically as explained in [Notes on function arguments].