Dissolution Kinetics of a BCS Class II Active Pharmaceutical Ingredient: Diffusion-Based Model Validation and Prediction.

In this work, a diffusion-theory-based model has been devised to simulate dissolution kinetics of a poorly water-soluble drug, ibuprofen. The model was developed from the Noyes-Whitney equation in which the dissolution rate term is a function of the remaining particulate surface area and the concentration gradient across the boundary layer. Other dissolution parameters include initial particle size, diffusion coefficient, material density, and diffusion boundary layer thickness. It is useful for predicting nonsink circumstances under which pure API polydisperse powders are suspended in a well-mixing tank. The model was used to compare the accuracy of simulations using spherical (single dimensional characteristic length) and cylindrical particle (multidimensional characteristic lengths) geometries, with and without size-dependent diffusion layer thickness. Experimental data was fitted to the model to obtain the diffusion layer thickness as well as used for model validation and prediction. The CSDs of postdissolution were also predicted with this model, demonstrating good agreement between theory and experiment.

Introduction

In the pharmaceutical industry, dissolution processes are of interest for many different reasons. For instance, for solid/semisolid and suspension dosage forms, it is only dissolved API which can be ultimately absorbed into the body, and it is fundamental and essential to investigate the dissolution behavior of API in solutions to understand the overall drug absorption in biological systems.[1,2] Particle size and morphology of a drug are important factors that affect thermodynamic solubility and oral absorption of drug molecules. Correlations between particle properties of an API and the dissolution profile of its dosage form are quite complex in general, as the relationship may depend not only on the dissolution of the API itself but also on the disintegration of the dosage form and interactions between the API and the excipients.[3]

Dissolution also plays an important role in industrial crystallization from solution such as in solvent-mediated polymorphic transformations[4−6] as well as in fines removal.[7−9] Characteristics of the final crystallized particles and downstream processability highly depend on the achieved particle size and shape distribution as a part of such processes (PSSD).[10] Cyclic processes, in which dissolution stages are included, have been reported to modify the PSSD systematically, as dissolution can add a degree of freedom to manipulate the crystal size and shape.[11−13]

Kinetic models of the dissolution process are valuable tools to assist with the characterization of experimental work, as such models provide useful information to design and optimize processes as well as allow the development and testing of feedback control strategies.[7,10] The dissolution rate is affected by many factors, such as the CSD and morphology, the agitation of the suspension, and undersaturation levels.[14] In addition, the solvent, additional additives, temperature, and pH values will also have significant impact on the dissolution kinetics. Suitable mathematical model frameworks can be summarized into a number of distinct categories as follows:

Diffusion theory on the basis of the Noyes–Whitney equation,[1,3,15] in which the dissolution rate is proportional to the concentration gradient and the surface area of the solid.

Population balance equation (PBE), which describes the evolution of the particle population.[6,9,10] In this framework, the dissolution rate is defined as the rate of change of particle size by single or multiple characteristic particle dimensions.

Dynamic modeling which calculates the API shape evolution on different molecular faces.[16] The dissolution rate is determined by the change in the perpendicular distance from a central point within the crystal to different faces.

Hydrodynamic modeling using computational fluid dynamics (CFD),[17] in which the velocity of the fluid and the thickness of the diffusion layer are examined in detail. Such hydrodynamic conditions and velocity profiles can be well described within a standard vessel and paddle apparatus.

Other frameworks include quasi-steady state models (QSM),[2] surface detachment controlling models,[18] and surface complexation mechanism models,[19] among others.

Although various attempts have been made in the literature to model the dissolution kinetics of APIs to fully predict the required release profile of the pharmaceutical product, so far no one has succeeded in developing a universal approach which directly correlates to the API powder specification (i.e.: its CSD, morphology, etc.) and dissolution release data. In this paper, the weight fraction of different particle sizes in a polydisperse powder together with Noyes–Whitney parameters (diffusion coefficient, solubility, density of the drug, boundary layer thickness, and dissolution volume) has been used to predict dissolution kinetics of an API. Such a simulation can be a valuable tool with respect to gaining a better understanding of the dependence of dissolution on particle size and geometry.

The purpose of this contribution is to address the modeling of dissolution processes of both polyhedral and needle-like crystals on the basis of diffusion fundamentals. First, the theoretical simulations are calculated assuming a spherical geometry and that the particles adhere to a normal distribution. Next, diffusion-based kinetic models are developed from experimental data using a well-established parameter estimation technique. This approach is applied to the dissolution of the polyhedral compound IBU in a phosphate buffer. To the best of the author’s knowledge, no multidimensional dissolution rate study has been previously presented for ibuprofen in the literature. The models developed are validated by a range of experiments covering different hydrodynamic mixing conditions, CSDs, and morphologies. In addition, the prediction of the dissolution under differing conditions of initial crystal seed mass loading is simulated using the model. The novelty of this work lies mainly in the fact that a multi-dimensional dissolution model is developed on the basis of diffusion theory to predict dissolution kinetics under nonsink circumstances, with the consideration of size-dependency and geometry.

Theory

According to the Noyes–Whitney equation, the rate of dissolution of a solid is dependent upon its solubility, the concentration of solute in solution at a particular time, diffusivity, and the surface area of the solid. For a drug of low aqueous solubility, the particle size and the resulting surface area could have a significant effect on the rate of dissolution and therefore affect the bioavailability of the drug.[20] As dissolution is the reverse process of crystal growth, it involves two steps:[21] (1) surface reaction and disintegration of the surface species and (2) mass transfer of these species into the bulk solution across the diffusion layer, as presented in Figure .

Figure 1

Dissolution kinetics model.

It is typically assumed that the overall process is limited by the slower rate of these steps. Assuming that the surface reaction rate rs and the mass transfer rate rm are proportional to S, the surface area of the undissolved crystals per unit volume of solution, and using the difference of concentration as a driving force, two equations were obtained as follows[22]where C* is the solubility of the solute at a certain temperature, Cs is the concentration of the solute at the surface of crystals, and Cb is the concentration of the bulk solution. ks and km are surface reaction and mass transfer rate constants.

If ks ≫ km, the dissolution process is dominated by the mass transfer mechanism, whereas if ks ≪ km, the surface reaction mechanism controls the process

It is also assumed that under quasi-steady (pseudo) state conditions, the dissolution rate r is equal to the surface reaction rate rs and the mass transfer rate rm, namelywhere m is the mass of the undissolved crystals and V and t are the volume of the solution and the dissolution time, respectively.

A combination of eqs –3 giveswhere k is the dissolution rate constant. The mass of the undissolved crystals can be obtained using a mass balance as follows, if Cb = 0 at t = 0where m0 is the initial mass of crystals. Combining eqs and 5,

For a diffusion-controlled process, a Noyes–Whitney type expression can be used to describe the dissolution processwhere rd is the diffusion rate, which equals to the mass transfer rate, and D is the diffusion coefficient. The diffusion layer thickness h is either a constant, when making the size-independent assumption (Figure A), or a function of L, the characteristic length of particles, when the size-dependent assumption is being made (Figure B), h = f(L).

Figure 2

Schematic of (A) size-independent and (B) size-dependent diffusion layer (yellow area) for spherical particles.

A number of dosage forms are used to disperse systems, and they can be classified according to the size of the domains of the dispersed (or internal) phase. There are three types: (1) coarse dispersions, >1 μm (2) colloidal dispersions, 1 nm–1 μm (3) molecular dispersions, <1 nm. Based on the distribution of particles sizes in the dispersed phase, coarse or colloidal dispersions can be further classified into monodisperse or polydisperse systems. In a monodisperse system, particles are of uniform size. In contrast, for a polydisperse system, the sizes of particles vary.[23] Ibuprofen crystals prepared in this work may be described as being a part of a coarse, polydisperse system, which must be taken into consideration when studying the dissolution kinetics and preparing suitable models.

As a part of the model development approach adopted, particles are initially considered as monodisperse spheres (i.e.: having the same single dimensional characteristic particle length). The model is then expanded to enable simulation of the particles’ polydisperse nature (distribution of particle size) and nonspherical geometries (multidimensional characteristic particle lengths). This polydisperse model is the model which is ultimately compared with experimental data.

Fundamental Dissolution Model

In a monodisperse system, all particles are considered to have the same characteristic length L. To determine the diffusion layer thickness for modeling of dissolution kinetics, the size dependency is illustrated as follows.

Size-independent Diffusion Layer Thickness

If the particles are assumed to have a size-independent diffusion layer thickness, thenwhere kd = D/h is the diffusion constant.where ρ is the actual density of the crystals, k is the volume shape factor, and N0 is the initial number of particles.

Assuming that N0 does not change with time, the surface area of undissolved crystals per unit volume of the solution at any time is given bywhere kA is the surface shape factor. For spherical crystals, kA = π and k = π/6.

Substituting eqs and 12 to eq yields

Size-dependent Diffusion Layer Thickness

For a size-dependent diffusion layer thickness model, it is assumed that there is a transition or critical particle length, Lc, bigger than which the particles are considered to dissolve with a constant diffusion layer thickness equal to (Lc/2). For particles with sizes smaller than Lc, the diffusion layer thickness is half of their own characteristic size. The selected method of diffusion layer thickness quantification, established previously in the literature,[15] which may be substituted into eq is provided in eq .

Model Expansion to a Polydisperse System

As a part of the dissolution of a polydisperse powder system, each particle size fraction could be considered as a monodisperse system, with all the particles in the relevant size bin dissolving at the same rate.[20]Figure shows the concept in a schematic form, in which the evolution of the particle size during dissolution is clearly illustrated.

Figure 3

Schematic of the dissolution of a polydisperse. L represents the characteristic length of seed crystals of i th particle size fraction at time t = j. Note that the number of particles in each particle size fraction is assumed to remain constant.

Substituting eqs , 12, and 7 to eq yields

In eq , hc is a constant diffusion layer thickness for size-independent situations, and hence, L is only a function of concentration. eq is a first-order ordinary differential equation, solved by numerical methods using a bespoke 4th Order Runge–Kutta solver developed in Matlab (to facilitate the simultaneous solution of eqs –25 at each time step). Initial conditions are taken directly from the experimental measurements corresponding to the scenario in the question being modeled.

The variables m, S, and C may be calculated at each time step (j) using L

Spherical vs Cylindrical Geometry

The shape and volume factors used to derive the preceding equations are based on a spherical geometry with constant values which do not change with time. In fact however, particles generally have nonspherical geometries. Because the dissolution rate is dependent on the surface area/volume ratio of any geometry, simulated dissolution rates of nonspherical particles using “uncorrected” particle characteristic sizes may be inaccurate when compared to experimental results. Introducing a cylindrical geometry characteristic into the model has the potential to describe particle shapes ranging from plates to needles by changing the shape factor.[15]

There are other potential advantages associated with assuming a cylindrical geometry. Particles with nonspherical geometries generally do not have a uniform diffusion layer thickness over their entire surfaces. Hence, dissolution may occur at different rates on various faces, giving rise to a change in the particle shape over time. A more rigorous model needs to consider the time dependency of the particle shape. By allowing more accurate estimations of the particle surface area, the performance of the shape simulation of dissolution kinetics will be improved. This being the case, a two-dimensional model should be introduced by taking two separate characteristic lengths of the particle into account.

As shown in Figure , for a single particle, the surface area S and volume V are given by

Figure 4

Schematic presentation of a cylindrical geometry particle.

Introducing the surface shape factor kA and the volume shape factor k to express S and V as functions of L only, we get

Substituting eq into eq yieldswhere L′/L is the aspect ratio of the particle.

Substituting eq into eq yields an expression for the shape volume factor

Similarly, if the surface area and volume are expressed as a function of L′ only

Substituting eq into eq yields

Substituting eq into eq yields

Dissolution rates are different on different faces of the particle, and hence, the particle shape will change over time, leading to a change of kA and k over the course of the dissolution process. As it is assumed that each face of the particle dissolves independently, the evolution of each dimension is given by eq . This being the case, two sets of the shape factors (kA and k, kA′ & k′) are needed in order to solve dL/dt and dL′/dt in eq . As the aspect ratio (L′/L) changes over time, the aspect ratio and hence all shape factors are updated on each time step as part of the numerical solution of eq .

Results and Discussion

Theoretical Simulation

In this part of simulation work, the CSDs of the particles are assumed to have a spherical geometry and a normal distribution, of which the population density f is calculated by eq . It is also assumed that no agglomerations or aggregations are present.where f is the population density, μ is mean particle size, and σ is standard deviation.

The CSDs of the spherical particles with (1) the same mean particle size (i.e. μ = 150 μm) but different standard deviations σ (i.e. 10, 20, 40, and 60) and (2) the same standard deviation (i.e. σ = 10) but different mean particle sizes (i.e. μ = 20, 50, 100, and 150 μm) are shown as number distributions in Figures A and 6A, respectively. The corresponding dissolution profiles are demonstrated visually in Figures B and, respectively. Equation was applied to simulate the dissolution, and D/h is assumed to be constant (i.e. D/h = 1), to simplify the calculation process.

Figure 5

Theoretical simulation (A) CSDs of spherical particles (μ = 150 μm and σ = 10, 20, 40, and 60); (B) corresponding dissolution profiles of ibuprofen.

Figure 6

Theoretical simulation (A) CSDs of spherical particles (σ = 10 and μ = 20, 50, 100, and 150 μm); (B) corresponding dissolution profiles of ibuprofen.

As illustrated in Figure , for particles within the same mean size (i.e. μ = 150 μm), it can be shown that the dissolution rate does not change significantly as the standard deviation σ increases from 10, 20, 40, to 60. Particles with a wider span dissolved slightly slower, as the overall surface area/volume ratio is smaller due to the presence of larger particles. Following this observation, it is proposed that the dissolution rate is less likely to be affected or determined by the span of the particles’ CSD for such particles.

However, in contrast as shown in Figure , for particles within the same standard deviation (i.e. σ = 10), the dissolution rate remarkably decreased as the value of μ varied from 20, 50, 100, to 150 μm. This leads to the conclusion that the mean of the particles’ CSD has a significant effect on the dissolution rate, which was previously shown and confirmed by published experimental dissolution work.[14]

To enable further investigation of these trends and to understand the relationship between these variables, Figure presents the simulated total dissolution times of particles with specific combinations of the mean particle size and standard deviation. Several different linear relationships may be observed in this figure. In general, the data can be separated into two distinct sets, displayed as hollow dots and solid dots in Figure .

Figure 7

Linear relationships between the mean size of particles μ and the total dissolution time t under different standard deviations σ (1, 5, 10, 20, 40, and 60).

The data with standard deviations (σ) of 1 (black), 5 (light blue), 10 (orange), and 20 (red), shown using hollow dots in Figure , represent particles with narrow distribution spans. All of these systems showed similar linear relationships, with slopes of 0.1975 (σ = 1, R2 = 0.99914), 0.1908 (σ = 5, R2 = 0.99955), 0.1796 (σ = 10, R2 = 0.99914), and 0.1609 (σ = 20, R2 = 0.99778). It may be observed here that the slope values decrease slightly with an increase in standard deviation.

The data sets illustrated using solid dots in Figure have wider distribution spans (with standard deviations of σ = 10, 20, 40, and 60). A similar linear relationship phenomena between the data sets is observed, with the systems having slopes of 0.0787 (σ = 10, R2 = 0.99618), 0.0521 (σ = 20, R2 = 0.99729), 0.0340 (σ = 40, R2 = 0.97114), and 0.0317 (σ = 60, R2 = 0.98606), respectively. Once more, the slopes have very similar values and decrease slightly with an increase in distribution standard deviation.

It may also be noted that two additional linear relationships are observed in Figure , namely, for particles with mean sizes of between 100 and 150 μm with standard deviations of 10 and 20. To account for this observation in this paper, particles with a standard deviation of 10–20 are treated as being a part of the first set (hollow markers in Figure ) if their mean size is less than 100 μm, while they are treated as being a part of the second set (filled markers in Figure ) if their size is above 100 μm.

As shown in Figure , groupings of linear trendlines are observed depending on the mean and standard deviation values of the distributions being dissolved. It is therefore possible upon further examination that more parallel linear relationships could be discovered. The identification of these linear relationships and further expansion of Figure may help to provide useful information when developing simplified correlations relating the total dissolution time to particle properties.

Determination of Diffusion Coefficient D

Since the 1950s, the diffusivities of dissolved solids in water have been correlated by various authors.[24−26] The correlations published by these authors have proven to be reasonably reliable as they yield similar results but none of them are identical.

Othmer–Thakar equation[24]

Scheibel equations:[25]

For V1 > V2

For V2 > V1

Wilke–Chang equation[26]

The Othmer–Thakar equation was subsequently redeveloped to propose new constants, and the revised Othmer–Thakar equation based on diffusivities of 87 different dissolved substances in dilute aqueous solution is shown below

In addition, a revised association parameter in the Wilke–Chang equation for water of 2.26 instead of 2.6 is likewise recommended. In these equations, the molecular properties of water are μ2 = 0.691, cp at 300 K, M2 = 18.015 g/mol, and V2 = 18.03 mL/mol and the molecular properties of IBU include M1 = 206.280 g/mol and V1 = 200.3 mL/mol.

The Stokes–Einstein relation has long been regarded as one of the hallmarks of transport in liquids. It predicts that the self-diffusion constant D is proportional to (τ/T)−1, where τ is the structural relaxation time and T is the temperature.

Stokes–Einstein relation[27]where k = 1.37 × 10–16, T = 300 K, v = 0.01 P, and s = 4.299 × 10–10 m.

The diffusion coefficient of IBU in water at 37 °C obtained by different methods is listed in Table . The averaged diffusion coefficient of IBU D = 8.88 × 10–6 cm2/s was used in the following simulations through this work.

Table 1

Diffusion Coefficient of IBU in Water at T = 37 °C

D, cm2/s (×106)
Model	Othmer–Thakar	Wilke–Chang	Scheibel	revised Othmer–Thakar	revised Wilke–Chang	Stokes–Einstein
Diffusivity	8.74	9.45	8.91	9.80	8.81	7.59
average D	8.88

Dissolution Rate Experimental Verification

To estimate the dissolution rate of ibuprofen in a phosphate buffer (pH = 7.20), the data collected from nine experiments previously published by Gao et al.[14] is used (the data’s main characteristics are listed in Table ) and the details can be found in Supporting Information. All of these experiments were operated at an undersaturation ratio of σ = −1, in which the driving force is very high, leading to rapid dissolution of the seed crystals. The temperature was kept constant at 37 °C. The experiments were performed using various agitation rates (from 100 to 150 rpm), and seed crystals with different CSDs and morphologies were used.

Table 2

List of Experiments Used to Estimate the Dissolution Kinetics of Ibuprofen in Phosphate Buffer (pH = 7.20) at T = 37 °Ca

data set #	morphology	L (μm)	D[4,3] (μm)	R (rpm)	–σ	fit./val.
1	polyhedral	<150	130	100	1	val.
2	polyhedral	150–300	287	100	1	fit.
3	polyhedral	300–500	462	100	1	val.
4	polyhedral	500–850	631	100	1	fit.
5	polyhedral	150–300	287	150	1	val.
6	polyhedral	300–500	462	150	1	fit.
7	polyhedral	500–850	631	150	1	val.
8	polyhedral	>850	1322	150	1	fit.
10	needle-like	150–300	181	100	1	val.

L is the sieved size fraction; D[4,3] is the volume weighted mean particle size; R is the agitation rate; −σ is the undersaturation ratio. In the column fit./val., fit. indicates whether an experiment was used for fitting the parameters of the dissolution kinetics in this paper or just for validation. Val. indicates that data was only used for model validation, while fit indicates that the data was used for model fitting.

Data sets 1–8 and 10 listed in Table correspond to the data collected from dissolution experiments previously published.[14] The experimental protocols as well as an explanation of the experimental set-up and procedure and the equipment employed for characterizing both the solid and liquid phase are described in detail in the original publication.

The dissolution profiles and the pre- and postdissolution CSDs from the experimental data sets 2, 4, 6, and 8 (covering 2 different hydrodynamic mixing conditions: 100 and 150 rpm) were used in this paper in order to regress the diffusion layer thickness. These data sets were selected for this regression in order to ensure that a wide range of hydrodynamic mixing conditions were considered due to the fact that the agitation rate is likely to be an important factor when determining the diffusion rate.[14] It should also be noted that this data used in the regression consisted only of crystal samples with polyhedral morphologies, which involve relatively smaller surface area/volume aspect ratios compared to needle-like crystals. The remaining data listed in Table was retained for model validation.

Data sets 1, 3, 5, and 7 were generated using similar conditions to data sets 2, 4, 6, and 8, respectively but with differing CSDs. Data set 10 uses needle-like seed crystals.

Finally, with the fitting parameters from the dissolution kinetics estimation, the model could be extended to predict the seed crystals with different conditions under the similar hydrodynamics, which will be discussed in detail in Section .

Determination of Diffusion Layer Thickness

The diffusion layer thickness h was determined using the following approach. First, an initial critical particle length, Lc, was estimated from the literature.[15,20] The diffusion layer thickness, h, in eq was then calculated from this Lc value using eq (e.g.: it was set equal to Lc/2 when the size-independent assumption is being applied). The CSD evolution during the dissolution process was then simulated using eq with initial conditions taken directly from the experimental predissolution data. Using the CSD at each time point, the total mass of remaining crystals m and the concentration of the bulk solution Cb were calculated through eqs –21, respectively. The calculated concentrations of the bulk solution Cbcal were then directly compared with experimental values Cbexp at each time point in order to calculate the total sum square error SSE for the data setwhere n is the number of experimental data points.

The diffusion layer thickness h was then optimized in order to minimize this SSE using the inbuilt “fminsearch” in Matlab. Four of the total nine datasets were used for this regression/optimization. The optimized h values were then validated using the remaining five dataset (data sets are labeled as being used for either regression or validation in Table ). This validation is discussed in further detail in Section .

Diffusion layer thicknesses h at agitation rates of 150 rpm were determined by the method discussed above (using both spherical and cylindrical geometry assumptions), with respect to data sets 6 and 8 as listed in Table and shown in Figure . Figure demonstrates good agreement between the regressed and experimental data.

Table 3

Regression of h With the Dissolution Kinetic Model (R = 150 rpm)

experimental condition	geometry	size dependency	h (μm)	time (s)	SSE x10–20
run 6: L = 300–500 μm	spherical	independent	32.87	1136	19.06
D[4,3] = 462 μm		dependent	32.74	1137	13.93
R = 150 rpm	cylindrical	independent	25.57	1115	11.96
morphology = polyhedral		dependent	25.71	1116	2.057
run 8: L > 850 μm	spherical	independent	30.65	2248	40.93
D[4,3] = 1322 μm		dependent	30.66	2250	20.60
R = 150 rpm	cylindrical	independent	23.65	2214	30.93
morphology = polyhedral		dependent	23.84	2217	9.636

Figure 8

Experimental dissolution profiles of fine (o) and coarse (Δ) ibuprofen when R = 150 rpm. Simulated curves were drawn using spherical geometry with (A) size-independent h and (B) size-dependent h, and cylindrical geometry with (C) size-independent h and (D) size-dependent h.

In the spherical geometry simulations, the values of the size-independent diffusion layer thickness were found to be h = 32.87 μm and h = 30.65 μm, with sum square errors (SSEs) of 1.906 × 10–19 and 4.093 × 10–19, respectively (Table ). The values of the size-dependent diffusion layer thickness regressed from the same two data sets were found to be quite similar at h = 32.74 μm and h = 30.66 μm (SSE = 1.393 × 10–19 and 2.060 × 10–19). SSE values are very low in both scenarios showing good agreement between experimental and regressed results. However, the SSEs of the latter simulation are slightly smaller, indicating that a size-dependent diffusion layer thickness gives marginally better modeling results. These same trends are also demonstrated in the cylindrical geometry dissolution simulation (listed in Table ), where the values of the size-independent diffusion layer thicknesses are 25.57 and 23.65 μm, with SSE values of 1.196 × 10–19 and 3.093 × 10–19, respectively, while in comparison, the values of size-dependent diffusion layer thicknesses are 25.71 and 23.84 μm with corresponding SSE values of 2.057 × 10–20 and 9.636 × 10–20.

The total dissolution time is found to be shorter when using the cylindrical geometry compared to the corresponding spherical geometry assumption. This is expected due to the fact that the dissolution rate is dependent on the surface area according to eq , and a sphere has the smallest surface area/volume ratio of any geometry, giving it a slower dissolution rate and a longer total dissolution time compared with cylindrical geometry particles.

It can be also seen in Table that the simulated dissolution rate when using h = f(L,t) is similar to that obtained when holding h constant. In addition, both approaches achieved extremely low SSE values which are very similar in magnitude. Hence, the size-independent diffusion layer thickness was averaged to be (25.57 + 23.65)/2 = 24.61 μm and the size-dependent diffusion layer thickness was averaged to be (25.71 + 23.84)/2 = 24.78 μm. For larger particles, the dissolution rate is slower and h also changes more slowly, leading to a relatively constant h despite changes of L over time in the process h = (f,t). The diffusion layer thickness was therefore considered to be approximately uniform equal to 25 μm when R = 150 rpm for the remainder of this work.

The same procedure as described above was used to regress the diffusion layer thickness at an agitation rate equal to 100 rpm. The regression was carried out by fitting the model with data sets 2 and 4, as listed in Table (also displayed graphically in Supporting Information Figure S2). Compared with the results obtained at R = 150 rpm, the corresponding diffusion layer thickness at R = 100 rpm is thicker due to the slower diffusion rate. Simulation using the cylindrical geometry assumption once more gives slightly better regression results than those based on a spherical geometry. This result may be justified by examining SEM images (Supporting Information Figure S1) of the crystals being dissolved in this analysis. The crystals have needle and polyhedral morphologies which are described better by the cylindrical modeling approach adopted compared to using a spherical assumption. The size-independent diffusion layer thickness had an average value of (33.08 + 33.65)/2 = 33.37 μm, and the size-dependent diffusion layer thickness had an average value of (32.75 + 33.84)/2 = 33.30 μm. In line with the assumption made as a part of the R = 150 rpm analysis, the diffusion layer thickness h at R = 100 rpm is approximated to be equal to 33 μm.

Table 4

Regression of h With the Dissolution Kinetic Model (R = 100 rpm)

experimental condition	geometry	size dependency	h (μm)	time (s)	SSE x10–20
run 2: L = 150–300 μm	spherical	independent	42.53	948	54.48
D[4,3] = 287 μm		dependent	42.06	946	24.33
R = 100 rpm	cylindrical	independent	33.08	935	31.96
morphology = polyhedral		dependent	32.75	930	8.036
run 4: L = 500–850 μm	spherical	independent	41.79	2141	40.93
D[4,3] = 631 μm		dependent	41.84	2134	20.60
R = 100 rpm	cylindrical	independent	33.65	2015	30.93
morphology = polyhedral		dependent	33.84	2007	9.636

Prediction of the Postdissolution CSD

With the regression of the dissolution layer thickness h, the CSDs of the postdissolution process could be recalculated, with a typical example shown in Figure . The CSDs were plotted as volume (%) v.s. particle size, which is consistent with the experimental results measured by experimental equipment used. The population density n of the simulation is calculated by

Figure 9

Comparison between experimental and simulated postdissolution CSDs of fine particles.

It may be noted that in the experimental study,[14] there is a slight right skew of the CSD of ibuprofen after dissolution as displayed in Figure (blue dashed line). This is also confirmed using the CSD details of ibuprofen pre- and postdissolution (Supporting Information Table S4), showing that the mean particle size gets larger following the dissolution process. This may be explained by the experimental design adopted and the Noyes–Whitney’s theory. In the experimental procedure, 1 g of ibuprofen crystal seeds were dispersed in the solution. However, only about half of the amount could dissolve due to solubility limits under experimental conditions. The purpose of this approach is to allow examination of particle size evolution, and hence, sufficient crystals must remain following dissolution in order to measure the CSD (using the Malvern Mastersizer). According to the Noyes–Whitney equation (eq ), smaller particles dissolve at a faster rate compared to larger particles, as they have a greater surface area per unit volume. In this case, as an excess amount of ibuprofen is added into the system and as the smaller particles dissolve at a faster rate than larger particles, the particle size distribution demonstrates a slight right skew. In addition, it is observed that the height of the distribution increases after dissolution in Figure . This may be explained as follows: the y-axis of the figure is presented as a normalized volume in which the total volume fraction is equal to 1 (or 100%). Therefore, the height of the distribution represents the volume fraction of the corresponding particle size and the higher height demonstrates that the percentage of relative larger particles increases as smaller particles dissolve faster than the larger ones.

For the prediction of the postdissolution CSD, the simulated curve (red dash line in Figure ) tends to slightly shift to the right, exhibiting good qualitative agreement between experimental data and model prediction.

Validation of the Dissolution Model

The predictive capabilities of the model were tested by using the validation experiments listed in Table . The diffusion layer thicknesses h are estimated to be h = 33 μm (R = 100 rpm) and h = 25 μm (R = 150 rpm), respectively, as outlined previously. The values have been applied to the model using the cylindrical geometry approach, using both size-independent and size-dependent assumptions. The fitted outputs are dissolution profiles (IBU concentration vs dissolution time) and postdissolution CSDs (volume vs particle size).

The validation involves different hydrodynamic mixing conditions (agitation rate), seed CSDs, and morphologies to cover a wide range of the experimental conditions, as contained in data sets 1, 3, 5, 7, and 10 (i.e.: the data sets which were not used for the initial regressions). The dissolution profiles together with the pre- and postdissolution CSDs of these experiments and the corresponding model predictions could be found in Supporting Information Figures S3–S7.

Overall, there is a good agreement between the experimental data and simulated results as shown in Table . A slight consistent overprediction of the postdissolution CSDs exists which may be explained by the fact that the simulated CSDs are monomodal distributions; however, the corresponding experimental CSDs contain a slight bimodal shape. Such errors are to be expected as the generic cylindrical particle model utilized in this work is only an approximation of the true particle shape. In addition, SEM images (Supporting Information Figure S1) from the experimental data in the published work show some aggregations or agglomerations in the particles after crystallization,[14] which will induce some degree of error when comparing simulated and experimental results.

Table 5

List of the Validation of the Dissolution Kinetic Model

data set	size dependency	aspect ratio	h (μm)	time (s)	SSEx10–8
1	independent	3	33	973	9.508
	dependent			952	7.230
3	independent	3	33	1175	0.524
	dependent			1162	0.267
5	independent	3	25	905	8.653
	dependent			902	6.248
7	independent	3	25	>2000	0.025
	dependent			>2000	0.016
10	independent	10	33	1070	6.303
	dependent			974	4.055

It may be observed that the simulated results align closely with experimental results at different agitation speeds, confirming the model’s applicability for predicting the dissolution behaviors of ibuprofen particles.

Prediction with the Dissolution Model

Using the regressed diffusion layer thickness h, the diffusion-based dissolution model could be applied widely to predict the dissolution behavior of ibuprofen particles under different conditions. The model is used to predict the dissolution profiles with different initial total seed crystal mass present, as shown in Table and Figure . It may be observed that the dissolution rate is higher with more ibuprofen particles initially introduced into the tank, leading to a shorter total dissolution time. This is because the greater the mass of initial seed crystals, the larger the total surface area/volume in the solvent system, leading to a quicker dissolution kinetics behavior.

Table 6

Prediction of the Dissolution Kinetics With Different Initial Dissolution Mass

	total dissolution time (s)
dissolution mass	R = 100 rpm, h = 33 μm	R = 150 rpm, h = 25 μm
m = 0.5	1758	1532
m = 0.75	1254	1051
m = 1.0	935	780
m = 2.0	480	386
m = 3.0	327	261

Figure 10

Prediction of the dissolution kinetics with different initial dissolution mass: (A) R = 100 rpm; (B) R = 150 rpm.

Simulated conditions

L = 150–300 μm, R = 100 rpm, initial total dissolution mass = 1 g, morphology = polyhedral, length ratio = 3, h = 33 μm.

L = 150–300 μm, R = 150 rpm, initial total dissolution mass = 1 g, morphology = polyhedral, length ratio = 3, h = 25 μm.

Conclusions

A diffusion-based model, developed from the Noyes–Whitney equation in which the nonsink dissolution term is a function of the remaining surface area and the concentration gradient across the boundary layer, has been developed to simulate dissolution kinetics of a poorly water-soluble drug, ibuprofen. Dissolution parameters include initial particle size, diffusion coefficient, density, and diffusion layer thickness. The model was established by using spherical geometry (single dimensional characteristic particle length) and cylindrical geometry (multidimensional characteristic particle lengths) assumptions. Data from previously published work was fitted to the model to obtain the diffusion layer thickness, considering both size-dependent or size-independent scenarios. The model was then validated using unused dissolution experimental data, demonstrating good agreement between experimental and predicted values. Postdissolution CSDs were also predicted using the model, once more demonstrating reasonable agreement with experimental data, with a slight overprediction of the distribution.