Hierarchical Bayesian approach applied to the formulation of sustained-release suppositories and dissolution profile modeling

The hierarchical Bayesian modeling approach was used to select the appropriate empirical kinetics model of sustained release and to optimize the in vitro dissolution rate of the sustained-release suppository by controlling the composition of Eudragit L-100 and Eudragit S-100 in the experimental mixture. Thirteen formulations of suppositories were prepared with 2 g (10%) mixture of Eudragit® R-100 and S-100 according to a personalized mixture experimental design. The cumulative release of active ingredient was measured at five times (20, 50, 80, 160, and 235 minutes). The best model was selected using Rsq (Adjust) and akaike information criterion for standard method and by using the weight of widely applicable information criterion (WAIC) and leave-one-out (LOO) cross-validation for the Bayesian approach. Frequentist approach gave three best model depending on the formulation. Compared to this, the Bayesian method was able to define a single model, which is the first-order model. The relative probability of this model is 0.97, 0.99 based on the WAIC, and LOO, respectively. The relationship between K1 (Release rate constant) and the quantities of the two Eudragits is quadratic, for Eudragit_L, Qrelease (%) = 0.0031X 2–0.0026 X + 0.0069 and X is the Eudragit L100 and K1 (Rate release) = 0.41 minutes −1. The Bayesian method allowed finding the most adequate model among several models that can be generated by the standard frequentist approach.


INTRODUCTION
Mathematical models play a crucial role in the interpretation of drug release mechanism (Gouda et al., 2017;Grassi et al., 2006;Paarakh et al., 2018;Varma et al., 2004). The drug release kinetics for given formulations has been usually fitted to predefined semi-empirical models using non-hierarchical frequentist modelling approach. The appropriate model for each formulation is selected based on the highest Rsq (Adjust) and the lowest Akaike Information Criterion (AIC) (Gouda et al., 2017;Paarakh et al., 2018;Victor and Francis, 2017). Such a situation may lead to variation in the best-fitting model among formulations, in particular when the number of points in the curve is limited. This is not meaningful biologically given the interdependency between data and the model convergence may be problematic (Piray et al., 2019). Hierarchical Bayesian modelling can solve these issues by expressing the model over all formulations and consequently addressing both model fitting and model comparison within the same framework (Piray et al., 2019). Hierarchical Bayesian modelling approach assumes that the formulations share the same best model and differ only in some or all parameters values of the model (Andrews et al., 2002;Millar, 2009;Wiecki et al., 2013). The differences between the two statistical concepts have their root in differing conceptions of probability. In practice, frequentist approach quantify the properties of data-derived quantifies, as for Bayesian approach quantify the incertitude of measure (VanderPlas, 2014).
For pharmaceutical development, the application of Bayesian statistics remains modest and mainly affects the field of clinical trials and not the field of formulation. This will be one of the strongest points of our formulation study (Natanegara et al., 2014;Price and LaVange, 2014) .
The use of polymers in the drug delivery system increased during the last years as a result of the advantages offered by these molecules (Andrews et al., 2002;Malekmohammadi et al., 2019;;Zuo et al., 2014). For sustained-release system, Eudragits, good safety profile molecules, may give the formulator the possibility of choosing between different categories depending on the purpose of the formulation such as oral, colon, rectal, and even transdermal (Joshi, 2013;Nikam et al., 2011;Thakral et al., 2013). For clonic targeting, both Eudragit L-100 and S100 are widely used for molecule delivery as they offer muco-adhesiveness and pH-dependent release (Khan et al., 1999;Lee et al., 2020). Only one study has been reported in the scientific literature on the use of these polymers in rectal drug delivery (Baloglu et al., 2002). In this study, authors prepared and evaluated the kinetics of sustained-release suppositories of polyethylene glycol (PEG) 400, PEG 4000, Witepsol, and Eudragit L-100. In our study, the choice of these Eudragits comes from the fact that the physiological conditions of the colon (especially the distal colon) and those of the rectum were similar for the pH (pH = 7-8); hence, the interest of these polymers called pH-Dependent Drug Delivery Systems in the formulation of Hydrophilic suppositories (Hua, 2019;Lee et al., 2020). By optimizing the mixture composition of these polymers (EL-100 and ES-100) and the total polymers content of the matrix, a specific release rate for the active molecule can be obtained by modifying the diffusivity through the gelled of the matrix structure (Di Colo et al., 2002;Misra, 2014;Nikam et al., 2011;Surti et al., 2020).
The novelty of this study is the combination of two Eudragit polymers in the formulation of sustained release suppositories and the application of the Bayesian hierarchical approach to the modeling of the release kinetics. The first objective of this study is to use the hierarchical Bayesian modeling approach to select the appropriate empirical kinetics model of Barbital sustained release. The second objective is to optimize the in vitro release rate of sustained-release suppository by controlling the composition of E-L100 and E-S100 in the experimental mixture.

Design of experiments and statistical analysis
In the present study, we used a personalized design of experiment with constraint and with process factor, which was time. The quantities of the two Eudragits varied from 0 to 2 g each, with upper constraint of 2 g for the sum of the two Eudragits L and S, based on preliminary tests. Beyond this limit, there would be little flow, poor filling of the metal molds, and an unacceptable appearance of the prepared suppositories. The quantity of Macrogol 6000 and Barbital were set to 17.4 and 0.6 g, respectively (Table 1). For each of the 13 formulations, three experiments was made and five sampling performed at 20, 50, 80, 160, and 235 minutes.

Preparation of suppositories by fusion molding
The preparation consisted in melting the quantity of PEG 6000 of each formulation at 46°C in the stainless steel capsule, then adding the other components under agitation (risk of gelation of the mass) and without incorporation of air bubbles. The prepared suppositories was scraped and stored in glass bottles.

Chemical and pharmaco-technical assay of the prepared suppositories
The prepared suppositories have been evaluated following the USP monographs for appearance, weight variation, content uniformity of dosage form, hardness, disintegration time, and dissolution assay. The Barbital content was determined using a UV-visible spectrophotometer. The hardness was performed with the Erweka AR 400 hardness tester (Erweka Apparatebau-GmbH Germany). The disintegration was performed in a 6.8 pH buffer solution at 37°C (+/−0.5) using the U.S.P. tablet disintegration apparatus (SOTAX DT 3, Heusenstamm, Germany).

Dissolution testing
The suppository was placed in the basket, lowered to the precise level in 1L vials filled with 500 ml phosphate buffer (pH 7.2), and conditioned at 37°C ± 0.5. The basket apparatus was set at 50 rpm. Three millilitre of the dissolution medium was withdrawn with a syringe at defined time intervals (t20, t50, t80, t160, and t235 minutes). The volume withdrawn was replaced by its equivalent in phosphate buffer. The optical density of the filtrate (filters 0.24 μm) was measured at 245 nm with the UVvisible spectrophotometer (Spectrophotometer 6305 UV-VIS 198 at 1000 nm Jenway, Bibby Scientific France SAS) (Onyeji et al., 1999;Özgüney et al., 2007;Ranjita and Kamalinder, 2010).

Modeling dissolution kinetics
For modified release suppositories, the release of the active principle is due to different parameters such as disintegration, diffusion of the active principle through a membrane, swelling, and diffusion through a polymer, erosion of a matrix (Costa and Faisant et al., 2002;Varma et al., 2004). Various kinetics models were used to describe the release kinetics (Faisant et al., 2002;Peppas and Narasimhan, 2014).
The zero-order model (Equation 1) is applied to pharmaceutical forms that do not disintegrate and that release the product over a long period of time. These are of matrix modifiedrelease forms with poorly soluble active ingredients, coated forms and osmotic forms (Moodley et al., 2012). The first-order model (Equation 2) was proposed by Wagner for forms whose release depends on pH (Basak et al., 2008). The Higuchi model (Equation  3) is applicable if diffusion is the only mechanism responsible for the release of the active ingredient (Craig, 2002). The Korsmeyer-Peppas model (Equation 4) is used to classify drug release from polymers as fickian and non-fickian (Güneri et al., 2004). The Hixson-Crowell model (Equation 5) describes the release of systems in which the surface area and diameter of the matrix change over time (Dash et al., 2010).
The so-called Zero Order: The first-order model: The Higuchi kinetics model: The Korsemyer-Peppas kinetics model: Hixson-Crowell model (5) where Q is the cumulative percent of drug release of Barbital, Q 0, Q H , Q KP , and Q HC are the intercepts of the models, k 0 , K 1 , K HC , and K KP are the rate of release, n is the diffusional exponent indicating the drug-release mechanism. When the values of n equal to 0.5, the drug release, follow the Fickian diffusion law. When 0.5 < n < 1, the drug release follow a non-Fickian diffusion. When n = 1.0, the drug release is caused by swelling of the polymer (Costa and Lobo, 2001;Zhang et al., 2010).

Statistical analysis
Frequentist approaches study the behaviour of the data, i.e., the probability of repetition of events. Bayesian approaches are based on the measurement of uncertainty.

Frequentist approach
Data of each Formulation were fitted individually to the five above-mentioned models using DD-Solver which is an add-in program in Microsoft Excel. DD-Solver uses the nonlinear leastsquares fitting technique, and models parameters were determined based by minimizing the sum of squares errors (Nakata, 2010;Zuo et al., 2014). The initial values of the model were set individually before the iterative optimization. The best model for each formulation is selected based on several criteria such as R-sq Adjusted and AIC. Both Rsq and AIC present a limitation in the model selection and inference process when the number of parameters are different between candidate models (Nakata, 2010).
In addition, the DD-solver does not allow for hierarchization of one or all model parameters. Given that, the selected model (bestfitting model) may vary from formulation to another in particular when the number of points in the curve is limited (Nakata, 2010).

Bayesian approach
Data of all formulations were fitted to the five above mentioned models using a hierarchical Bayesian model using the PyMC3 packages of Python (Patil et al., 2010;Salvatier et al., 2016). This approach assumes that the formulations share the same best model and differ only on values of some or all parameters of the model, for each of the five models. Prior parameters of each model were as follows: the variation was assumed zero and model parameters were set to the initial value obtained when the model was fitted to all data.
The trace plot, the posterior distribution of the parameters over each Markov chain Monte Carlo (MCMC) and the convergence of the chains were used to examine the performance of the models (Bates and Campbell, 2001;Cowles and Carlin, 1996). MCMC is a class of algorithms which can efficiently characterize even high dimensional posterior distributions through drawing of randomized samples such that the points are distributed according to the posterior. For each model, chain MCMC of 2000 cycles were considered with a burn-in period of first 200 cycles.
We used two methods for Bayesian model selection such as the widely applicable information criterion (WAIC) and Leaveone-out cross-validation (LOO) Watanabe and Opper, 2010). The methods evaluate models based on their predictive validity. LOO method examines each left out piece of data, the performance of the model. WAIC uses the entire posterior distribution of models to assess model performance (Gelman et al., 2014;Luo and Al-Harbi, 2017;Vehtari et al., 2015). Widely available information criterion is calculated as Gelman et al.
(2014) as: The lppd is the log point-wise predictive density and can be approximated as: lppd = ∑log (E post p (yi|θ)) and p WAICn = ∑var post log(p(yi|θ)) The simpler way to perform model selection is to rank models based on their WAIC or LOO value and to choose the model with the lowest one, However, this simple method do not take in to account the uncertainty of WAIC of each model given the standard error of WAIC. As such, we performed a model selection calculating an average weight for each model relative to the best model (the lowest WAIC or LOO): where wi is the weight of the model i, which corresponds to its relative probability relative to all tested models. For the model selection, we have five hierarchical models cited above, in addition we add the corresponding pooled models for each category of the model, which aim to test the hypothesis of the significant effect of the quantities of the two Eudragits on the rate of release. In other words, the selection of hierarchical model corresponds to a significant effect of Eudragits quantities on release rate; inversely the selection of pooled model corresponds to the lack of effect of Eudragits quantities on the kinetic of release.
The effect of Eudragits on the release rate was assessed by investigating the relationship between the parameters of the selected model and the weight of EL-100 or ES-100 in the formulation using the non-linear regression technique in Scipy package.

RESULTS AND DISCUSSION
All suppositories satisfy the requirement of the European pharmacopeia content uniformity, weight variation, hardness, and disintegration time. Although, the uniformity of content varied from 96.1% to 101.5%, and the weight variation from 2.23 to 2.67 g. The hardness of all formulations was superior to 4.8 kg and the disintegration time ranged between 30 and 203 minutes.

In vitro release modeling
In vitro dissolution tests have shown that the prepared suppositories have an sustained-release effect, even if it does not contain Eudragit (formulation No. 13 with PEG 6000 only). The percentage of Barbital dissolved in 4 hours (235 minutes) was more than 90% for four formulations (F3, F5, F12, and F13), more than 80% for formulations F2, F6, F7, and F8 and less than 80% for five formulations (F1, F4, F9, F10, and F11).

Frequentist model selection approach
In vitro release characteristics of sustained-release suppositories were evaluated to determine the best-fitted kinetic model ( Table 2). The release rates of all formulations were best described with first-order model, which presented an Rsq (Adjust) > 0.90 and AIC < 16 model (Fig. 1). Based on the Rsq (Adjust) ,the first-order model was the best model for formulations 1, 2, 3, 4, 5, 8, and 13, the Higuchi model is best for formulation 9, the Korsmeyer-Peppas model was best for formulations 6, 10, and 11 and the Hixson-Krowell model was best for both formulations 7 and 12. Based on AIC, the first-order model was best for formulations 1, 2, 3, 4, 5, 8, and 13, Korsmeyer-Peppas is the best model for formulations 9, 10, and 11 and finally the best model for formulations 6 and 7 was Hixson-Krowell (Zhang et al., 2010;Zuo et al., 2014). This showed that the release for the same formulations can have different best models in vitro, which is one of the weaknesses of the Frequentist approach. In fact, the findings from this approach made no clear biopharmaceutical meaning.

Bayesian model selection approach
The first order model presented the lowest value of WAIC and LOO. The standard errors (SE) varied from 3 to 15. The weight of the model was 0.97 and 0.99 based on both WAIC and LOO, respectively. The relative probability (weight) of the selected model was greater than 0.5, indicating that no competing model was possible. The LOO and WAIC estimate the same predictive performance criteria. The selection of the pooled model represented the H0 hypothesis, i.e., K 1 is not affected by the quantities of the two Eudragits, while the selection of the hierarchical model represented the alternative hypothesis, i.e., K 1 is affected by the quantities of the two Eudragits. According to results, the alternative hypothesis H1 was accepted, the quantities of the two Eudragits have a significant effect on K 1 . (Tables 3  and 4). The Bayesian estimation of the parameters of the selected model was completed based on the posterior normal distribution of the estimate. The Markov chain has stabilized as indicated by the trace plot (Fig. 2). The posterior distribution of the parameter of the best model (K 1 ) was normal for all formulations. The estimated K1 varied from 0.006 to 0.012 and it SE from 0.0001 to 0.001. The SE of the model was 6.154 based on WAIC and 6.4 based on LOO, which confirms the strength of the results. LOO Table 2. The Rsq (Adjust) and AIC for the release kinetics of Barbital using different mathematical models obtained by DD-Solver. and WAIC estimated the same predictive performance criteria and were asymptotically equal (Gelman et al., 2014;Vehtari et al., 2016;Vehtari and Ojanen, 2012). The hierarchical method gives a single stable model (First-order model) for the same series of formulations whereas the standard method (frequentist approach) gives several models. The K 1 of the first-order model estimated by frequentist and Bayesian method were strongly significantly related (Fig. 3). The Bayesian method gave a good estimate of observed data (Fig.  4) and was used for further analysis (VanderPlas, 2014). The relationship between K 1 and the quantities of the two Eudragits was quadratic (Equation 10). The release rate (K 1 ) was low for Eudragit L values up to 0.5 g (2.5%) and increases until it reaches its maximum at 2 g (10%) of this Eudragit, on the contrary for Eudragit S, K 1 decrease from 0 g to reach it minimum at 1.5 g (7.5%) (Fig. 5).
with K = −2b/a = 02*0.0026/0.0031 = 0.41 By this equation which relates the rate of dissolution to the quantity of EL-100, and therefore the mixing of the two Eudragits (total =10%). Drug release from suppository bases generally depends upon the drug solubility in the base and chemical composition of the base. Barbital is a drug in the sodium salt form with a high affinity for hydrophilic bases, i.e., PEG 6000. PEG bases are also known to have a solubility effect, which may partly explain the higher drug release rates of t20 minutes formulations. Both Eudragits create a swellable matrix through which water-soluble active ingredients are released. The mixture advantage became in the difference in solubility of the two polymers in a range of pH =7 to 8, i.e., with EL-100 Intestinal soluble-fluid from (pH 6) and ES-100 is intestinal soluble-fluid    Colo et al., 2002;Hong et al., 2013). The firstorder kinetics model as validated by the Bayesian method as the most adequate model to describe and predict the release of active ingredient, this result is due to the dissolved matrix as a function of the pH of the rectal region (Dash et al., 2010). Gibaldi and Feldman were the first to propose first-order kinetic models of liberation in 1967, followed by Wagner in 1969, for this model the rate of liberation depends on concentration (Costa and Lobo, 2001).

CONCLUSION
In this study, Bayesian approach was applied to drug release modeling and formulation of sustained-release suppositories containing a mixture of EL-100 and ES-100. This approach allowed finding a unique and the most adequate model, which was the first-order model, when the standard frequentist approach gave three best models. At the same time, the Bayesian approach is a tool for predicting the optimal composition for the formulation of a suppository with prolonged effect. In drug development field, this represented an original and new way to master the kinetics modeling of sustained-release products in order to minimize the number of experiments, save time, and save expensive chemicals.

ACKNOWLEDGMENTS
We thank Dr. Lahcen Benomar (University of Quebec in Abitibi-Temiscamingue) for his help with statistical analyses. Special thanks are also extended to Dr. Mohamed Taha Moutaoufik (University of Regina) and Read Elfarjani ( TELUQ University) for their valuable comments and suggestions. We thank the reviewers and the editor for relevant comments that helped to improve the content of this article.