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)Numeric vector of time points for the dissolution profiles. See Details.
Numeric vectors of the target mean dissolution profile
(dp) and its respective CV at time points tp (dp.cv). See Details.
Character strings of model names. Currently only "Weibull"
and "first-order" models are supported.
A list with model parameters. If missing, a list of
random model.par will be generated. See Details.
Integer seed value for reproducibility. If missing, a random seed will be generated for reproducibility purpose.
An integer indicating the number of individual profiles to be simulated.
Character strings indicating the product name of the simulated profiles. If missing, a random name with 3 letters and 4 digits will be generated.
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.
Logical. If TRUE, simulated profiles will always
increase with time. Only applicable when the approach based on
multivariate normal distribution is used. See Details.
Logical. If TRUE, basic information of the simulation
will be printed on screen.
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().
Logical. If TRUE, a a dissolution versus time plot will be
included in the output.
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().
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.
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.
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.
The first-order model is expressed as $$f_t = f_\mathrm{max} \left(1 - % e^{-k\left(t - t_\mathrm{lag}\right)}\right),$$ and the Weibull model was expressed either as $$f_t = f_\mathrm{max} \left(1 - % e^{-\left(\frac{t - t_\mathrm{lag}}{\mathrm{MDT}}% \right)^\beta}\right)$$ or $$f_t = f_\mathrm{max} \left(1 - % e^{-\frac{(t - t_\mathrm{lag})^\beta}{\alpha}}\right)$$ where \(f_\mathrm{max}\) is the maximum dissolution, \(\mathrm{MDT}\) is the mean dissolution time, \(t_\mathrm{lag}\) is the lag time, \(\alpha\) and \(\beta\) 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 \(\alpha = \mathrm{MDT}^\beta\).
Individual model parameter were simulated according to the exponential error model $$P_i = P_\mu e^{N\left(0, \sigma^2\right)}$$ where \(P_i\) and \(P_\mu\) denote the individual and population model parameters; \(N\left(0, \sigma^2\right)\) represents the normal distribution with mean zero and variance \(\sigma^2\) (\(\sigma = \mathrm{CV}/100\)).
model.parThe 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.
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.
# NOT RUN {
# 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)
# }
Run the code above in your browser using DataLab