A new modelling approach for dissolution of polydisperse powders

A new modelling approach for dissolution of polydisperse powders is developed within the framework of the classical Noyes – Whitney/Nernst – Brunner analysis. Its distinguishing feature is that the underlying continuous particle-size distribution is retained. Two different but related dependencies of the diffusion-layer thickness on particle size are considered. First, a power-law dependence that interpolates between a thickness that is proportional to (or equals) the particle radius (obtained when the exponent equals 1) and a constant thickness (obtained when the exponent is 0). Second, a piecewise linear function such that the thickness equals the particle radius for sufficiently small particles and is constant for larger ones. The modelling approach is exemplified by consideration of a lognormal particle-size distribution. Highly accurate closed-form expressions for the fraction of dissolved drug are obtained for dissolution under sink conditions (which are exact if the diffusion-layer thickness is radius-independent). Moreover, it is demonstrated that any result derived under sink conditions can be reused to determine the fraction of dissolved/absorbed drug under non-sink conditions, using the concept of a retarded time. Comparison with literature data and experiments are used to validate the modelling approach and to demonstrate its usefulness in a practical context.


Introduction
Drug solubility and dissolution are critically important for oral bioavailability of drugs (Amidon et al., 1995) and predictive dissolution modelling therefore plays a central role during formulation, development and testing of pharmaceutical products intended for oral delivery (Zaborenko et al., 2019).Dissolution is highly relevant for other delivery routes as well, since delivery of solid drug requires a dissolution step before drug can be transported across epithelia.A notable example is inhaled drug delivery, where dissolution can dictate the bioavailability and pharmacokinetics of drugs administered in solid form (Hastedt et al., 2022).
Dissolution is generally (but not universally; see e.g. the article by Shekunov ( 2017)) considered to be diffusion-controlled (Siepmann and Siepmann, 2013).The dissolution rate can then be rationalized in terms of an effective hydrodynamic diffusion layer (for brevity henceforth referred to as the diffusion layer).The modifications of the classical Noyes-Whitney equation (Noyes and Whitney, 1897) proposed by Nernst (1904) and Brunner (1904) are consistent with this view.The diffusion-layer thickness depends on the hydrodynamic conditions in the vicinity of the dissolving particle(s).There is also a geometric factor, which is important especially for small (≲10μm) particles, as the analysis of Wang andFlanagan (1999, 2002) demonstrates.Different models of the diffusion-layer thickness have been compared by Sugano (2008).
It is now widely recognized that polydispersity needs to be considered in order to obtain reliable estimates of dissolution rate, both in an oral (Johnson and Swindell, 1996) and inhaled setting (Amini et al., 2021;May et al., 2014).The most widely used approach, pioneered by Johnson and co-workers (Hintz and Johnson, 1989;Johnson, 2003;Johnson and Swindell, 1996;Lu et al., 1993), relies on binning the underlying continuous particle-size distribution.In other words, the inherently continuous distribution is replaced by a discrete one.An alternative approach, mostly utilized in chemical engineering, relies on population balance modelling as proposed by Leblanc andFogler (1987, 1989) and elaborated further upon by Giona et al. (2002).In addition, hierarchical models with confinement have recently been proposed (Salehi et al., 2020;Wang et al., 2015).Although this analysis is elegant, it introduces additional complexity that precludes an analytical treatment.
The purpose of this article is to demonstrate that the modelling approach of Johnson and co-workers (Hintz and Johnson, 1989;Johnson, 2003;Johnson and Swindell, 1996;Lu et al., 1993) can be readily generalized to continuous particle-size distributions without recourse to population-balance modelling.This article thus generalizes and further develops the results obtained in our previous work (van der Zwaan et al., 2022).Highly accurate (albeit generally approximate) analytical results are obtained for lognormal particle-size distributions when particles dissolve under sink conditions.For dissolution under non-sink condition, the fraction of dissolved/absorbed drug is obtained as the solution of a single nonlinear first-order ordinary differential equation.

Starting points
We assume that dissolution of monodisperse powders can be sufficiently well approximated by the Noyes-Whitney equation (Brunner, 1904;Nernst, 1904;Noyes and Whitney, 1897;Siepmann and Siepmann, 2013), here written as where the minus sign is included because M represents the mass of drug that remains in solid form.The remaining symbols have their usual meaning, i.e.D is the diffusion coefficient of the drug in the dissolution medium, h HDL is the thickness of the hydrodynamic diffusion layer (diffusion layer for short), A is the surface area of the solid-liquid interface, C s is the solubility and C is the bulk concentration of the drug.We will consider two different but related dependencies of h HDL on particle size (radius) r: A power-law dependence of the form h HDL ∝r α , where α is a constant, and a piecewise linear function of the form h HDL = min(r, h), where h is the maximal diffusion-layer thickness.The powerlaw dependence was assumed for mathematical convenience rather than necessity.It naturally interpolates between a radius-independent diffusion-layer thickness, obtained for α = 0, and a thickness that is pro-portional to (or equals) particle size (radius), obtained for α = 1.The piecewise linear function approximates the dependence proposed by Wang andFlanagan (1999, 2002) and results in comparable dissolution profiles (Johnson, 2012).It is reasonable to assume that transport occurs through a nearly stagnant layer when particle size is sufficiently small, implying that h HDL = r for spherical particles (Manca and Rovaglio, 2003).Moreover, experimental evidence suggests that the dependence of h HDL on particle size becomes less pronounced for larger particle sizes; see Andersson et al. (2022) and references therein.All particles are assumed to have a similar shape that is retained during dissolution, so that geometric arguments can be used to relate the particle size to the particle area and volume/mass.Thus, we assume that A∝r 2 and M∝r 3 .This argument is independent of particle shape (as long as it does not change), but we will nevertheless for convenience refer to r as the radius.
All drug is considered to be present in solid form initially.

Monodisperse powders, sink conditions
For monodisperse powders dissolving under sink conditions (i.e.C = 0), the assumed proportionalities A∝r 2 and M∝r 3 and the Noyes--Whitney Eq. (1) give When h HDL ∝r α , Eq. ( 2) takes the form where the proportionality constant has been written as K ′ α /(1 +α) for convenience.Eq. ( 3) can be rewritten as and with R the initial particle size (radius), the solution can be immediately obtained as (5) In particular, α = 0 (h HDL = constant) gives a linear decrease, i.e. r = R − K ′ 0 t.For α = 1 (h HDL ∝r), a square-root dependence is obtained, i.e.
The time t d required for complete dissolution is obtained from Eq. ( 5) as These results are illustrated in Fig. 1 for the particular case that the initial particle size (radius) R = h.To enable an easy comparison between the curves, we have let K ′ 0 = K/h and K ′ 1 = 2K; see the discussion following Eq.( 7) below.
When h HDL = min(r, h), Eq. ( 2) takes the form where K is a dissolution-rate constant.If the initial particle size (radius) R ≤ h, the solution is obtained from Eq. ( 5) with α = 1 and K ′ 1 = 2K.For R > h, the initial decrease in r is obtained from Eq. ( 5) with α = 0 and Thereafter, the decrease in r is obtained from Eq. ( 5) with α = 1 and the following substitutions: , R→h and t→t − t c .In summary, the solution of Eq. ( 7) is obtained as This behaviour is illustrated in Fig. 2, which displays r/h vs Kt/h 2 for a number of values of R/h between 0.5 and 3.5.
To determine the time required for complete dissolution, we note that Eq. ( 6) is applicable for R ≤ h (with α = 1 and K ′ 1 = 2K).For R > h, we reach r = h at t = t c , and the remaining time is obtained from Eq. ( 6) In summary, we obtain In any case, the fraction of drug that remains in solid form is obtained as where M 0 is the initial mass of drug and r is given by Eq. ( 5) or Eq. ( 8) depending on the choice of model.The fraction of drug that remains in solid form for a radius-independent diffusion-layer thickness (h HDL = h) and a radius-dependent thickness (h HDL = r) are compared in Fig. 3.A polynomial approximation of r 3 , derived in Appendix A, is also illustrated in the figure.

Polydisperse powders, sink conditions
We next consider a polydisperse powder with a (number) particlesize distribution f(R), dissolving under sink conditions.When sink conditions prevail, each particle dissolves independently and the fraction of solid drug can be determined by superposition.For a given time t, particles with initial size R ≤ R d have dissolved, where R d increases with time.Hence, R d represents the initial radius for those particles that become completely dissolved at time t.Specifically, the fraction of solid drug is obtained as where is the third moment of the particle size distribution.Eq. ( 11) was stated in slightly different form in our prior work (van der Zwaan et al., 2022) and a related expression appears in the work by Leblanc and Fogler (1987); cf.their Eq.( 25).

Non-sink conditions: Mass balance
To address dissolution under non-sink conditions, we consider the situation depicted in Fig. 4, in which dissolution occurs in a donor compartment with volume V d .Dissolved drug enters a receptor compartment (i.e. is absorbed) at a rate proportional to the concentration in the donor compartment.Taking absorption into account, mass balance can be expressed as where K ′ abs is a mass-transfer coefficient.It proves convenient to use the initial "concentration" S 0 = M 0 /V d of solid drug in the donor compartment (remember that M 0 is the initial mass of drug) to write the dependent variables in non-dimensional form.Letting c = C/S 0 and noting that M/M 0 = s, the fraction of drug that remains in solid form, Eq. ( 15) becomes where K abs = K ′ abs /V d is an absorption-rate constant.We remark that 1/K abs was referred to as the diffusion and/or permeation time and denoted by t diff in our prior work (Frenning et al., 2020;van der Zwaan et al., 2022).Noting that K ′ abs C represents the absorption rate, the fraction of absorbed drug is obtained as Using Eq. ( 17), Eq. ( 16) can be integrated to produce where the initial values c(0) = 0 and s(0) = 1 have been used.Eq. ( 18) could also be obtained from the fact that all fractions must add to unity, c + s + u = 1, noting that c = (1/K abs )du/dt, as Eq. ( 17) demonstrates.This type of model reduction has previously been used by us (Frenning et al., 2005;van der Zwaan et al., 2022) and others (Butcher et al., 2008).It can be noted, however, that the source term 1 − s needs to be specified before a solution of Eq. ( 18) can be obtained.This is done in the following section(s).

Non-sink conditions: Retarded time
For dissolution under non-sink conditions, the analogue of Eq. ( 4) takes the form where c s = C s /S 0 could be regarded as a non-dimensional solubility or an inverse dose number (Hastedt et al., 2022;Oh et al., 1993).The derivation of Eq. ( 19) parallels that of Eq. ( 4), but C s − C is substituted for C s in Eq. (2) since C is not set equal to 0. An analogous result would be obtained for other dependencies of h HDL on r, such as the piecewise linear one investigated in this work.Eq. ( 19) demonstrates that the change in particle size caused by dissolution under non-sink conditions during a small time interval dt is the same as the change caused by dissolution under sink conditions during a time interval dt ′ given by This observation allows us to define a retarded time t ′ as where Eq. ( 17) has been used.Hence results obtained for sink conditions do, when evaluated at the retarded time, correspond to those for nonsink conditions, i.e.
A single nonlinear ODE results if this expression is inserted in Eq. ( 18).This result was derived in our previous work (van der Zwaan et al., 2022), but only for a radius-independent diffusion-layer thickness.

Discrete distribution
To illustrate the procedure and to verify the results obtained, we first consider a discrete distribution of the form where f i is the number fraction in bin i and δ(R − R i ) is the Dirac delta function, which makes it possible to express a discrete distribution in continuous form (see, e.g., Arfken et al., 2013).For simplicity we here assume that h HDL ∝r α so that the particle size (radius) is obtained from Eq. ( 5).Inserting Eq. ( 23) in Eq. ( 11) and using the properties of the Dirac delta function, we obtain Here μ 3 = ∑ n i=1 f i R 3 i and the summation is to be taken over size fractions that have not already dissolved, i.e. for which ) 1/(1+α) according to Eq. ( 13).

Lognormal distribution
As a more realistic example, we consider a polydisperse powder with a lognormal particle-size distribution, where R 1 is a scale parameter (the median) and σ is a shape parameter (the natural logarithm of the geometric standard deviation of the distribution).We first assume that h HDL ∝r α and that α = 0. Hence Eq. ( 5) where the binomial theorem has been used to expand the cube of The fraction of drug remaining in solid form is calculated via Eq.( 11), i.e. as the integral of r 3 f(R) from R d to ∞, normalised by μ 3 .
It proves convenient to introduce the function where β is a constant, is the β th moment of the lognormal distribution and erfc(•) is the complementary error function (Abramowitz and Stegun, 1965).With r 3 given by Eq. ( 26), r 3 f(R) can be readily integrated with the aid of Eq. ( 27).Since R d = K ′ 0 t for α = 0, c.f. Eq. ( 13), we obtain the exact result which was obtained in a slightly different form in our previous work (van der Zwaan et al., 2022).
Finally, we assume that h HDL = min(r,h).It then proves convenient to subdivide the integral of r 3 f(R) into three pieces, one for each instance of Eq. ( 8).Each integral can be (approximately) evaluated as described above for h HDL ∝r α (cf.Appendix B).The final result takes the form where The coefficients b i are given by For n = 5, expanded expressions are provided in Appendix B.
To illustrate these results, model calculations were performed for two different particle-size ranges.One of these has a mass-median diameter (MMD) of 2μm that is typical for inhaled medicines.The other has an MMD of 50μm as would be relevant for oral drugs.In both cases, a range of standard deviations σ between 0.1 and 1.0 were investigated.The resulting volume-based particles-size distributions are illustrated in Fig. 5a and b (notice that the volume-based median R 1V = MMD/2).To convert to number-based particle-size distributions, used in the calculations, we use the Hatch-Choate conversion equation (Heintzenberg, 1994) Model calculations were performed for the piecewise linear model, h HDL = min(r, h), with h = 10μm (Sugano, 2008) to yield the results displayed in Fig. 6.A monodisperse powder would, according to Eq. ( 9), have completely dissolved for Kt/h 2 = 0.005 in (a) and for Kt/h 2 = 2 in (b).It is evident that polydispersity can significantly prolong this time, especially for small particles with a radius-dependent diffusion-layer thickness.Johnson and Swindell (1996) have presented simulation data (fraction absorbed after 6h) obtained for three different absorption-rate constants (0.001, 0.01 and 0.1min − 1 ), four different values of solubility (0.001, 0.01, 0.1 and 1.0g/cm 3 ), four doses (1, 10, 100 and 250mg) and four particle size distributions (with MMDs of 10, 25, 50 and 100μm).Assuming a spherical particle shape, the dissolution rateconstant was calculated as K = DC s /ρ where D = 5 × 10 − 6 cm 2 /s is the diffusion coefficient, C s is the solubility (as before) and ρ = 1.3g/cm 3 is particle density.The thickness of the hydrodynamic diffusion layer was set equal to particle size for radii smaller than h = 30μm and equal to h otherwise.The dissolution volume V d = 250mL.

Comparison with mathematical models from literature
We first consider the smallest particle size (MMD = 10μm), for which r < h for the vast majority of the particles.Hence we assume that h HDL ∝r α with α = 1, as was previously done by Butcher et al. (2008).
The fraction of solid drug remaining under sink conditions is then given by Eq. ( 24) with α = 1 and K ′ 1 = 2K.Specifically, the volume-based discrete distribution used by Johnson and Swindell (1996) was employed (Table 1).To convert from a volume-based to a number-based distribution, we used the fact that The fraction of absorbed drug after 6h was calculated via Eqs.( 18) and ( 22).As shown in Fig. 7, we do for all absorption rates and values of solubility obtain nearly identical results as those reported by Johnson and Swindell (1996).This is not surprising, since our approach in this case is equivalent to the one used by Butcher et al. (2008) who also obtained virtually identical results.
We next consider lognormal particle-size distributions with the same parameters as used by Johnson and Swindell (1996).Their geometrical standard deviation of 2.0 corresponds to our σ = ln(2.0)≈ 0.693.The Hatch-Choate conversion equation (Heintzenberg, 1994) is used to obtain R 1 = MMDe − 3σ 2 /2.As seen in Fig. 8a, the results obtained for a MMD of 10μm are virtually identical with those obtained by Johnson and Swindell (1996).Similar but not identical results are obtained for the other particle sizes as well (Fig. 8b -d).In particular, our results indicate a somewhat slower absorption for intermediate and high absorption rates, especially for a solubility of 0.1g/cm 3 (blue and grey triangles in Fig. 8c and d).The reason for this is that Johnson and Swindell (1996) considered the diffusion-layer thickness to be proportional to particle size also for large particles, but adjusted the dissolution-rate constant so that the initial dissolution rate would be correct (K eff = KR/h was substituted for K, which results in an initial dissolution rate that is proportional to K eff /R = K/h).However, as a consequence, the dissolution rate will be overestimated for h < r < R, since K eff /r = (K/h)(R/r)〉K/h.This effect is only seen when R > h and it is most evident when dissolution rather than absorption is rate limiting, as for the mentioned parameter values.

Comparison with experimental data from literature
To demonstrate the usefulness of the modelling approach in an applied context, we consider dissolution data obtained for three inhaled drug products containing budesonide (BUD), beclomethasone dipropionate (BDP) and fluticasone propionate (FP) (van der Zwaan et al., 2022).In vitro dissolution was assessed by a Transwell dissolution method according to the procedure described elsewhere (Franek et al., 2018).Phosphate-buffered saline (PBS) with 0.5 % sodium dodecyl sulphate (SDS) was used as dissolution medium.Parameters characterising the particle-size distribution (MMD and Span), transport across the Transwell membrane (K abs ) and solubility (c s ) were determined in separate experiments (Table 2).The Span is defined as the difference between the 90th and 10th percentile of the particle-size distribution,

Table 1
Discrete particle-size distribution used by Johnson and Swindell (1996).For each size fraction i, the particle diameter 2R i and volume (or mass) fraction ( f V ) i are given.where z 0.1 ≈ 1.28 is the 0.1th quantile of the standard normal distribution (see Appendix C).
The dissolution-rate constant K was determined by fitting the theoretical fraction absorbed drug to the experimental data (Fig. 9).For simplicity, it was assumed that h HDL ∝r α , with α = 1, which is justified since very few large particles exist.Hence, the fractions dissolved under Fig. 7. Comparison between our results (fraction absorbed after 6h) for a discrete distribution and those reported by Johnson and Swindell (1996) for a mass-median particle size of 10μm.Absorption-rate constants are distinguished by colour.0.001 (red), 0.01 (blue) and 0.1min − 1 (grey) and values of solubility by symbol shape: 0.001 (circles), 0.01 (squares), 0.1 (triangles) and 1.0g/cm 3 (diamonds) (colour online).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 8.Comparison between our results (fraction absorbed after 6h) for continuous lognormal distributions and those reported by Johnson and Swindell (1996) for the indicated mass-median diameters.Absorption-rate constants are distinguished by colour: 0.001 (red), 0.01 (blue) and 0.1min − 1 (grey) and values of solubility by symbol shape: 0.001 (circles), 0.01 (squares), 0.1 (triangles) and 1.0g/cm 3 (diamonds) (colour online).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)non-sink (sink) conditions were obtained from Eq. ( 29) evaluated at the retarded time t ′ (normal time t).To enable an easy comparison with the results presented in our previous work (van der Zwaan et al., 2022), we also report a characteristic time for dissolution, denoted t diss and defined so that the initial dissolution rate equals M 0 /t diss .This yields ( ds dt . In particular, for α = 0 the coefficient a 1 = − 3 so that simplification with the aid of Eq. ( 28) results in ) . As seen in Table 2, t diss is about 3 times larger for FP than for the other drugs.

Conclusions
This work has expanded upon a new modelling approach for dissolution of polydisperse powders, developed within the framework of the classical Noyes-Whitney equation.Its main novelty lies in the fact that the underlying continuous particle-size distribution is retained, contrary to the most commonly used methods.The analysis has been performed for two different but related dependencies of the diffusion-layer thickness on particle size.First, a power-law dependence and, second, a piecewise linear function such that the thickness equals the radius for sufficiently small particles and a constant for larger ones.A lognormal particle-size distribution was considered.For dissolution under sink conditions, highly accurate closed-form expressions were obtained for the fraction of dissolved drug, which are exact if the diffusion-layer thickness is radius-independent.It is demonstrated that results derived under sink conditions can be reused to determine the fraction dissolved and the fraction absorbed drug under non-sink conditions, using the concept of the retarded time.However, the solution is in this case obtained from a nonlinear first-order ordinary differential equation that needs to be solved numerically.The validity and usefulness of the modelling approach has been demonstrated by comparison with literature data and experiments.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix B
Fraction absorbed for the lognormal particle-size distribution.To evaluate the integral that results when Eq. ( 5) is inserted into Eq.( 11), first note that where x = K ′ α t/R 1+α is a non-dimensional quantity with a value between 0 and 1.The approximate Eq. (A1) thus yields After multiplication with f(R), this expression can be integrated term-wise with the aid of Eq. ( 27).Noting that R d = ( K ′ α t ) 1/(1+α) , cf.Eq. ( 13), we obtain Eq. (30) in the main text.When h HDL = min(r,h), it proves convenient to subdivide the integrals into three pieces, one for each instance of Eq. ( 8), as follows: Here, R c = h +Kt/h represents the initial size (radius) of those particles for which the dissolution kinetics changes at time t (from a radiusindependent to a radius-dependent h HDL ).
To evaluate the integral appearing in Eq. (B3), we proceed in the same manner as in the derivation of Eq. ( 30).Making the substitutions α = 1 and K ′ 1 = 2K and noting that integration is from min(R d , h) to h, we obtain Eq. ( 32) in the main text.To evaluate the integral appearing in Eq. (B4), we first note that where x = 2(R c − R)/h is a non-dimensional quantity with a value between 0 and 1. Letting X = 2R/h and Y = 2R c /h, the approximate Eq. (A1) thus yields where the binomial theorem has been used.Changing the order of summation, the above expression may be restated as

Fig. 1 .
Fig. 1.Illustration of the decrease in particle radius for h HDL ∝r α (α = 0 and 1) and R = h.When R = h, dissolution is twice as fast for α = 1 than for α = 0.

Fig. 2 .Fig. 3 .
Fig. 2. Decrease in particle radius for h HDL = min(r, h) for a number of initial radii R/h between 0.5 and 3.5.The dashed red line highlights the change in kinetics that occurs when r/h = 1.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. Dissolution under non-sink conditions in a donor compartment followed by absorption across a membrane into an acceptor compartment.

Fig. 6 .
Fig. 6.Fraction of dissolved drug under sink conditions for lognormal particle-size distributions with mass-median diameters of (a) 2μm and (b) 50μm and 10 values of σ between 0.1 and 1.0.The diffusion-layer thickness h = 10μm.
normalized by the median.The scale parameter σ was in this work calculated from the Span using the equation