sim.dp: Simulate Dissolution Profiles
In bootf2: Simulation and Comparison of Dissolution Profiles

Description Usage Arguments Details Value Examples

Function to simulate dissolution profiles based on mathematical models or multivariate normal distribution.

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)

`tp`	Numeric vector of time points for the dissolution profiles. See Details.
`dp, dp.cv`	Numeric vectors of the target mean dissolution profile (`dp`) and its respective CV at time points `tp` (`dp.cv`). See Details.
`model`	Character strings of model names. Currently only `"Weibull"` and `"first-order"` models are supported.
`model.par`	A list with model parameters. If missing, a list of random `model.par` will be generated. See Details.
`seed`	Integer seed value for reproducibility. If missing, a random seed will be generated for reproducibility purpose.
`n.units`	An integer indicating the number of individual profiles to be simulated.
`product`	Character strings indicating the product name of the simulated profiles. If missing, a random name with 3 letters and 4 digits will be generated.
`max.disso`	Numeric value for the maximum possible dissolution. In theory, the maximum dissolution is 100%, but in practice, it is not uncommon to see higher values, such as 105%, or much lower values, especially for products with low solubility.
`ascending`	Logical. If `TRUE`, simulated profiles will always increase with time. Only applicable when the approach based on multivariate normal distribution is used. See Details.
`message`	Logical. If `TRUE`, basic information of the simulation will be printed on screen.
`time.unit`	Character strings indicating the unit of time. It should be either `"min"` for minute or `"h"` for hour. It is mainly used for and making plot and generating `model.par` and `dp.cv` when they are not provided by the user. @seealso `calcf2()`.
`plot`	Logical. If `TRUE`, a a dissolution versus time plot will be included in the output.
`plot.max.unit`	Integer. If the number of individual units is no more than this value, the mean and all individual profiles will be plotted; otherwise, individual profiles will be represented by boxplots at each time point. Therefore, to avoid overplotting, this value should not be too large. @seealso `calcf2()`.

Simulation approaches

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. Additional details on the mathematical models are given below.
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. Additional details on this approach are given below.

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 mathematical models

The first-order model is expressed as

f(t) = fmax (1 - exp(-k(t - tlag))),

and the Weibull model was expressed either as

f(t) = fmax(1 - exp(-((t - tlag)/MDT)^β))

f(t) = fmax (1 - exp(-(((t - tlag)^β)/α))),

where fmax is the maximum dissolution, MDT is the mean dissolution time, tlag is the lag time, α and β are the scale and shape factor in Weibull function, and k is the rate constant in the first-order model. Obviously, The two Weibull models are mathematically equivalent by letting α = MDT^β.

Individual model parameter were simulated according to the exponential error model

P(i) = P(μ) exp(N(0, σ^2)),

where P(i) and P(μ) denote the individual and population model parameters; N(0, σ^2) represents the normal distribution with mean zero and variance σ^2 (σ = CV/100).

How to supply `model.par`

The argument model.par should be supplied as a named list of model parameters as explained above, and their respective CV for simulation of individual parameters. Therefore, for the first-order model, three pairs of parameters should be specified: fmax/fmax.cv, k/k.cv, and tlag/tlag.cv; and for Weibull model, four pairs: fmax/fmax.cv, tlag/tlag.cv, beta/beta.cv, and either alpha/alpha.cv or mdt/mdt.cv, depending on the mathematical formula used. CV should be given in percentage, e.g., if CV for beta is 30%, it should be specified as beta.cv = 30, not beta.cv = 0.3. The order of the parameters does not matter, but the name of the parameters must be given exactly same as described above. For example:

model.par = list(fmax = 100, fmax.cv = 5, k = 0.6, k.cv = 25, tlag = 0, tlag.cv = 0) for the first-order model.
model.par = list(fmax = 100, fmax.cv = 5, tlag = 5, tlag.cv = 10, mdt = 15, mdt.cv = 20, beta = 1.5, beta.cv = 5), or
model.par = list(fmax = 100, fmax.cv = 5, tlag = 5, tlag.cv = 10, alpha = 60, alpha.cv = 20, beta = 1.5, beta.cv = 5) for Weibull models.

Notes on multivariate normal distribution approach

When this approach is used, depending on dp/dp.cv, there is no guarantee that all individual profiles increase with time; near the end of the time points, some individual profiles may decrease, especially when the dissolution is close to the plateau and the variability is high. This can happen to real life data, especially for those products with drug substances that are unstable in dissolution medium. To force increase for all profiles, set ascending = TRUE. Depending on the data, the program may take long time to run, or may even fail.

A list of 3 to 5 components:

sim.summary: A data frame with summary statistics of all individual dissolution profiles.
sim.disso: A data frame with all individual dissolution profiles.
sim.info: A data frame with information of the simulation such as the seed number and number of individual profiles. If modelling approach is used, the data frame will contain model parameters as well.
model.par.ind: A data frame of all individual model parameters that were used for the simulation of individual dissolution profiles. Available only if the modelling approach is used, i.e., when dp is missing.
sim.dp: A plot. Available only if the argument plot is TRUE.

# time points
tp <- c(5, 10, 15, 20, 30, 45, 60)

# using all default values to simulate profiles
d1 <- sim.dp(tp, plot = FALSE)

# summary statistics
d1$sim.summary

# individual profiles
d1$sim.disso

# simulation information
d1$sim.info

#individual model parameters
d1$mod.par.ind

# using Weibull model to simulate 100 profiles without lag time
mod.par <- list(fmax = 100, fmax.cv = 2, tlag = 0, tlag.cv = 0,
                mdt = 20, mdt.cv = 5, beta = 2.2, beta.cv = 5)
d2 <- sim.dp(tp, model.par = mod.par, seed = 100, n.units = 100L,
             plot = FALSE)

# using multivariate normal distribution approach
# target mean profile with same length as tp
dp <- c(39, 56, 67, 74, 83, 90, 94)

# CV% at each time point
dp.cv <- c(19, 15, 10, 8, 8, 5, 3)

# simulation
d3 <- sim.dp(tp, dp = dp, dp.cv = dp.cv, seed = 1234, plot = FALSE)