SODF2: IV Dissolution similarity of tablets made at two manufacturing sites

Author

Dave LeBlond, Consultant in CMC Statistical Studies, Wadsworth, IL 60083

Published

September 2, 2020

Introduction

In this second example, we compare the dissolution profiles of a tablet product (SODF2) made at two different manufacturing sites, referred to as sites 1 and 2. We treat site 1 as the approved (or reference) site and site 2 as an unapproved (or test) site. We illustrate Bayesian model dependent tests of similarity based on both $F_{2}$ and multi-parameter space similarity regions.

Data description

We index the SODF2 manufacturing site with $m$ = 1 or 2. $N L 1$ = 8 lots were obtained from manufacturer 1 and $N L 2$ = 5 lots from manufacturer 2. Twelve units (tablets) from each lot were tested ( $N U 1$ = 96 and $N U 2$ = 60). Percent dissolution measurements were obtained from each unit at 1, 2, 4, 8, 16, 32, and 45 minutes. Thus, there were $N 1$ = 8x12x7 = 672 and $N 2$ = 5x12x7 = 420 total observations from these sites. The vectors $l o t m_{n}$ , $u n i t m_{n}$ , $y m_{n}$ , and $x m_{n}$ identify the unique lot, unique unit, dissolution time (0 to 45), and observed % dissolution, respectively, for observation $n$ from site $m$ .

Model description

Some key assumptions of this model

Conditional independence of sampling and prior events (i.e., all “~”s are iid).
Homogeneity of error variance, $σ^{2}$ as a function of time and mean % dissolution. This is an important assumption that is partially justified in this example. However, usually variance is lower near zero and 100% dissolution and higher near 50% dissolution.
Lot and unit (tablet) effects can be approximated as linear effects on the 3 underlying Weibull coefficients ( $M$ , $\log (T)$ , and $\log (b)$ ).
Covariance of lot effects are negligible. Actually this assumption has little justification. however, with only 6 lots the covariances are not well determined. So we assume they are zero.
The dissolution profile of all units (tablets) can be adequately described by a 3-parameter Weibull model. The Weibull model does have a good success history. However, there are other models that may be appropriate in various applications.

Hierarchical likelihood sampling model

The model is identical to that used in the accompanying SODF1 example, except that data from each manufacturing site are included. The data from each manufacturing site are fitted independently except that the same analytical standard deviation, $σ$ , is assumed for both since the dissolution testing laboratory was the same for both products.

The variable names have been necessarily modified from the SODF1 example and are slightly different than those described in BMIPR Chapter 23. The variable names include an “ $m$ ” suffix which identifies the manufacturing site. The model for each site is given below.

${y m}_{n} \overset{i i d}{\sim} N ({w b m}_{n}, σ^{2})$

${w b m}_{n} = {θ m}_{1, u n i {t m}_{n}} \cdot [1 - \exp {{(- \frac{{x m}_{n}}{\exp ({θ m}_{2, u n i {t m}_{n}})})}^{\exp ({θ m}_{3, u n i {t m}_{n}})}}]$

${(\begin{matrix} M m \\ \log (T m) \\ \log (b m) \end{matrix})}_{u n i t_{n}} = {(\begin{matrix} {θ m}_{1} \\ {θ m}_{2} \\ {θ m}_{3} \end{matrix})}_{u n i t_{n}} = (\begin{matrix} {μ m}_{1} \\ {μ m}_{2} \\ {μ m}_{3} \end{matrix}) + {(\begin{matrix} {β l m}_{1} \\ {β l m}_{2} \\ {β l m}_{3} \end{matrix})}_{l o {t m}_{n}} + {(\begin{matrix} {β u m}_{1} \\ {β u m}_{2} \\ {β u m}_{3} \end{matrix})}_{u n i {t m}_{n}}$

${θ m}_{u n i {t m}_{n}} = μ m + {β lm}_{l o {t m}_{n}} + {β um}_{u n i {t m}_{n}}$

${β lm}_{i} \sim M V N_{3} (0, [\begin{matrix} {σ l m}_{1} & 0 & 0 \\ 0 & {σ l m}_{2} & 0 \\ 0 & 0 & {σ l m}_{3} \end{matrix}]), i = 1 : N L m$

${β um}_{i} \sim M V N_{3} (0, Σ_{u n i t m}), i = 1 : N U m$

$Σ u n i t m = [\begin{matrix} σ u m_{1} \\ σ u m_{2} \\ σ u m_{3} \end{matrix}] \cdot [\begin{matrix} 1 & C m_{12} & C m_{13} \\ C m_{12} & 1 & C m_{23} \\ C m_{13} & C m_{23} & 1 \end{matrix}] \cdot [\begin{matrix} σ u m_{1} & σ u m_{2} & σ u m_{3} \end{matrix}] = σ u m \cdot Cm \cdot σ u m^{T}$

Prior model

Initially, more informative priors were employed, but these were gradually relaxed. The vague priors seemed to have little impact on the posterior estimates. Ultimately it was found that essentially “flat” priors gave essentially the same estimates as vague priors. It is possible that informative priors would be more appropriate scientifically and also might lead to better convergence properties.

$σ \sim U n i f (0, + \infty)$

${μ m}_{i} \sim U n i f (- \infty, + \infty), i = 1 : 3$

${σ l m}_{i} \sim U n i f (0, + \infty), i = 1 : 3$

${σ u m}_{i} \sim U n i f (0, + \infty), i = 1 : 3$

$Cm \sim L K J c o r (1)$

Stan model code

# Create model for Stan
modelString = "
data {
    int<lower=1> K;

int<lower=0> N1; 
int<lower=1> NL1;
int<lower=1> NU1;
int<lower=1,upper=NL1> lot1[N1];
int<lower=1,upper=NU1> unit1[N1];
vector[N1] y1;
vector[N1] x1;

int<lower=0> N2;
int<lower=1> NL2;
int<lower=1> NU2;
int<lower=1,upper=NL2> lot2[N2];
int<lower=1,upper=NU2> unit2[N2];
vector[N2] y2;
vector[N2] x2;
}

transformed data {
  vector[K] zero;
  zero = rep_vector(0., K);
}
parameters {
real<lower=0> sigma;

cholesky_factor_corr[K] L1;
vector<lower=0>[K] sig_l1;
vector<lower=0>[K] sig_u1;
vector[K] beta_l1[NL1]; 
vector[K] beta_u1[NU1]; 
vector[K] mu1;

cholesky_factor_corr[K] L2;
vector<lower=0>[K] sig_l2;
vector<lower=0>[K] sig_u2;
vector[K] beta_l2[NL2]; 
vector[K] beta_u2[NU2]; 
vector[K] mu2;
}

transformed parameters {
vector[N1] wb1;
vector[K] theta1;
vector[N2] wb2;
vector[K] theta2;
matrix[K,K] C1 = L1 * L1';
matrix[K,K] C2 = L2 * L2';
matrix[K,K] Sigma_unit1 = quad_form_diag(C1, sig_u1);
matrix[K,K] Sigma_unit2 = quad_form_diag(C2, sig_u2);
for (n in 1:N1){
theta1 = mu1 + beta_l1[lot1[n]] + beta_u1[unit1[n]];
wb1[n] = theta1[1]*(1 - exp(-(x1[n]/exp(theta1[2]))^exp(theta1[3])));
}
for (n in 1:N2){
theta2 = mu2 + beta_l2[lot2[n]] + beta_u2[unit2[n]];
wb2[n] = theta2[1]*(1 - exp(-(x2[n]/exp(theta2[2]))^exp(theta2[3])));
}
}

model {
// Priors
//sigma ~ normal(0, 1);
//mu1[1] ~ normal(100, 320);
//mu1[2] ~ normal(0, 64);
//mu1[3] ~ normal(0, 64);
//sig_l1 ~ normal(0, 64);
//sig_u1 ~ normal(0, 64);
L1 ~ lkj_corr_cholesky(1);
//mu2[1] ~ normal(100, 320);
//mu2[2] ~ normal(0, 64);
//mu2[3] ~ normal(0, 64);
//sig_l2 ~ normal(0, 64);
//sig_u2 ~ normal(0, 64);
L2 ~ lkj_corr_cholesky(1);

// Likelihood
for (i in 1:NL1) beta_l1[i] ~ normal(0, sig_l1); 
beta_u1 ~ multi_normal_cholesky(zero, diag_pre_multiply(sig_u1,L1));
y1 ~ normal(wb1, sigma);
for (i in 1:NL2) beta_l2[i] ~ normal(0, sig_l2); 
beta_u2 ~ multi_normal_cholesky(zero, diag_pre_multiply(sig_u2,L2));
y2 ~ normal(wb2, sigma);
}
" # close quote for modelString

MCMC execution time

The code is set up to request 16,000 MCMC draws ( 4 chains of 4,000 draws each after warm up). This can take as long as 2 hours on a 4 core machine. Considerably fewer draws may be adequate. Some parameters appear poorly estimated and result in disappointing n_eff and Rhat (e.g., $σ u 1_{3}$ ). However, the model is likely somewhat over-specified. However I feel that convergence was adequate with reliable estimation and prediction uncertainty.

I recommend modifying the code to execute the stan() function and then save the StanFit object to disk. Then in subsequent executions (if needed) to adjust the code to read the StanFit object from disk.

As always, patience, persistence, perseverance, and due-diligence are in order to assure MCMC convergence. I have found it is helpful to purchase a rubber ducky and keep it nearby as you work. It can be a source of comfort and inspiration in times of stress.

[1] "2022-09-23 16:39:25 BST"

[1] "2022-09-23 16:39:30 BST"

Posterior statistical summary

[1] "Table 23.5: Posterior summary for two manufacturers of SODF2."

                 mean      se_mean           sd         2.5%         25%         50%         75%        97.5%       n_eff     Rhat
sigma      2.28017066 7.025449e-04 6.206215e-02  2.163255698  2.23742922  2.27871696  2.32099255   2.40476285  7803.78790 1.001072
sig_l1[1]  2.41475494 1.065102e-02 9.251326e-01  1.280218921  1.80139038  2.21653073  2.77658711   4.73453700  7544.41175 1.000619
sig_l1[2]  0.08747982 7.597702e-04 4.146739e-02  0.028806398  0.06055267  0.08070320  0.10552064   0.18841095  2978.85166 1.000249
sig_l1[3]  0.04950760 2.180364e-04 2.108883e-02  0.022171220  0.03542551  0.04539870  0.05846073   0.10212481  9355.07447 1.000367
sig_l2[1]  5.44518040 1.504359e-01 4.471154e+00  2.089350362  3.27963694  4.40631096  6.19202696  14.97376988   883.35754 1.004873
sig_l2[2]  0.39509036 8.253042e-03 3.268176e-01  0.152495805  0.23531278  0.31329146  0.44060345   1.13635348  1568.13228 1.003097
sig_l2[3]  0.06345259 1.007974e-03 5.697027e-02  0.004187357  0.03025603  0.05088474  0.07904362   0.20203897  3194.46283 1.001520
sig_u1[1]  1.35736328 4.726255e-03 2.195613e-01  0.939765040  1.20944367  1.35082682  1.49913433   1.79824680  2158.12957 1.002936
sig_u1[2]  0.15717111 1.324626e-04 1.256546e-02  0.134534303  0.14826410  0.15647380  0.16523001   0.18367945  8998.50843 1.000126
sig_u1[3]  0.02102480 1.281161e-03 1.136480e-02  0.001411619  0.01254504  0.02087509  0.02897995   0.04398608    78.68937 1.046081
sig_u2[1]  2.57945803 2.820568e-03 3.209238e-01  2.017568358  2.35490036  2.55983303  2.78010014   3.27259280 12945.85528 1.000201
sig_u2[2]  0.14196343 1.380837e-04 1.534317e-02  0.114786489  0.13130494  0.14076951  0.15153301   0.17492592 12346.54103 1.000104
sig_u2[3]  0.08866669 2.035069e-04 1.537214e-02  0.059796450  0.07825088  0.08823686  0.09844477   0.12028818  5705.72502 1.000291
C1[1,1]    1.00000000          NaN 0.000000e+00  1.000000000  1.00000000  1.00000000  1.00000000   1.00000000         NaN      NaN
C1[1,2]    0.70624704 1.913693e-03 1.016140e-01  0.487182324  0.64336911  0.71382389  0.77706867   0.88431367  2819.43746 1.003651
C1[1,3]    0.14853170 9.950780e-03 3.834284e-01 -0.599322329 -0.12749011  0.15211476  0.43214202   0.85330828  1484.75306 1.002639
C1[2,1]    0.70624704 1.913693e-03 1.016140e-01  0.487182324  0.64336911  0.71382389  0.77706867   0.88431367  2819.43746 1.003651
C1[2,2]    1.00000000 4.578939e-19 5.131565e-17  1.000000000  1.00000000  1.00000000  1.00000000   1.00000000 12559.43236 0.999750
C1[2,3]    0.45586639 7.418190e-03 3.243305e-01 -0.358520931  0.28393129  0.50176243  0.69114552   0.93554032  1911.52397 1.001038
C1[3,1]    0.14853170 9.950780e-03 3.834284e-01 -0.599322329 -0.12749011  0.15211476  0.43214202   0.85330828  1484.75306 1.002639
C1[3,2]    0.45586639 7.418190e-03 3.243305e-01 -0.358520931  0.28393129  0.50176243  0.69114552   0.93554032  1911.52397 1.001038
C1[3,3]    1.00000000 4.784209e-19 5.283211e-17  1.000000000  1.00000000  1.00000000  1.00000000   1.00000000 12194.82910 0.999750
C2[1,1]    1.00000000          NaN 0.000000e+00  1.000000000  1.00000000  1.00000000  1.00000000   1.00000000         NaN      NaN
C2[1,2]    0.36807151 1.197018e-03 1.336909e-01  0.085009391  0.28188831  0.37566302  0.46306393   0.60617986 12473.90228 1.000010
C2[1,3]   -0.14239190 2.094871e-03 1.911487e-01 -0.497730810 -0.27689751 -0.14864541 -0.01537035   0.24669587  8325.83606 1.000362
C2[2,1]    0.36807151 1.197018e-03 1.336909e-01  0.085009391  0.28188831  0.37566302  0.46306393   0.60617986 12473.90228 1.000010
C2[2,2]    1.00000000 6.909260e-19 8.677163e-17  1.000000000  1.00000000  1.00000000  1.00000000   1.00000000 15772.20630 0.999750
C2[2,3]   -0.51004753 1.404582e-03 1.401675e-01 -0.748983319 -0.61065500 -0.52228036 -0.42029909  -0.20380040  9958.65003 1.000241
C2[3,1]   -0.14239190 2.094871e-03 1.911487e-01 -0.497730810 -0.27689751 -0.14864541 -0.01537035   0.24669587  8325.83606 1.000362
C2[3,2]   -0.51004753 1.404582e-03 1.401675e-01 -0.748983319 -0.61065500 -0.52228036 -0.42029909  -0.20380040  9958.65003 1.000241
C2[3,3]    1.00000000 9.048801e-19 7.444731e-17  1.000000000  1.00000000  1.00000000  1.00000000   1.00000000  6768.86636 0.999750
mu1[1]    95.57517761 1.668159e-02 9.558253e-01 93.715785423 95.01577926 95.56262236 96.12219737  97.48783404  3283.08443 1.002009
mu1[2]     1.88329995 5.539943e-04 3.701418e-02  1.808017780  1.86141766  1.88330704  1.90494729   1.95778940  4464.01413 1.001386
mu1[3]     0.54360908 2.900954e-04 2.056420e-02  0.502194135  0.53178729  0.54368020  0.55568525   0.58442917  5025.06724 1.000151
mu2[1]    95.75226750 1.241873e-01 3.212510e+00 89.807536471 94.36059805 95.72847901 97.08067699 101.22728362   669.16771 1.006042
mu2[2]     1.56043740 7.020145e-03 2.266532e-01  1.132841981  1.46427397  1.55927301  1.65740267   1.99766039  1042.39302 1.001799
mu2[3]     0.50412983 7.496008e-04 3.943456e-02  0.425603647  0.48464568  0.50360980  0.52302921   0.58509991  2767.54094 1.001870

Comparison of observed and predicted % dissolution

The dissolution profiles for each unit are shown in Figure 23.5 below. Each of the 13 panels includes dissolution profiles for 12 tablets with symbols indicating individual observed dissolution values and lines joining corresponding Weibull model based predicted posterior mean values.

FIGURE 23.5: SODF2 data: Dissolution profiles for each batch. Batch numbers are shown above each panel.

Posterior probability of similarity relative to the spheroidal $F_{2}$ similarity region

As described in the text, only the first $N_{t i m e s} = 5$ time points (1, 2, 4, 8, and 16 minutes) are employed for the $F_{2}$ similarity test. The condition of similarity is

$F_{2} = 50 \cdot \log_{10} (\frac{100}{\sqrt{1 + \frac{1}{N_{t i m e s}} \sum_{i = 1}^{N_{t i m e s}} {(w b 2_{i} - w b 1_{i})}^{2}}}) > 50.$

Letting $Δ w b_{i} = w b 2_{i} - w b 1_{i}$ , this condition is equivalent to

$R = \sqrt{\sum_{i = 1}^{N_{t i m e s}} Δ w b_{i}^{2}} < \sqrt{99 \cdot N_{t i m e s}}$

where $R$ is the radius of an $N_{t i m e s}$ dimensional hypersphere centered at the origin (no difference in true mean dissolution profiles).Notice that the definition of similarity now depends critically on the number of dissolution time points included in the calculation with larger values of $N_{t i m e s}$ allowing greater differences at some time points to be compensated by smaller differences at others. In this sense a definition of similarity based on $F_{2}$ seems incoherent. Never-the-less, it continues to serve as a regulatory standard in leu of bio-equivalence.

FIGURE 23.7: IV dissolution similarity of SODF2 made by two manufacturers. Kernal density estimate of the posterior distribution of $F_{2}$

Figure 23.7 provides a kernal density estimate of the posterior distribution of $F_{2}$ . The estimated posterior probability of similarity, the event $F_{2} > 50$ , is thus 0.396, obtained as a simple counting exercise.

There is no universal regulatory standard for the minimum posterior probability of similarity, not should their be. From a scientific and risk management perspective, this minimum should depend on the consequences of non-similarity and decision errors in any particular application. However, this result seems low. Figure 23.6 below shows that the grand mean dissolution difference between sites at 5 minutes is greater than 10% which is consistent with the low posterior probability of similarity.

FIGURE 23.6: Comparison of the mean dissolution profiles of manufacturing sites 1 and 2. Each symbol plots the grand mean percent dissolution over all batches and units for each site and time point.

Posterior probability of similarity relative to a hyper-rectangular similarity region

An alternative definition of similarity is based on the proposal of Sathe et al (1996) who defined a hyper-rectangular region $S R 1$ of similarity as.

$S R 1 = (- 3 {\hat{S}}_{1} \leq Δ μ_{1} \leq 3 {\hat{S}}_{1}) \cap (- 3 {\hat{S}}_{2} \leq Δ μ_{2} \leq 3 {\hat{S}}_{2}) \cap (- 3 {\hat{S}}_{3} \leq Δ μ_{3} \leq 3 {\hat{S}}_{3})$

where ${\hat{S}}_{p}$ are univariate standard deviation estimates that include both between and within lot variation estimated for the pth parameter from reference site 1. We can estimate these limits from the posterior draws as follows:

${\hat{S}}_{p} = \frac{1}{N_{d r a w s}} \sum_{i = 1}^{N_{d r a w s}} \sqrt{{(σ l 1_{p}^{[i]})}^{2} + {(σ u 1_{p}^{[i]})}^{2}}, p = 1, 2, 3$

Here, the superscript $[i]$ indicates the $i^{t h}$ posterior draw of the between-lot or within-lot variances respectively for the $p^{t h}$ parameter from the reference site 1.

FIGURE 23.8: Scatter plot matrix of joint posterior sample of Weibull parameter differences. The red rectangles correspond to bivariate shadows of the $S R 1$ similarity region. The black ellipses are bivariate marginal 90% confidence regions of the plotted MCMC draws.

Again, using the simple counting exercise to determine the proportion of posterior draws lying within the hyper-rectangular similarity region $S R 1$ , we find the posterior probability of similarity based on the $S R 1$ definition to be 0.8870625. This differs considerably from the posterior probability based on the $F_{2}$ definition, but there is no reason to expect agreement because these definitions of similarity are different. In particular, the factor of 3 used in the $S R 1$ definition is quite liberal.

A criticism of the $S R 1$ definition of similarity is that it is based on estimated (plug-in) values of ${\hat{S}}_{p}$ . Thus the definition of similarity will be different for each new set of data. This does not meet the normative requirements for a coherent similarity definition outlined by Eaton et al (2003). A definition of similarity based on the required properties of underlying model parameters would be preferable from a perspective of safety, efficacy, and life-cycle quality control. A simple modification of the above procedure would be to set fixed marginal limits on each of the three Weibull parameters based on the fit-for-purpose requirements of the pharmaceutical product. In principal, the above procedure can accommodate any similarity region whose shape can be expressed algorithmically. This is because the posterior probability of similarity is estimated using a simple counting exercise. It is only necessary to determine whether each MCMC draw is, or is not, contained within the defined region.