sim.dp {bootf2} | R Documentation |
Simulate Dissolution Profiles
Description
Function to simulate dissolution profiles based on mathematical models or multivariate normal distribution.
Usage
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)
Arguments
tp |
Numeric vector of time points for the dissolution profiles. See Details. |
dp , dp.cv |
Numeric vectors of the target mean dissolution profile
( |
model |
Character strings of model names. Currently only |
model.par |
A list with model parameters. If missing, a list of
random |
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 |
message |
Logical. If |
time.unit |
Character strings indicating the unit of time. It should
be either |
plot |
Logical. If |
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 |
Details
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 usingtp
,model
, andmodel.par
. The parameters defined bymodel.par
are considered as the population parameter; consequently, the calculateddp
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, todp
. The difference betweensim.mean
anddp
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, thenn.units
of individual dissolution profiles will be simulated using multivariate normal distribution. The mean of all simulated individual profiles,sim.mean
, will be identical todp
. In such case, it is recommended thatdp.cv
, the CV at time pointstp
, should also be provided by the user. Ifdp.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 thedp.cv
is provided or generated automatically, the CV calculated using all individual profiles will be equal todp.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 = 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
).
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.
Value
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., whendp
is missing. -
sim.dp
: A plot. Available only if the argumentplot
isTRUE
.
Examples
# 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)