Validation of kinetic modeling of progesterone release from polymeric membranes.

Mathematical modeling in drug release systems is fundamental in development and optimization of these systems, since it allows to predict drug release rates and to elucidate the physical transport mechanisms involved. In this paper we validate a novel mathematical model that describes progesterone (Prg) controlled release from poly-3-hydroxybutyric acid (PHB) membranes. A statistical analysis was conducted to compare the fitting of our model with six different models and the Akaike information criterion (AIC) was used to find the equation with best-fit. A simple relation between mass and drug released rate was found, which allows predicting the effect of Prg loads on the release behavior. Our proposed model was the one with minimum AIC value, and therefore it was the one that statistically fitted better the experimental data obtained for all the Prg loads tested. Furthermore, the initial release rate was calculated and therefore, the interface mass transfer coefficient estimated and the equilibrium distribution constant of Prg between the PHB and the release medium was also determined. The results lead us to conclude that our proposed model is the one which best fits the experimental data and can be successfully used to describe Prg drug release in PHB membranes.

Introduction

Mathematical modeling in drug release systems, are of utmost importance in their development and optimization. Its importance lies on predict drug release rates as well as their diffusion from the polymer matrix and to elucidate the physical transport mechanisms involved. The practical benefit of an adequate mathematical model is the possibility to foresee the design parameter effects on drug release profiles.

Furthermore, the in vitro study can provide information about the polymer-drug interaction and could be useful as a preliminary stage to predict in vivo behavior [1]. That is why mathematical models are steadily increasing in importance in academic and industrial fields, with huge future potential.

However, modeling of drug release is complicated because the properties change over time, as for example the shape of the device, which influence on the dissolution and diffusion of the drug. Nevertheless, numerous delivery systems were characterized using partial differential equations to explain their behavior, using analytical or numerical resolution methods [[2], [3], [4]].

Since the pioneer work of Higuchi [5], several empirical and semi-empirical [[6], [7], [8], [9], [10], [11], [12], [13], [14]], as well as mechanistic realistic models [[3], [13], [15], [16], [17], [18]] have been developed. The first ones are explicit equations of drug release amount as function of time, while the latter needs to be solved numerically to obtain the amount and rate of drug released.

However, usually the empirical models cannot describe the entire drug release profile, especially those derived from the power law. Furthermore, the predictive capacity of empirical/semiempirical models is often low due to the mathematical treatment is descriptive unlike mechanistic mathematical theories which consider real phenomena such as diffusion, dissolution, erosion, etc. [[10], [15]].

Generally, drug release occurs in three phases. The initial one can be a burst period, where the dissolved drug can pass easily to the release medium, or a lag time. Then the polymer matrix controls the release mechanism. Finally, the release rate of the drug decreases as it is depleted [19]. A better understanding of controlled-release mechanisms and improved development of technologies will increase the availability of pharmaceutical products [[20], [21]].

In this work we validate a new model, recently developed and published by our research group [22], for progesterone release from PHB membranes, derived from a second order kinetic expression. This kinetic can lump together the main stages involved in release processes. This new model fits experimental data from t = 0 (amount of drug released, M = 0) to t → ∞ (amount of drug released at equilibrium, M). Moreover, the rate of drug released is easily found.

To validate the model, a statistical analysis was conducted to compare the fitting of our model with six different empirical or semi-empirical models. The Akaike information criterion, which considers the number of experimental data and the number of parameters in a particular model, was used to find the equation that yields the best fit. External mass transfer coefficient and equilibrium constant were determined as well.

Materials and methods

Materials

Powder PHB with a molecular weight around 524,000 g/mol was generously provided by BIOCYCLE®, PHB Industrial S.A. (Brazil) with a purity of 99.5% and moisture content below 0.3%. Chloroform, by Cicarelli (Argentina), was used as solvent and Prg as drug (Farmabase, Rovereto, Italy).

Membrane synthesis and characterization

Membranes were prepared by the solution-casting technique. A detailed experimental procedure for the membrane preparation was reported in a previous paper of our research group [22]. Briefly, PHB was dissolved in chloroform at 60 °C for 4 h under reflux and the drug was introduced by direct dispersion. The solution was poured in glass Petri dishes, allowing solvent evaporation at room temperature. Drug content ranged from 23 to 41 wt%. The 41 wt% is the maximum amount of Prg that could be added to the PHB to obtain membranes with homogeneous distribution of Progesterone crystals. Membranes thicknesses (around 110 µm) were taken from cross section SEM images [22].

A complete physical, morphological and chemical characterization of PHB-Prg membranes was performed and PHB-Prg interactions were corroborated. Prg crystals were distributed throughout the membrane thickness, indicating that the drug incorporation in the film was effective [22].

Release experiments

The in vitro release data measurements were performed using pre-weighed pieces of 3x4 cm progesterone loaded membranes of 110 µm thickness on average, placed in contact with 100 ml of release medium in 250 ml beakers (pH 6.8 phosphate buffer solution), at 32 °C and with continuous horizontal stirring in a water bath shaker with controlled temperature. Samples of 3 cm3 were withdrawn and the released amount of progesterone was determined using UV–visible spectroscopy (UV-Visible 2100 C) at 245 nm. The sample volume was then immediately returned to the original solution, to keep constant the total volume, not being replaced by fresh medium. This procedure fulfills batch process conditions (constant mass of drug in the system).

Mathematical analysis

In general, steps involved in drug release processes include drug dissolution, diffusion through the polymeric matrix, eventually polymer swelling or erosion, and transference to the receptor solution at the membrane-fluid interface. A mathematical model can be proposed only when the physical aspects of the involved phenomenon have been properly established. Evidently, hypotheses and assumptions that have to be made influence on the correspondence between the mathematical model and the phenomenon.

The number of models that can be associated with a specific phenomenon depends on researchers' imagination. However, there is no model able to describe all the issues of a given phenomenon, and actually, in most cases it is not necessary. Moreover, when more general the mathematical model, the calculation expressions become more difficult, complicating their practical application.

It should be also noted that not all the aspects of a particular phenomenon have always the same relative importance. In other words, while in some cases the release kinetic behavior is determined by swelling or erosion of the polymer, in other cases the diffusion and dissolution of the drug may play the major role. In fact, different mechanisms can be occurring at the same time or in stages during the release process. It is important to establish these mechanisms for the successful design and manufacture of controlled release systems and to identify potential failure modes.

Usually, as the model approaches more to reality, it becomes more complex. For this reason, simple equations are used to represent the so-called empirical models, which aim to the description of the macroscopic behavior of a phenomenon, without considering the microscopic aspects. They become useful when first studying a phenomenon, or to compare qualitatively different sets of data obtained in the laboratory.

We have proposed a mathematical mechanism that follows a lumped second-order kinetic, model that we will validate in the present contribution, so it is imperative to consider the model equations [22], which are transcribed below. The drug release rate is directly proportional to the square of the amount of drug available in the membrane at each moment (Eq. 1). We found that this model satisfactorily describes processes in which there are various steps involved, and when the drug concentration in the release medium increases steadily [[22], [23], [24], [25]].

where M and M are the total amount of drug released at time t, and the amount of drug feasible to be released at equilibrium, respectively. By elementary integration of this differential equation, between the initial condition (t = 0, M = 0) and any other (t = t, M = M), Eq. 2 is obtained:

Since, M and k are constant in each experimental run, a general model for the amount of drug released as function of time is:where

Therefore, using Eq. 4 or considering t → ∞ in Eq. 3, results:

By simply deriving M versus time in Eq. 3, crucial information as the drug release rate, can be obtained:

Eq. 3 and Eq. 6 can be applied from t = 0 to t → ∞.

Validation of the new model

Typical engineering tools, such as mathematical models, can be very useful to predict the performance of controlled release systems or to measure some important related parameters, such as the diffusion coefficient of the drug [[3], [13], [15], [16]]. For their application, all phenomena governing the release kinetics should be clearly understood. To validate our model, six different models were applied to assess their capacity in fitting the experimental data. The first mathematical model, based on the diffusion front approach, is that of Higuchi [5]. His model, initially conceived for planar systems, was then extended to different geometries and porous systems [26]. The model is based on the hypotheses that a) initial drug concentration in the matrix is much higher than drug solubility, b) drug diffusion takes place only in one dimension (edge effects must be negligible), c) solid drug particles are much smaller than system thickness, d) matrix swelling and dissolution are negligible, e) drug diffusivity is constant, and f) perfect sink conditions are always attained in the release environment [[10], [26]].

If a drug is evenly dispersed in a non-degradable polymeric matrix, such as in the case of some membranes, it has to dissolve and then diffuse throughout the polymer to be released in the medium, since the polymer will not erode, or will do it in a long time.

The Higuchi equation, which describes the release transport when it is a diffusion-controlled process, establishes a direct relationship between the release rate and the square root of time [27]:where M is the amount of the released drug until time t, A is the release area, D is the drug diffusion coefficient in the polymer matrix, C0 is the initial drug concentration in the matrix whereas C is the drug solubility in the polymer matrix. Interestingly, this model shows that M depends on square root of time and coincides with Fick's solution when less than 60% of the drug is released [11].

Ritger & Peppas described an empirical and simple equation for the first 60% of the release curve [7], and according to the diffusional exponent value it can be elucidated which release mechanism took place. This model is the so called power law presented by Peppas and coworkers [[6], [7]]:where a is a constant and n is the diffusional exponent related to the drug release mechanism. It should be noticed that this equation is usually presented as the ratio of , where is the amount of drug released at infinite time. However, to compare with our model, the value of was directly included in the equation's parameters in Eq. 8 and in the following models.

According the geometry, n value in Eq. 8 is equal to 0.5 for a thin film, 0.45 for a cylinder, and 0.43 for a sphere, when Fickian diffusion takes place [26]. If the n value is higher, non-Fickian release takes place.

When the release exponent (n) is equal to 0.5, an equivalent Higuchi equation (Eq. 7) is obtained (Eq. 9):

Due to its simplicity, this equation is widely used in the pharmaceutical area, however many assumptions were made to arrive to this expression. For this reason, an incorrect diffusion mechanism can be assumed if they are not considered [28].

Values of 0.5 < n ≤ 1.0 are indicative of superposition of diffusion and swelling controlled drug release, and an anomalous transport is observed. When n = 1.0, it corresponds to a zero-order release mechanism.

A model that account for the coupled effect of Fickian diffusion and polymer relaxation contribution, is that based on Peppas-Sahlin [29], equation (Eq. 10):where a and b are the kinetic constants related to the Fickian and non-Fickian diffusional contribution, respectively. Regardless of the device geometry used, the exponent for the polymer relaxation transport mechanism is twice the diffusion Fick mechanism. As can be seen, the two phenomena controlling the release can be considered as additives and n is the pure Fick diffusion exponent.

A simpler expression of Eq. 10, is that where the exponent n is set in 0.5 [[8], [29]]:

When the effect of the external (interface) mass transfer resistance is significant, the model that takes into account the coupled effects of drug diffusion through the polymeric membrane and the interface transport, is represented by Eq. 12 [[9], [11], [12], [13], [14]]:

Another simple and useful semi-empirical model [9], for the amount of drug released in a slab devise is:

Note that a, b, c and n, are positive parameters, otherwise would go against mass conservation principle [[11], [30]].

Statistical analysis

In mathematical modeling, the main problem is to discern which model is the best one among several yielding good fits. One possibility would be to use the method of the sum of the squared residuals (SSR) to find the model that best explains the experimental data (minimal value for the SSR). The release tests results were analyzed for the seven models by a nonlinear regression analysis, through the Polymath 6.0 program, and the SSR values were calculated.

However, if the model has too many parameters, a small SSR value could be obtained, so this parameter by itself does not weigh the complexity of the model. That is why it is necessary to use a discriminatory criterion.

Different approaches can be followed for this purpose, but due to its simplicity, the Akaike's method [[30], [31]] is the more convenient. Assuming that the random errors follow a Gaussian distribution, the Akaike number AIC (Akaike Information Criterion) is defined in Eq. 14 and it was used in the discrimination analysis [[30], [32], [33], [34]].where N is the number of experimental data, SSR is the square residuals sum and the number of parameters in the model is represented by p. The model that best represents statistically the drug release mechanism is the one with minimum AIC. The AIC criterion considers both, the number of parameters and experimental data in the model. Therefore, the AIC criterion is better than simple comparison of the SSR values of the different models.

Results and discussion

Progesterone release

In a previous article of our research group, the lumped second order mathematical model was presented and used to fit the release data for 23, 29 and 33 wt% Prg in PHB membranes [22]. In this contribution, the maximum Prg loading (41 wt%), where the Prg crystals are still uniformly distributed across the membrane thickness, was included.

Fig. 1 shows the progesterone mass release from PHB-Prg membranes, with different Prg content, as function of time. It can be seen that the amount of hormone released, is influenced by the original Prg load in the membrane; i.e. higher Prg load, lower drug release. This behavior can be attributed to an excess of progesterone crystals in the membrane, and particularly at the interface. Dissolution of these crystals in the release medium is hindered by the very low Prg solubility in the polymer and receptor solution, although drug solubility in PHB is higher than in the release medium, as it will be discussed later. Thus, the surface available for drug release decreases as Prg loading increases. When the dissolved and non-dissolved Prg coexist within the polymeric matrix, the dissolved drug is the only available for diffusion [22].

Fig. 1

Fit of the proposed model to experimental Prg release data. Symbols are the mean value experimental data and their sizes represent the standard deviation. Lines represent the theoretical release predictions with nonlinear regression fit developed in this work (Eq. 3).

Fitting the experimental data with our model through the non-linear regression analysis (Polymath 6.0 program), the lines in Fig. 1 were obtained, showing a good fit of the model (Eq. 3) to experimental data. Thus, parameters a and b are determined from which M values were calculated (Eq. 5).

The same non-linear regression analysis was carried out with the six models considered. The values of the parameters for each model as well as the R2, SSR and AIC values are reported in Table 1, Table 2, Table 3, Table 4, for the four Prg contents in the PHB membranes.

Table 1

Estimated parameters, SSR and AIC values for the models, when membranes were loaded with 23 wt%Prg.

Model	M∞	a	b	n	c	R2	SSR × 102	AIC
Eq. 31.9020.013570.00714	1.902	0.01357	0.00714	–	–	0.998	0.214	−51.32
Eq. 8	–	0.05032	–	0.58600	–	0.996	0.477	−44.13
Eq. 9	–	0.07678	–	–	–	0.986	1.923	−33.56
Eq. 10	–	0.02843	−1.489E−4	0.75928	–	0.999	0.025	−60.38
Eq. 11	–	0.06231	0.00119	–	–	0.995	0.325	−40.47
Eq. 121.3751.375380.97633	1.375	1.37538	0.97633	–	0.00836	0.998	0.325	−45.56
Eq. 131.3331.33339	1.333	1.33339	–	–	0.00920	0.996	0.488	−43.89

Table 2

Estimated parameters, SSR and AIC values for the models, when membranes were loaded with 29 wt%Prg.

Model	M∞	a	b	n	c	R2	SSR × 102	AIC
Eq. 31.8130.012840.00708	1.813	0.01284	0.00708	–	–	0.999	0.168	−53.49
Eq. 8	–	0.04638	–	0.5918	–	0.988	1.430	−34.23
Eq. 9	–	0.07283	–	–	–	0.976	2.947	−29.72
Eq. 10	–	0.01775	−6.983E-5	0.8668	–	0.999	0.013	−65.80
Eq. 11	–	0.05883	0.00116	–	–	0.985	1.816	−32.08
Eq. 121.2681.268460.99375	1.268	1.26846	0.99375	–	0.00911	0.999	0.145	−52.82
Eq. 131.2601.25993	1.260	1.25993	–	–	0.00931	0.999	0.182	−52.77

Table 3

Estimated parameters, SSR and AIC values for the models, when membranes were loaded with 33 wt%Prg.

Model	M∞	a	b	n	c	R2	SSR × 102	AIC
Eq. 31.7670.011630.00658	1.767	0.01163	0.00658	–	–	0.996	0.719	−50.29
Eq. 8	–	0.05950	–	0.52138	–	0.991	1.693	−40.86
Eq. 9	–	0.06676	–	–	–	0.990	1.854	−41.86
Eq. 10	–	0.03028	−1.703E−4	0.71091	–	0.997	0.542	−51.39
Eq. 11	–	0.06496	0.00012	–	–	0.990	1.815	−40.10
Eq. 121.3511.350830.96746	1.351	1.35083	0.96746	–	0.00686	0.994	1.184	−42.80
Eq. 131.3191.31890	1.319	1.31890	–	–	0.00761	0.992	1.540	−41.90

Table 4

Estimated parameters, SSR and AIC values for the models, when membranes were loaded with 41 wt%Prg.

Model	M∞	a	b	n	c	R2	SSR × 102	AIC
Eq. 31.6310.010700.00656	1.631	0.01070	0.00656	–	–	0.998	0.360	−57.88
Eq. 8	–	0.05432	–	0.52302	–	0.9883	1.887	−35.70
Eq. 9	–	0.06148	–	–	–	0.987	2.047	−36.89
Eq. 10	–	0.02403	−1.215E−4	0.74680	–	0.998	0.385	−55.15
Eq. 11	–	0.05991	0.00010	–	–	0.987	2.017	−38.94
Eq. 121.2331.232730.97591	1.233	1.23273	0.97591	–	0.00712	0.996	0.628	−49.77
Eq. 131.2131.21281	1.213	1.21281	–	–	0.00766	0.995	0.791	−49.23

Table 1, Table 2, Table 3, Table 4 show that the model given by Eq. 10 is the one that gives the smallest AIC values practically for all Prg loads. However, parameter b is always negative in Eq. 10. Obviously, the minimum AIC value does not guarantee the model reliability, since the parameters also must assume reasonable values (in the physics sense). Therefore Eq. 10 is not physically consistent, as previously indicated, and should be discarded. Comparing with the other five models, our proposed model (Eq. 3) is the one with minimum AIC value, and therefore it is the one that statistically fits better the experimental data obtained for all the Prg loads tested.

Even though mathematical models can always be getting better, it is important to remember that a widely applicable complicated model may be accurate but not useful. In this context, the simplest model fitting properly the experimental data can be the best one when R&D purposes are pursued. Equation 3 is both simple and accurate; it only has two parameters and could be applied in the whole time range.

Normalized release rate

In a previous paper [22], we have presented a valuable parameter, the normalized release rate per unit of effective surface area (R, mg/cm.min), given by Eq. 15where A is the membrane effective surface area, calculated through Eq. 16,

A is the total membrane surface area and (v/v) is the volume fraction of PHB in the membranes, considering that due to the low Prg solubility in the PHB, the density of both PHB and Prg could be considered as pure compounds (1.300 g/cm3 and 1.166 g/cm3 respectively). The (v/v) values are 0.750, 0.687, 0.645 and 0.563 for 23, 29, 33 and 41%Prg repectively.

The R value must be the same, regardless of Prg load. Fig. 2 shows that the R values calculated as function of time, for the four Prg loads used, gives a unique line that fits all the experimental data points. These explain exactly what we stood previously.

Fig. 2

RRef vs. time. Symbols are the mean value experimental data and their sizes represent the standard deviation.

External mass transfer coefficient and equilibrium constant

Prg release rate should be equal to the rate of drug transfer through the external membrane effective surface area and the solution. This relation is given by Eq. 17.where C and C are the Prg concentration at the interface and in the bulk fluid, respectively, and kc is the external mass transfer coefficient.

At the beginning of the process (t = 0), C = C and C = 0. Here, C is the solubility of Prg in the release medium. Then, Eq.17 gives:

Therefore, according to Eq.6:

Considering that C = 0.021 mg/cm3 [22], the mean value of kc is (6.1 ± 0.1) × 10 cm/sec. This value agrees with kc values found in systems with low rate of mixing, characteristic of the horizontal shaker [[14], [35], [36]].

Another interesting information that can be obtained is the equilibrium distribution constant (K) (partition coefficient) of Progesterone between the PHB and the release medium. This thermodynamic quantity is estimated as the drug solubility's ratio in the membrane material and the release medium. The approximate value of K is:where is the Prg solubility in PHB [22], (equivalent to C in the Higuchi equation, Eq. 7). Values of K for Prg distribution between lipids/buffer solutions ranging from 1200 to 2000 were presented for Heap et al. [37]. Also a partition coefficient as high as 6918 (log K = 3.84) for Prg between Dulcoo's PHS (pH 7.40) and octanol, was reported [38].

Conclusions

The new model, developed for drug release processes and derived from a second order kinetic expression, was validated. This model was compared with several mathematical models, for progesterone release from PHB membranes. The model is simple; it only has two parameters and can describe the entire drug release profile, even for t = 0, unlike the power law expression. Via the Akaike information criterion (AIC), it was demonstrated that our model is the one that best fits experimental results. External mass transfer coefficient and equilibrium distribution constant were determined as well.

Declaration of interest

The authors declare that there is no conflicts of interest. The authors alone are responsible for the content and writing of this article.