Dissolution of an ensemble of differently shaped poly-dispersed drug particles undergoing solubility reduction: mathematical modelling doi: http://dx.doi.org/10.5599/admet.841 297 ADMET & DMPK 8(3) (2020) 297-313; doi: http://dx.doi.org/10.5599/admet.841 Open Access : ISSN : 1848-7718 http://www.pub.iapchem.org/ojs/index.php/admet/index Original scientific paper Dissolution of an ensemble of differently shaped poly-dispersed drug particles undergoing solubility reduction: mathematical modelling Michela Abrami 1 , Lucia Grassi 2 , Rosario di Vittorio 1 , Dritan Hasa 3 , Beatrice Perissutti 3 , Dario Voinovich 3 , Gabriele Grassi 4 , Italo Colombo 1 and Mario Grassi 1* 1 Dept. of Engineering and Architecture, Trieste University, Via Alfonso Valerio, 6/A, Trieste, I-34127 Italy 2 Liceo Scientifico G. Galilei, Trieste, Via Mameli 4, I-34139 Italy 3 Dept. of Chemical and Pharmaceutical Sciences, Trieste University, Piazzale Europa 1, Trieste, I-34127, Italy 4 Dept. of Life Sciences, Cattinara University Hospital, Trieste University, Strada di Fiume 447, Trieste, I-34149 Italy *Corresponding Author: E-mail: mario.grassi@dia.units.it; Tel.: +39-040-558-3435; Fax: +39-040-569823 Received: April 30, 2020; Revised: June 18, 2020; Published: July 14, 2020 Abstract The aim of this theoretical paper is to develop a mathematical model for describing the dissolution process, in a finite liquid environment, of an ensemble of poly-dispersed drug particles, in form of sphere, cylinder and parallelepiped that can undergo solubility reduction due to phase transition induced by dissolution. The main result of this work consists in its simplicity as, whatever the particular particles size distribution, only two ordinary differential equations are needed to describe the dissolution process. This, in turn, reflects in a very powerful and agile theoretical tool that can be easily implemented in electronic sheets, a widespread tool among the research community. Another model advantage lies on the possibility of determining its parameters by means of common independent techniques thus enabling the evaluation of the importance of solid wettability on the dissolution process. Keywords Recrystallization; mathematical modelling; dissolution; particles; poly-dispersion; bioavailability Introduction The great variety of controlled drug delivery systems is designed in reason of the drug to be delivered and the clinical target that, in turn, determine the choice of administration routes, among which the most important are the oral, injectable, inhaled, transmucosal, transdermal and the implantable one [1]. Despite the specific delivery system and the administration route considered, the in vivo drug fate is usually characterized by some common steps such as the release from the delivery system and the subsequent tissue absorption and distribution, metabolism and elimination (L-ADME processes) [2]. Specifically to the first step, the release kinetics can be ruled by different phenomena including swelling, erosion, drug dissolution, drug transport (due to diffusion and/or convection), drug interaction with the delivery system, initial drug distribution inside the delivery system and some delivery system geometrical characteristics [3,4]. Independently from the administration route and drug, a key factor for the success and reliability of every delivery system is drug bioavailability, defined as the rate and extent to which the active drug is http://dx.doi.org/10.5599/admet.841 http://dx.doi.org/10.5599/admet.841 http://www.pub.iapchem.org/ojs/index.php/admet/index mailto:mario.grassi@dia.units.it Mario Grassi et al. ADMET & DMPK 8(3) (2020) 297-313 298 absorbed from a pharmaceutical form thus becoming available at the site of drug action [5,6]. In turn, bioavailability depends strongly on drug permeability through cells membrane and drug solubilisation in physiological fluids. This is clearly pointed out by the bio-pharmaceutic classification according to which drugs can be subdivided into four classes: class I (high solubility and high permeability), class II (low solubility, high permeability), class III (high solubility, low permeability) and class IV (low solubility and low permeability) [7]. While permeation implies drug partitioning between a polar aqueous phase and an a- polar phase (cellular membranes), unless active mechanisms rule drug permeation, solubilisation implies the drug dissolution process. This phenomenon must not be confused with the release process as dissolution is defined as “the mixing of two phases with the formation of a new homogeneous phase (i.e. the solution)” [8]. Dissolution, in particular, assumes relevant importance for class II drugs that, interestingly, represent about 40 % of the marketed drugs [9,10] and 70-90 % of the new chemical entities [9–12]. Indeed, most of the drugs are optimized solely on the basis of their pharmacological activity and not for what concerns bioavailability. Examples of commonly marketed drugs that are poorly soluble in water (less than 100 g/cm 3 [13]) include non-steroidal anti-inflammatory drugs (NSAIDs), anticholesterol, antimycotics, antibiotics, anticonvulsants, chemotherapeutics, antivirals, β-blockers, calcium channel blockers and immunosuppressants [4,14-18]. The above considerations make clear why the drug dissolution process is actively studied in the pharmaceutical field from both an experimental and a theoretical viewpoint [4,19]. Experimentalists, typically, prefer to perform Dissolution Rate Test (DRT) [20,21], implying the dissolution of an ensemble of poly-dispersed drug particles in water or in a physiological fluid. In contrast, theorists propend for Intrinsic Dissolution Rate test (IDR), implying the dissolution from the flat surface of a drug tablet fixed to a rotating shaft [4]. While DRT is much closer to the in vivo drug performance, IDR is much simpler to be modelled. Indeed, IDR takes place in a much more controlled frame as fluid hydrodynamics around the flat rotating surface is well understood thanks to the elegant approach of Levich [22]. However, nowadays, the necessity of designing more and more sophisticated delivery systems stimulated theorists to move towards experimentalists. Therefore, interesting studies about the hydrodynamic conditions taking place in DRT [23- 24] and the direct observation of particles size reduction during DRT [19] have been undertaken. Additionally, due to the relevance in the pharmaceutical field, theorists started considering also possible variation of solubility occurring upon dissolution, whose kinetics, as later discussed, essentially depends on the surface available for dissolution, on mass transport resistance occurring at, and around, the solid/liquid interface and on drug solubility in the liquid phase. Indeed, it is well known that the contact between the drug and the liquid environment can lead to polymorphic transformations, as in the case of anhydrous theophylline that becomes hydrated [25] or nicergoline [26], or to amorphous drugs that recrystallize in the most stable crystalline form as it occurs for temazepam [14], nimesulide [27] and posaconazole [28]. Typically, phase transformations imply a reduction of solubility that, unavoidably, reflects in a reduction of drug bioavailability. Thus, these evidences underline the important role played by solubility on dissolution kinetics and increase the complexity of the already complex problem regarding solubility determination [29–31]. Historically, the first fundamental approach aimed at the description of DRT was that of Hixson and Crowell [32-34]. For the first time, they accounted for surface reduction upon dissolution of spherical particles and build up the famous cubic law. Then, the elegant model of Pedersen and co-workers accounted also for spherical particles poly-dispersion [35-38]. Interestingly, this model reduces to the Hixson-Crowell one in the case of monodispersed spherical particles. Since then, many other models were built up. Among others, we can remember the one of Thormann and co-workers that considers also the ADMET & DMPK 8(Y) (2020) 297-313 Mathematical modelling of dissolution rate doi: http://dx.doi.org/10.5599/admet.841 299 possible drug degradation after dissolution [39], the model of Hirai and co-workers that does not account explicitly for the shape of particles but focuses the attention on a law able to describe the time dependence of the dissolution surface [40] and the model of Guo and co-workers focusing on polymorphic transformation induced by dissolution (rifampicin, from form II to I) [41]. The aim of this paper is to build up a mathematical model able to merge a reasonably accurate description of the physical phenomena occurring in DRT with the practical need of experimentalists who demand theoretical tools able to guide experiments design and to allow a reliable data interpretation. In particular, the model will account for drug solubility reduction upon dissolution, particle poly-dispersion, particle geometry (spherical, cylindrical and parallelepiped) and the presence of a finite release environment. Mathematical Modelling Dissolution can be considered as a consecutive process made up by five steps [8,42,43]: 1) contact of the solvent with the solid surface (wetting), which implies the production of a solid/liquid interface starting from solid/vapor one, 2) breakdown of intermolecular bonds in the solid phase (fusion), 3) molecules transfer from the solid phase to the solid/liquid interface (solvation), 4) diffusion of the solvated molecules through the unstirred boundary layer surrounding the solid surface (diffusion), 5) convective transport of solvated drug molecules into the well stirred bulk solution (convection). The first four steps, that can be viewed as the sum of four energy steps, represent the total resistance that the drug molecules have to overcome in order to move from the solid phase to the solution one (dissolution). Obviously, the higher the dissolution energy required (i.e. the higher the mass transfer resistance) the lower the dissolution kinetics is. In order to connect the four steps, we need to know the time evolution of the drug concentration profile in the unstirred boundary layer surrounding the solid surface, whose presence is unavoidable and whose thickness depends on the relative velocity among particles and external fluid, on fluid kinematic viscosity, on particles dimension and on drug diffusion coefficient in the unstirred layer [24]. For this purpose, recourse can be made to the Fick’s second equation:   C D C t      (1) where C is drug concentration, t is time, D is the drug diffusion coefficient in the unstirred layer and  is the gradient vector whose components analytical expression depends on the particular reference system chosen (Cartesian, spherical, cylindrical). The first model assumption relies on the hypothesis that mass transport inside the unstirred layer is one dimensional (the direction is that perpendicular to the solid surface) while the second one consists of the rapid attainment of pseudo-stationary conditions in the unstirred layer (this hypothesis is supported by the numerical solution of Eq. (1) assuming usual values for D (~ 10 –10 m 2 s –1 ) [4] and stagnant layer thickness ≤ 20 m). Thus, Eq. (1) becomes:   0D C   (2) Eq. (2) has to be solved considering the following initial and boundary conditions: Initial: C(ξ) = 0 ξmin < ξ ≤ ξmax (3) Boundary:     min min( )m sD C k C C        n (4) C(ξmax) = Cb (5) http://dx.doi.org/10.5599/admet.841 Mario Grassi et al. ADMET & DMPK 8(3) (2020) 297-313 300 where  is the one dimensional spatial coordinate, (max -min) is the thickness of the unstirred layer (), n is the surface normal versor, Cs is drug solubility in the dissolution liquid while km is the interface mass transfer coefficient mainly depending on the dissolution surface wetting properties and representing dissolution step 1. Eq. (3) affirms that, at the beginning, the unstirred layer does not contain drug molecules while Eq. (4) states that the drug flux leaving the solid surface depends on km and on the difference between drug solubility and drug concentration at the solid-liquid interface on the liquid side. Finally, Eq. (5) imposes that drug concentration equates Cb at the unstirred layer – bulk liquid interface. Eq. (2) solution, in Cartesian, cylindrical and spherical coordinates reads, respectively:  s b b d m ( ) 1 1 C C C C k k              Cartesian (6)   max maxm mb s b d min d min ( ) ln / ln k k C C C C k k                         cylindrical (7)   2 maxm mb s b min max min d d C( ) = 1 / k k C C C k k                          spherical (8) where Cb is the drug concentration in the bulk liquid at time t, Cs is drug solubility in the liquid phase and kd is the mass transport coefficient (= D/) accounting for the fourth dissolution step (drug diffusion through the unstirred liquid layer). Eqs. (6) – (8) make clear that while in the Cartesian case (flat dissolution surface) a linear concentration profile develops in the unstirred layer, not linear concentration profiles take place when the dissolution surface is not flat. In addition, Eqs. (6) – (8) allow to evaluate the drug concentration C0 in  = min, i.e., at the solid-liquid interface (liquid side): m m 0 s b s d d ( ) / (1 ) k k C C C C k k    Cartesian (9)   max maxm m0 b s b d min min d min ln / ln k k C C C C k k                         cylindrical (10)   m m0 b s b min max min d d C = / k k C C C k k            spherical (11) Eqs. (9) – (11) affirm that at t = 0, when Cb = 0, C0 is a fraction of Cs, while C0 = Cs after a very long time when Cb = Cs, i.e. C0 increases with time up to Cs. When no wettability problems occur (km → ∞), C0 is always equal to Cs while it is equal to Cb = 0 when km → 0 (very poorly wettable solids). By means of Eqs. (6) – (8), it is possible writing the ordinary differential equation accounting for drug concentration increase in the liquid phase: V dCb dt =-S (D ∂C ∂ξ )| ξ=ξmax = SK(Cs-Cb) , (12) where V is the liquid volume and K is an overall mass transport coefficient assuming different analytical expressions depending on the particular coordinate system considered: Cartesian cylindrical spherical m d 1 K = 1 / 1 /k k min max max min m d max min / K = ln( / )1 1 ( / 1)k k         in ax max d m min / K = 1 1 ( ( )) m m k k      (13) ADMET & DMPK 8(Y) (2020) 297-313 Mathematical modelling of dissolution rate doi: http://dx.doi.org/10.5599/admet.841 301 Eq. (12) states that the increase of the drug mass in the fluid phase is equal to the drug flow leaving the solid surface available for dissolution. In addition, Eq. (13) says that when no wettability problems arise (km → ∞) K coincides with kd, provided that the thickness of the unstirred layer  is very small (max ≈ min) in the case of cylindrical and spherical geometry (obviously, when km → 0, K → 0). Although K varies with particle dimension as kd depends on particle dimension [24], in order to simplify the scenario, the third hypothesis of the present model considers the K independence on particle dimensions so that an average value has to be considered. As a matter of fact, the great advantage of this usual hypothesis consists in an easy description of the dissolution process characterizing an ensemble of poly-dispersed particles of different shape (parallelepiped, cylinder, sphere). Indeed, it allows assuming that the solid amount dissolved per unit time and surface is the same whatever the surface delimiting the solid particle. Accordingly, in the case of a parallelepiped of dimensions Xi, Yi, Zi, we have: P i i i i i i i i i i i i i d d( ) d( ) d( ) d( ) d d d d d M XY Z X Y Z Y Z X Z XY t t t t t        , (14)    i i i ii i i i s b s b d( ) d( ) d( ) d( ) ρY    -2     ---       -2  =       d d d d X X Y ZK Z Y Z K C C C C t t t t        , (15) where Mi P and  are the parallelepiped mass and density, respectively. As the displacement  of the dissolution front is the same for each one of the parallelepiped surface, from Eq. (15) we have:  i i i0i s b d( ) d( ) d( )d( ) 1 1 1 -2 --- d 2 d 2 d 2 d i X Y Z K X X C C t t t t           (16) In the case of a cylinder of dimensions Ri, Li, we have: c 2 2i i i i i i i i d d( ) d( ) d( ) 2π π d d d d M R L R L LR R t t t t       (17)    i ii i i i s b s b d( ) d( ) 2π  -2π      ---       d d R R K LR RLK C C C C t t        (18)    2 2i ii i s b s b d( ) d( ) π    -2π           ---       2 d d L L K R R K C C C C t t        (19) where Mi c is the cylinder mass. Consequently, it follows:  i ii 0i i 0i s b d( ) d( )d( ) 1 - ; -2 --- = d d 2 d R L K R R L L C C t t t           (20) Finally, for a sphere of radius Ri, we have:     3 s i 2 2i i i i i s b s b 4 d π d d( ) d( )3 4π    -4π  --     - d d d d R M R R K K R R C C C C t t t t                 (21) where Mi s is the sphere mass. Consequently, it follows:  ii 0i s b d( )d( ) - --- = d d R K R R C C t t        . (22) The most important message of Eqs. (14) – (22) relies on the possibility of determining particles dimensions upon dissolution, regardless of the geometry, by the evaluation of just the dissolution front displacement . http://dx.doi.org/10.5599/admet.841 Mario Grassi et al. ADMET & DMPK 8(3) (2020) 297-313 302 In order to complete the model, it is necessary evaluating both Cs and Cb. Indeed, as discussed in the introduction, it is quite common that, upon dissolution, the drug undergoes a phase transformation (polymorphic or amorphous – crystalline) implying a solubility reduction. Typically, this phenomenon is described by a first order reaction [44] occurring at the solid-liquid interface and leading to the following expression for the Cs temporal reduction: Cs = Csf+(Cs-in-Csf)e (-krt) , (23) where Csf and Cs-in are, respectively, the final and initial values of solubility while kr is the recrystallization constant and t is time. As a matter of fact, Eq. (23) accounts for the dissolution step 2 as solubility is directly connected with the crystal network breakdown attitude that is quantified by its melting temperature and enthalpy [45]. In order to account also for a possible recrystallization in the bulk phase (occurring when Cb(t) > Cs(t)), it is necessary considering the following equation: dMc dt = krb V(Cs(t) - Cb(t)) , (24) where Mc is the amount of recrystallized solid and krb is the bulk recrystallization constant that can differ from kr. Obviously, while Eq. (24) works only when Cb(t) exceeds Cs(t), the initial value for Mc is set to zero. The determination of Cb relies on a global mass balance ensuring that the initial solid mass (M0) must be equal, at any time, to the sum of the undissolved mass, the solubilized drug present in the bulk solution and Mc (the very small thickness of the boundary layer renders negligible the drug amount contained in it):     i=N i=N 0 pi i ci=1 0 pi pi b c b i=1 ( ) = ( ) -- M N V M t M N V C t V M t C t V          , (25) where Npi represents the number of particles characterized by a volume Vpi and N is the number of classes into which the particle size distribution is subdivided in. It is interesting to notice that Eq. (25), de facto, accounts for the presence of a finite volume of the liquid phase (V). Indeed, when V → ∞, Cb will be always zero. For its mathematical attitude in describing particle size distributions, the Weibull equation was considered to describe the initial particles size distribution [46]: i min(-(2 ) ) i i 0 =1- V W e V      , (26) where V0 is the total volume of the solid drug, Vi is the volume occupied by all the particles sharing the same characteristic dimension i (radius in the case of spheres and cylinders, X dimension in the case of parallelepipeds), min is the minimum value of i while  and  are two model parameters. Once particle geometry is fixed, Vi enables the determination of Npi by a simple geometric relation. Obviously, when the generic i goes to zero due to dissolution, Npi is set to zero as this particle class has disappeared and it no longer contributes to the dissolution phenomenon. Thus, Eq. (26) serves only to evaluate the initial number of particles belonging to class i th . In conclusion, the proposed model, establishing a connection among the dissolution steps, assumes that the dissolution kinetics can be essentially affected by steps 1, 2 and 4, attributing to step 3 a negligible role. In addition, whatever the number of classes constituting the particles size distribution, this model needs only two differential equations: one among the differential equations appearing in Eq. (16), (20) and (22), aimed at the prediction of the time evolution of the dissolution front position (), plus Eq. (24) accounting for the time evolution of the recrystallized drug amount from the solution. Indeed, once  is known, it is ADMET & DMPK 8(Y) (2020) 297-313 Mathematical modelling of dissolution rate doi: http://dx.doi.org/10.5599/admet.841 303 possible evaluating the dimension i for all the particle class particles resorting to the algebraic relation linking i to  for the different geometries (see the first relation appearing in Eqs. (16), (20) and (22)). Of course, when i ~ 0 due to the dissolution process, the particle class particles is no longer considered letting i = Npi = 0. As this model does not lead to an analytical solution, its iterative numerical solution (relaxation method, relaxation parameter  = 1, relative tolerance = 10 -3 [47]) was performed discretizing the two differential equations by means of the implicit Euler method. In order to ensure the numerical solution accuracy and stability, the time step was set to 0.25 s and the particles size distribution was subdivided into N = 200 classes. Results and Discussion As many parameters affect model output, it would be impossible discussing the effects of all of them and how they interact to affect the dissolution kinetics. Accordingly, we decided to focus on three aspects, namely particle shape (parallelepiped, cylinder, sphere), particle size distribution and the ratio (V + ) between the liquid phase volume (V) and the initial solid drug volume (M0/). Indeed, especially the last two, represent the most important parameters available for the designing of the experimental set up. In order to properly evaluate the effect of particles shape on dissolution kinetics, the parameters of the Weibull distribution, whose differential expression is: i min(-(2 ) ) (α-1)i 0 i min di d( / ) = (2 )(2 ) d V V W e           , (27) were fixed to get similar profile of Eq. (27) for parallelepiped, cylinder and sphere when plotted versus the volume competing to particles identified by the characteristic dimension i. In the case of sphere i = Ri and the corresponding volume is (4/3)Ri 3 , in the case of cylinder i = Ri and the corresponding volume is Ri 3 FR, while, in the case of parallelepiped, i = Xi so that the corresponding volume is Xi 3 FyxFzx. For the sake of simplicity, in the case of cylinders and parallelepipeds, it was assumed that the ratio FR between cylinder radius and length, and the ratios Fyx = Yi/Xi and Fzx = Zi/Xi for parallelepiped were the same for all particles. Accordingly, FR, Fyx and Fzx can be considered as shape factors characterizing, respectively, cylinders and parallelepipeds. Figure 1a shows the trend of Wdi for spheres, parallelepipeds and cylinders vs Vi assuming FR = 2 and Fyx = Fzx = 1 (cubical cylinder and cube). Although the three Wdi do not share exactly the same wideness in Vi, they are characterized by the same peak amplitude and position. In addition, as all of them are very narrow (10 m ≤ Ri ≤ 11 m, 16 m ≤ Xi ≤ 17 m and 9 m ≤ Ri ≤ 10 m for spheres, parallelepipeds and cylinders, respectively) they can be thought to be representative of approximately mono-dispersed size distributions. Relaying on these particle size distributions and assuming typical values for the other model parameters [4], Figure 1b shows model predictions considering different values of the ratio V + between liquid volume (V) and initial particles volume (V0=M0/) in the case of a drug undergoing recrystallization upon dissolution. It is clear that whatever V + , particle shape does not seem important as the dimensionless drug concentration in the liquid phase Cb + (= Cb/Cf; left vertical axis) is very similar for spheres (thick line), parallelepiped-cubes (thin line) and cubical cylinders (dashed line). In addition, when V + is high (≥ 1500), the classical oversaturation peak does not appear (the over saturation condition is never met as Cb + (t) is always lower that than Cs + (t)) and only for smaller V + it occurs, becoming evident for V + ≤ 270. Finally, it is interesting to underline that the Cb + peak always occurs when the time dependent dimensionless drug solubility Cs + (= Cs(t)/Cf; right vertical axis) crosses the Cb + trend. Indeed, from now on, the solution is in over saturation conditions and only after a very long time Cb + equates Cs + . http://dx.doi.org/10.5599/admet.841 Mario Grassi et al. ADMET & DMPK 8(3) (2020) 297-313 304 Figure 1. (a) Weibull differential distribution Wdi vs particles volume Vi for spheres, parallelepipeds and cylinders. Wdi parameters read: sphere –  = 2,  = 0.8 m, Rmin = 10 m -, parallelepiped -  = 2.35,  = 0.9 m, Xmin = 16 m, Fyx = Fzx = 1 (cubes) -, cylinder -  = 2,  = 0.8 m, Rmin = 9 m, FR = 2 (cubical cylinders). (b) Model predictions relative to the particle size distribution shown in Figure 1a. Cb + (left vertical axis) and Cs + (right vertical axis) are, respectively, the dimensionless drug concentration and solubility in the liquid phase (normalized with respect of the drug final solubility Cf) while t + is a dimensionless time defined by t+= (tK) √V0 3⁄ , being K the dissolution constant and V0 the particles volume before dissolution. V + is the ratio between liquid volume (V) and V0. Other model parameters, set according to [4], read:  = 1.5 g/cm 3 , Cs-in = 2*10 -2 g/cm 3 , Csf = 10 -3 g/cm 3 , kr = krb = 10 -2 s -1 , K = 10 -3 cm/s. In order to appreciate the effect of particle shape factor on dissolution kinetics, we considered parallelepipeds characterized by three different values of the shape factor Fzx (0.25, 1 and 10), the same value for Fyx = 1 (parallelepipeds with square basis and different heights) and the same Xi size distribution adopted in Figure 1a (parallelepiped case). The three distributions considered, depicted in Figure 2a, clearly show that all other parameters being equal, the decrease of Fzx implies smaller particles and this, in turn, reflects in an increased dissolution surface. Practically speaking, we are comparing almost monodispersed particles size distributions composed by thin platelets (Fzx = 0.25), cubes (Fzx = 1) and long rods (Fzx = 10). Figure 2. (a) Weibull differential distribution Wdi vs. particles volume Vi for parallelepipeds characterized by three values of the Fzx shape factor (0.25, 1, 10), the same shape factor Fyx = 1 and the parameters adopted in Figure 1a (parallelepiped case):  = 2.35,  = 0.9 m, Xmin = 16 m. (b) Model predictions relative to the particle size distribution shown in Figure 2a. Cb + is the dimensionless drug concentration in the liquid phase (normalized with respect of the drug final solubility Cf) while t + is a dimensionless time defined by t+= (tK) √V0 3⁄ , being K the dissolution constant and V0 the particles volume before dissolution. V + (= 150) is the ratio between liquid volume (V) and V0. Other model parameters, set according to [4], read:  = 1.5 g/cm 3 , Cs-in = 2*10 -2 g/cm 3 , Csf = 10 -3 g/cm 3 , kr = krb = 10 -2 s -1 , K = 10 -3 cm/s. 0 5000 10000 15000 20000 25000 4.0E-09 4.5E-09 5.0E-09 5.5E-09 6.0E-09 W d i (c m -1 ) Vi (cm 3) sphere parallelepiped cylinder a 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 0.2 0.4 0.6 0.8 Cs +Cb + t+ sphere parallelepiped cylinder C+b V + = 150 C+s C+b V + = 270 C+b V + = 750 C+b V + = 1500 b 0 5000 10000 15000 20000 25000 1.0E-09 1.0E-08 1.0E-07 W d i( cm -1 ) Vi(cm 3) Fzx = 1 Fzx = 10 Fzx = 0.25 a 0 2 4 6 8 10 12 0 0.2 0.4 0.6 0.8 Cb + t+ Fzx = 10 Fzx = 0.25 Fzx = 1 b ADMET & DMPK 8(Y) (2020) 297-313 Mathematical modelling of dissolution rate doi: http://dx.doi.org/10.5599/admet.841 305 Assuming V + = 150 and the same model parameters adopted in Figure 1b, model output is shown in Figure 2b. It is clear that, in this case, the effect of particle shape is no longer negligible as Cb + profile is undoubtedly affected by the increasing surface area available for dissolution induced by the different values of the shape factor Fzx. Thus, particle shape factor becomes important as soon as it reflects in a significant variation of the area available for dissolution. In addition, it turns out that, in the case of drugs undergoing recrystallization, the characteristics of the oversaturation peak can also depend on particles shape factor. With the aim of exploring also the effect of particles poly-dispersion, Figure 3a considers parallelepipeds in form of platelets (Fzx = 0.25, Fyx = 1), cubes (Fzx = 1, Fyx = 1) and rods (Fzx = 10, Fyx = 1) characterized by the same poly-dispersion of the Xi dimension. With respect to the distributions shown in Figure 2a, the wideness of the three distributions is now increased as Xi spans from about 1 m up to 100 m while in Figure 2a we have 16 m ≤ Xi ≤ 17 m. Figure 3b witnesses that also in this case, despite poly-dispersion, a clear difference in the Cb + kinetics emerges as, again, platelets are characterized by the smallest particles and rods correspond to the biggest one while cubes put in the middle as shown by the position of each distribution peak. Obviously, the Cb + temporal evolution is a little bit depressed with respect to that reported in Figure 2b as poly-dispersion implies a reduction of the surface area available for dissolution. Again, the importance of poly-dispersion is strictly related to the effects it implies on the surface area available for dissolution. Figure 3. (a) Weibull differential distribution Wdi vs particles volume Vi for parallelepipeds characterized by three values of the Fzx shape factor (0.25, 1, 10), the same shape factor Fyx = 1 and the following parameters:  = 2.35,  = 90 m, Xmin = 1 m. (b) Model predictions relative to the particle size distribution shown in Figure 3a. Cb + is the dimensionless drug concentration in the liquid phase (normalized with respect of the drug final solubility Cf) while t + is a dimensionless time defined by t+= (tK) √V0 3⁄ , being K the dissolution constant and V0 the particles volume before dissolution. V + (= 150) is the ratio between liquid volume (V) and V0. Other model parameters, set according to [4], read:  = 1.5 g/cm 3 , Cs-in = 2*10 -2 g/cm 3 , Csf = 10 -3 g/cm 3 , kr = krb = 10 -2 s -1 , K = 10 -3 cm/s. Case studies After having explored some important model characteristics, with the aim to come closer to experimentalist need, in the following sections the attention will be focused on three different drugs undergoing solubility reduction upon dissolution. The first drug considered is anhydrous theophylline (C7H8N4O2, Mw = 180.2; essentially neutral compound; it is a bronco-dilatator indicated mainly for asthma, bronchospasm, and COPD) that, upon dissolution in water, transforms into the more stable monohydrate form (C7H8N4O2●H2O, Mw = 198.2,  = 1.49 g/cm 3 helium pycnometry) [25]. This polymorphic transformation implies a solubility reduction (T = 25 °C and pH = 7) from Cs-in = 11.6 mg/cm 3 (anhydrous) to Csf = 6.1 mg/cm 3 (monohydrate). In addition, the 0 50 100 150 200 250 1.00E-12 1.00E-10 1.00E-08 1.00E-06 W d i( cm -1 ) Vi(cm 3) a 0 2 4 6 8 10 12 0 0.2 0.4 0.6 0.8 Cb + t+ Fzx = 10 Fzx = 0.25 Fzx = 1 b http://dx.doi.org/10.5599/admet.841 Mario Grassi et al. ADMET & DMPK 8(3) (2020) 297-313 306 analysis of IDR experiments reveals that K = 3.4*10 -3 cm/s and kr = 6*10 -3 s -1 [4]. As one of the key factors ruling dissolution kinetics is the surface area of theophylline particles, particular care has to be devoted to the determination of this parameter. First of all, the inspection of Figure 4a suggests that rods (~ Fyx = Fzx = 0.25) could approximate the shape of anhydrous theophylline particles. Then, on the basis of this consideration, the parameters ruling the Weibull distribution, Eq. (27), are chosen in order to meet the experimentally determined theophylline specific area (0.331 m 2 /g, mercury porosimetry, unpublished data). The following values come out from this procedure:  = 2,  = 115 m and 1 m ≤ Xi ≤ 150 m. Although this is not the unique choice of Eq. (27) parameters leading to the experimentally determined specific surface area, it surely represents a physically sound choice. Relying on this set of parameters, and assuming that the recrystallisation constant in the bulk liquid, krb, equates that on particles surface (kr), model predictions about anhydrous theophylline dissolution are performed assuming different values of V + as this is one of the most important and easy parameter to act on in order to properly design the experimental set up. Figure 4b indicates that the over saturation peak appears for V + ≤ 194 and its amplitude increases with V + reduction up to about 2, i.e. the value of the ratio Cs-in/Csf. The decrease of the dimensionless solubility Cs + , indicated by the dotted line in Figure 4b, always intersects the Cb + trend in its maximum and it is a measure of how fast the polymorphic transition takes place. In addition, Figure 4c allows to evaluate, for different V + values, the temporal evolution of the dimensionless recrystallized drug amount (Mc + ) jointly with the dimensionless drug amount that has not yet undergone dissolution (Ms + ). While Ms + shows a faster decrease for increasing V + , less regular is the behaviour of Mc + due to its sigmoidal shape. The second drug considered is griseofulvin (C17H17ClO6, Mw = 352.8;  = 1.49 g/cm 3 helium pycnometry; neutral compound; it is an antifungal drug used to treat a number of types of dermatophytoses), a typical class II drug that, upon milling, rapidly transforms in its amorphous state. The analysis of IDR tests reveals that, at 37 °C and pH = 7.5, the solubility of amorphous griseofulvin in water is 235 g/cm 3 while crystalline solubility is 60 g/cm 3 . It is worth mentioning that 60 g/cm 3 exceeds the true griseofulvin solubility, i.e. that competing to griseofulvin crystals obtained by double recrystallization of native drug from acetonitrile, that is 12 g/cm 3 . The reason for this discrepancy is due to the presence of many defects in the crystal structure of the commercially available griseofulvin so that it results to be a mixture of macrocrystals, nanocrystals and amorphous phase [4]. As it is well known that solubility increases with the reduction of crystal dimension [45], this discrepancy is not surprising. However, as, typically, DRT tests are performed on commercially available drugs, the solubility value of 60 g/cm 3 will be adopted for the following model simulations. Furtherly, the analysis of IDR test provides the following values for the dissolution K = 3.4*10 -3 cm/s and the recrystallization kr = 9*10 -4 s -1 constants [4]. Again, the determination of the griseofulvin size distribution is performed by looking at the morphology of griseofulvin particles (Figure 5a) and, then, by choosing a reasonable set of Weibull parameters able to lead to the experimentally measured specific surface area (0.954 m 2 /g [48]). Looking at Figure 5a, we can roughly assume that griseofulvin particles are cubes (Fyx = Fzx = 1) while a reasonable choice for the Weibull parameters is:  = 1.49,  = 13.6 m and 0.1 m ≤ Xi ≤ 100 m. Relaying on this set of parameters, and assuming that the recrystallisation constant in the bulk liquid, krb, equates that on particles surface (kr), model predictions about griseofulvin dissolution are performed assuming different values of V + as this is one of the most useful parameters aimed at the proper experimental set up design. Figure 5b underlines that, due to relatively small kr (= krb) value (at least in comparison with that competing to anhydrous theophylline), the over saturation peak is scarcely visible, although present (see in Figure 5a the intersection of the various Cb + trends with Cs + trend), whatever V + . In this case, indeed, we ADMET & DMPK 8(Y) (2020) 297-313 Mathematical modelling of dissolution rate doi: http://dx.doi.org/10.5599/admet.841 307 should more properly speak about an “over saturation plateau” whose entity, obviously, increases with V + decrease. A direct consequence of this behaviour is 1) the fast reduction of the drug amount that has not yet undergone dissolution (Ms + ) at time t, whose reduction kinetics increases with V + (Figure 5c) and 2) the limited amount of the dimensionless recrystallized drug amount (Mc + ) that never exceeds the 10 % of M0 even for the smallest V + considered (Figure 5c). The last drug studied is nimesulide (C13H12N2O5S, Mw = 308.5;  = 1.49 g/cm 3 helium pycnometry; weak acid compound characterized by pKa = 6.5 [49]; it is a non-steroidal anti-inflammatory drug), a typical class II drug that, upon milling, can be transformed in nanocrystals and amorphous state. Figure 4. (a) Picture of anhydrous theophylline particles. (b) Model predictions about theophylline DRT assuming different values of the ratio V + between liquid (V) and particles (V0) volume. Cb + (left vertical axis) and Cs + (right vertical axis) are, respectively, the dimensionless drug concentration and solubility in the liquid phase (normalized with respect to the drug final solubility Cf). (c) Dimensionless amount of drug (Ms + ) still solid after time t and amount of recrystallized drug (Mc + ) in the liquid phase. Both amounts are normalized with respect of the initial drug particles mass M0. Nimesulide is, among the three drugs considered, the less soluble in water being its solubility at 37 °C and pH ≤ 6 equal to 9 g/cm 3 and around 100 g/cm 3 at the same temperature but pH = 7.5 [50]. Fixing the attention on the pH ≤ 6 range, nimesulide solubility increases after co-grinding (in presence of a stabilizer) due to the transformation of macrocrystals into nanocrystals. Although solubility depends on nanocrystals dimension [45], we can set for nanocrystals solubility the value of 28 g/cm 3 that corresponds to four hours co-grinding as described in [27]. Figure 6a inspection reveals that nimesulide crystals look like needles so that their approximation by long cylinder (FR ~ 15) seems the best choice. In order to match the a 300 m 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 100 200 300 400 500 Cs + Cr + t(s) Cr +, V+ = 75 Cr +, V+ = 150 Cr +, V+ = 194 Cr +, V+ = 270 Cs + b 0 10 20 30 40 50 60 70 80 90 100 0 5 10 15 20 25 0 100 200 300 400 500 1 0 0 M s+ 1 0 0 M c + t(s) Ms +, V+ = 75 Ms +, V+ = 150 Ms +, V+ = 270 c http://dx.doi.org/10.5599/admet.841 Mario Grassi et al. ADMET & DMPK 8(3) (2020) 297-313 308 experimentally determined specific area of nimesulide powder after grinding (2.7 m 2 /g, mercury porosimetry; unpublished data), a reasonable choice for the Weibull distribution is:  = 4,  = 0.54 m and 0.3 m ≤ Ri ≤ 0.8 m. Finally, the analysis of the release test performed in [27] let to conclude that the dissolution K is equal to 1.8*10 -5 cm/s and that the recrystallization constant kr is equal to 1.3*10 -2 s -1 . If the small K value witnesses the very low nimesulide propensity to dissolve, this being due to its low wettability (static water contact angle ~ 70° [51]), the high kr value indicates that nanocrystals and the amorphous phase are very instable and tend to rapidly recrystallize in the more stable crystalline form (this is the reason why nimesulide has to be co-ground together with a stabilizer, typically polyvinylpyrrolidone). Relaying on the above parameters set, and assuming that the re-crystallisation constant in the bulk liquid, krb, equates that on particles surface (kr), model predictions about nimesulide dissolution are performed assuming different values of V + . Figure 6b shows the appearance of more and more pronounced oversaturation peak for reducing V + . In addition, the high kr value (in comparison to those competing to theophylline and griseofulvin), makes peak shape very evident as the re-crystallization process is fast. Clearly, when V + is sufficiently high, the peak disappears. The presence of a clearly identifiable peak, however, does not mean that a considerable amount of nimesulide goes in solution, as shown in Figure 6c. Indeed, dimensionless undissolved drug amount (Ms + ) is always close to 100 % M0 and the dimensionless re-crystallized drug amount (Mc + ) is always lower that 0.025% M0, whatever V + considered in the simulations. Figure 5. (a) SEM picture of griseofulvin particles. (b) Model predictions about griseofulvin DRT assuming different values of the ratio V + between liquid (V) and particles (V0) volume. Cb + (left vertical axis) and Cs + (right vertical axis) are, respectively, the dimensionless drug concentration and solubility in the liquid phase (normalized with respect to the drug final solubility Cf). (c) Dimensionless amount of drug (Ms + ) still solid after time t and amount of recrystallized drug (Mc + ) in the liquid phase. Both amounts are normalized with respect of the initial drug particles mass M0. a 2 m 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4 0 100 200 300 400 500 Cs + Cb + t(s) Cb +, V+ = 745 Cr +, V+ = 2980Cs + b 0 2 4 6 8 10 12 0 10 20 30 40 50 60 70 80 90 100 0 100 200 300 400 500 1 0 0 M s+ 1 0 0 M c + t(s) Ms +, V+ = 2980 Ms +, V+ = 745c ADMET & DMPK 8(Y) (2020) 297-313 Mathematical modelling of dissolution rate doi: http://dx.doi.org/10.5599/admet.841 309 Figure 6. (a) SEM picture of nimesulide particles (needles) co-ground with polyvinylpyrrolidone (bigger particles). (b) Model predictions about nimesulide DRT assuming different values of the ratio V + between liquid (V) and particles (V0) volume. Cb + (left vertical axis) and Cs + (right vertical axis) are, respectively, the dimensionless drug concentration and solubility in the liquid phase (normalized with respect to the drug final solubility Cf). (c) Dimensionless amount of drug (Ms + ) still solid after time t and amount of recrystallized drug (Mc + ) in the liquid phase. Both amounts are normalized with respect of the initial drug particles mass M0. Conclusions Assuming that: 1) mass transport inside the unstirred layer surrounding each particle is one dimensional, 2) that a rapid attainment of pseudo-stationary conditions in it takes place and that 3) the dissolution constant K is independent on particle dimensions, the proposed model provides a description of the dissolution kinetics (DRT) from an ensemble of solid drug particles, eventually undergoing re-crystallization in a finite release environment, by means of only two differential equations, regardless of the number of the classes into which the continuous particle size distribution is subdivided. Obviously, this represents a considerable advantage in computational and practical terms as it allows to easily implement the model also by means of electronic sheets that are widespread in the research community. In addition, this model allows to select particle shape among spheres, cylinder (characterized by different length/radius ratios) and parallelepiped (characterized by different Y/X and Z/X ratios), this driving the model close to the real situation. Despite this model requires a considerable number of parameters, the majority of them can be properly determined by means of common independent experiments (typically IDR test, mercury porosimetry, gas adsorption analysis (B.E.T.) and helium pycnometry) and only the dissolution constant K should be a 1 m 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 0 100 200 300 400 500 Cs + Cb + t(s) Cs + b 0 20 40 60 80 100 0 0.005 0.01 0.015 0.02 0.025 0 100 200 300 400 500 1 0 0 M s+ 1 0 0 M c + t(s) Ms +, V+ = 30, 60, 150, 745 c http://dx.doi.org/10.5599/admet.841 Mario Grassi et al. ADMET & DMPK 8(3) (2020) 297-313 310 determined by data fitting. Indeed, as it depends on the relative velocity among particles and the surrounding liquid phase, its theoretical determination would require a complex analysis of the hydrodynamic conditions taking place in the liquid phase. In addition, we have to remember that K depends also on the mass transfer coefficient connected to the wettability of the solid surface (km) (1/K is equal to the sum of the mass transfer resistances due to poor wettability (1/km) and due to the presence of the unstirred layer surrounding each particle (1/kd)). Consequently, the K experimental determination is, in our opinion, the simplest way to proceed also because a proper experimental set up, implying DRT tests performed at increasing values of the Reynolds number (Re) (increasing with the stirring velocity imposed in the liquid phase), allows to separately evaluate kd and km. Indeed, for high Re, K should be almost independent on Re as kd is proportional to Re 0.5 [24] and, thus, dkd/dRe → 0 for high Re. Consequently, K should be mainly dependent on solid wettability (lim𝑅𝑒→∞(1 𝐾⁄ ) = lim 𝑅𝑒→∞ (1 𝑘𝑑⁄ + 1 𝑘𝑚⁄ ) ≈ 1 𝑘𝑚⁄ ). Once km is known, kd could be determined in correspondence of the different Re considered in order to know the function kd(Re) as km is Re independent. Finally, km could be related to the work of solid immersion in the liquid phase. In conclusion, this model allows to study and design DRT experiments permitting to enucleate the rate- determining steps ruling the entire phenomenon. Acknowledgements: none Conflict of interest: The authors declare no conflict of interest. References [1] Handbook of pharmaceutical controlled release technology. D. L. Wise Editor. Marcel Dekker Inc; New York, NY, USA: 2000. [2] N. Bertrand N, J.C. Leroux. The journey of a drug-carrier in the body: an anatomophysiological perspective. Journal of Controlled Release 161 (2012) 152 – 163. [3] M. Grassi, G. Grassi. Application of mathematical modeling in sustained release delivery systems. Expert Opinion on Drug Delivery 11 (2014) 1299 – 1301. [4] M. Grassi, G. Grassi, R. Lapasin, I. Colombo. Understanding Drug Release and Absorption Mechanisms: A Physical and Mathematical Approach, CRC Press, Boca Raton, USA, 2007, 371-479. [5] K. C. Kwan, Oral bioavailability and first pass effects. Drug Metabolism and Disposition 5 (1997) 1329 – 1336. [6] R. Naidu, K.T. Semple, M. Megharaj, A.L. Juhasz, N.S. Bolan, S.K. Gupta, B.E. Clothier and R. Schulin. Bioavailability: definition, assessment and implications for risk assessment. Developments in Soil Science 32 (2008) 39 – 51. [7] G.L. Amidon, H. Lennernäs, V.P. Shah, J.R. Crison. A theoretical basis for a biopharmaceutics drug classification: the correlation of in vitro drug product dissolution and in vivo bioavailability. Pharmaceutical Research 12 (1995) 413 – 420. [8] J. Siepmann, F. Siepmann. Mathematical modelling of drug dissolution. International Journal of Pharmaceutics 453 (2013) 12 – 24. [9] T. Loftsson, M.E. Brewster. Pharmaceutical applications of cyclodextrins: basic science and product development. Journal of Pharmacy and Pharmacology 62 (2010) 1607 – 1621. [10] M. Davis, G. Walker. Recent strategies in spray drying for the enhanced bioavailability of poorly water-soluble drugs. Journal of Controlled Release 269 (2018) 110 – 127. ADMET & DMPK 8(Y) (2020) 297-313 Mathematical modelling of dissolution rate doi: http://dx.doi.org/10.5599/admet.841 311 [11] S. Bertoni, B. Albertini, N. Passerini. Spray Congealing: An Emerging Technology to Prepare Solid Dispersions with Enhanced Oral Bioavailability of Poorly Water Soluble Drugs. Molecules 24 (2019) 3471 1-21. [12] M.R. Gigliobianco, C. Casadidio, R. Censi, P. Di Martino. Nanocrystals of Poorly Soluble Drugs: Drug Bioavailability and Physicochemical Stability. Pharmaceutics 10 (2018) 134 2-29. [13] K.T. Savjani, A.K. Gajjar, J.K. Savjani. Drug Solubility: Importance and Enhancement Techniques. ISRN Pharmaceutics 2012 (2012) 195727 1 – 10. [14] M. Grassi, I. Colombo, R. Lapasin. Drug release from an ensemble of swellable crosslinked polymer particles. Journal of Controlled Release 68 (2000) 97–113. [15] D. Hasa, D. Voinovich, B. Perissutti, M. Grassi, A. Bonifacio, V. Sergo, C. Cepek, M. R. Chierotti, R. Gobetto, S. Dall’Acqua, S. Invernizzi. Enhanced Oral Bioavailability of Vinpocetine through Mechanochemical Salt Formation: Physico-Chemical Characterization and In Vivo Studies. Pharmaceutical Research 28 (2011) 1870-1880. [16] Y. Tsume, D. M. Mudie, P. Langguth, G. E. Amidon, G. L. Amidon. The biopharmaceutics classification system: subclasses for in vivo predictive dissolution (IPD) methodology and IVINC. European Journal of Pharmaceutical Sciences 57 (2014) 152 – 163. [17] S. Baghel, H. Cathcart, N. J. O’Reilly. Polymeric amorphous solid dispersions: a review of amorphization, crystallization, stabilization, solid-state characterization, and aqueous solubilization of biopharmaceutical classification system class II drugs. Journal of Pharmaceutical Sciences 105 (2016) 2527 2544. [18] S. Deshmukh, A. Avachat, A. Garkal, N. Khurana, J. M. Cardot. Optimization of a Dissolution Method in Early Development Based on IVIVC Using Small Animals: Application to a BCS Class II Drug. Dissolution Technologies November (2016) 25 – 41. [19] R. Laitinen, J. Lahtinen, P. Silfsten, E. Vartiainen, P. Jarho, J. Ketolainen. An optical method for continuous monitoring of the dissolution rate of pharmaceutical powders. Journal of Pharmaceutical and Biomedical Analysis 52 (2010) 181 – 189. [20] Z. Li, X He, S. Tian, G. Feng, C. Huang, M. Xun, Z. Wu, Y. Wang. Simultaneous Evaluation of Dissolution and Permeation of Oral Drug Solid Formulations for Predicting Absorption Rate–Limiting Factors and In Vitro–In Vivo Correlations: Case Study Using a Poorly Soluble Weakly Basic Drug. AAPS PharmSciTech 20 (2019) 321 1-13. [21] S. Saboo,N. A. Mugheirbi, D. Y. Zemlyanov, U. S. Kestur, L. S. Taylor. Congruent release of drug and polymer: A “sweetspot” in the dissolution of amorphous solid dispersions. Journal of Controlled Release 298 (2019) 68–82. [22] V. G. Levich. Physicochemical Hydrodynamics, Prentice Hall, Englewood Cliffs, N. J. 1962, pages 60 - 72. [23] S. Kindge, H. Wachtel, B. Abrahamsson, P. Langguth. Computational Fluid Dynamics Simulation of Hydrodynamics and Stresses in the PhEur/USP Disintegration Tester Under Fed and Fasted Fluid Characteristics. Journal of Pharmaceutical Sciences 104 (2015) 2956 – 2968. [24] D.M. D’Arcy, T. Persoons. Mechanistic Modelling and Mechanistic Monitoring: Simulation and Shadowgraph Imaging of Particulate Dissolution in the Flow-Through Apparatus. Journal of Pharmaceutical Science 100 (2011) 1102 – 1115. [25] J.H. de Smidth, J.G. Fokkens, H. Grijseels, D.J.A. Crommelin. Dissolution of theophylline monohydrate and anhydrous theophylline in buffer solutions. Journal of Pharmaceutical Science 75 (1986) 497 – 501. [26] R. Censi, M. R. Gigliobianco, C. Casadidio, P. Di Martino. Changes in the solid state of nicergoline, a poorly soluble drug, under different grinding and environmental conditions: effect on polymorphism and dissolution. Journal of Pharmaceutical Science 108 (2019) 929 – 948. [27] N. Coceani, L. Magarotto, D. Ceschia, I. Colombo M. Grassi. Theoretical and experimental analysis of drug release from an ensemble of polymeric particles containing amorphous and nano-crystalline drug. Chemical Engineering Science 71 (2012) 345 – 355. http://dx.doi.org/10.5599/admet.841 Mario Grassi et al. ADMET & DMPK 8(3) (2020) 297-313 312 [28] A. Elkhabaz, S. Sarkar, G.J. Simpson, L.S. Taylor. Characterization of Phase Transformations for Amorphous Solid Dispersions of a Weakly Basic Drug upon Dissolution in Biorelevant Media. Pharmaceutical Research 36 (2019) 174 1 - 17. [29] A. Avdeef. Solubility Temperature Dependence Predicted from 2D Structure. ADMET & DMPK 3 (2015) 298 – 344. [30] A. Avdeef, E. Fuguet, A. Llinàs, C. Ràfols, E. Bosch, G. Völgyi, T. Verbić, E. Boldyreva, K. Takács-Novák. Equilibrium solubility measurement of ionizable drugs – consensus recommendations for improving data quality. ADMET & DMPK 4 (2016) 117 – 178. [31] C.A.S. Bergström, A. Avdeef. Perspectives in solubility measurement and interpretation. ADMET & DMPK 7 (2019) 88 - 105. [32] A.W. Hixson, J. H. Crowell. Dependence of reaction velocity upon surface and agitation I. Theoretical consideration. Industrial and Engineering Chemistry 23 (1931) 923 – 931. [33] A.W. Hixson, J.H. Crowell. Dependence of reaction velocity upon surface and agitation II. Experimental procedure in study surface. Industrial and Engineering Chemistry 23 (1931) 1002 – 1009. [34] A.W. Hixson, J.H. Crowell. Dependence of reaction velocity upon surface and agitation III. Experimental procedure in study of agitation. Industrial and Engineering Chemistry 23 (1931) 1160 – 1169. [35] V.P. Pedersen, K.F. Brown. Dissolution profile in relation to initial particle distribution. Journal of Pharmaceutical Science 64 (1975) 1192 – 1195. [36] V.P. Pedersen, K.F. Brown. Size distribution effects in multiparticulate dissolution. Journal of Pharmaceutical Science 64 (1977) 1981 – 1986. [37] V.P. Pedersen, K.F. Brown. General class of multiparticulate dissolution models. Journal of Pharmaceutical Science 66 (1977) 1435 – 1438. [38] V.P. Pedersen, J.W. Myrick. Versatile kinetic approach to analysis of dissolution data. Journal of Pharmaceutical Science 67 (1978) 1450 – 1455. [39] U. Thormann, M. De Mieri, M. Neuburger, S. Verjee, P. Altmann, M. Hamburger, G. Imanidis. Mechanism of chemical degradation and determination of solubility by kinetic modeling of the highly unstable sesquiterpene lactone nobilin in different media. Journal of Pharmaceutical Sciences 103 (2014) 3139 – 3152. [40] D. Hirai, Y. Iwao, S. I. Kimura, S. Noguchi, S. Itai. Mathematical model to analyze the dissolution behaviour of metastable crustals or amophous druig accompanied with solid-liquid interface reaction. International Journal of Pharmaceutics 522 (2017) 58 – 65. [41] N. Guo, B. Hou, N. wang, Y. Xiao, J. Hiang, Y. Guo, S. Zong, H. Hao. In situ monitoring and modeling of the solution-mediated polymorphic transformation of rifampicin: from form II to form I. Journal of Pharmaceutical Sciences 107 (2018) 344 – 352. [42] D. Hasa, B. Perissutti, D. Voinovich, M. Abrami, R. Farra, S. M. Fiorentino, G. Grassi, M. Grassi. Drug Nanocrystals: Theoretical Background of Solubility Increase and Dissolution Rate Enhancement. Chemical and Biochemical Engineering Quarterly 28 (2014) 247 – 258. [43] A. Parmar, S. Sharma. Engineering design and mechanistic mathematical models: standpoint on cutting edge drug delivery. Trends in Analytical Chemistry 100 (2018) 15 – 35. [44] H. Nogami, T. Nagai, T. Yotsuyanag. Dissolution phenomena of organic medicinals involving simultaneous phase changes. Chemical Pharmaceutical Bulletin 17 (1969) 499 – 509. [45] G. Chiarappa, A. Piccolo, I. Colombo, D. Hasa, D. Voinovich, M. Moneghini, G. Grassi, R. Farra, M. Abrami, P. Posocco, S. Pricl, M. Grassi. Exploring the Shape Influence on Melting Temperature, Enthalpy, and Solubility of Organic Drug Nanocrystals by a Thermodynamic Model. Crystal Growth Design 17 (2017) 4072 – 4083. [46] B. G. Tenchov, T. K. Yanev. Weibull distribution of particle sizes obtained by uniform random fragmentation. Journal of Colloid Interface Science 111, (1986) 1 – 7. ADMET & DMPK 8(Y) (2020) 297-313 Mathematical modelling of dissolution rate doi: http://dx.doi.org/10.5599/admet.841 313 [47] S.C. Chapra, R.P. Canale. Numerical methods for engineers 3 rd edition, McGraw-Hill, Boston, 1998, pages 233 - 263. [48] K.Y. Chow, D.J.W. Grant. Surface analysis of griseofulvin powders by krypton adsorption: evaluation of specific surface area, BET constant C and polanyi adsorption potential. Powder Technology, 56 (1988) 209 – 223. [49] S. Singh, N. Sharda, L. Mahajan. Spectrophotometric determination of pKa of nimesulide. International Journal of Pharmaceutics 176 (1999) 261–264. [50] M. Grassi, N. Coceani, L. Magarotto. Modelling partitioning of sparingly soluble drugs in a two-phase liquid system. International Journal of Pharmaceutics 239 (2002) 157–169. [51] I. De Simone, N. Coceani, R. Farra, S. M. Fiorentino, G. Grassi, R. Lapasin, D. Hasa, B. Perissutti, M. Grassi, D. Voinovich. Study on polymer-surfactant interactions for the improvement of drug delivery systems wettability. Chemical and Biochemical Engineering Quarterly 26 (2012) 405–415. ©2020 by the authors; licensee IAPC, Zagreb, Croatia. This article is an open-access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/3.0/) http://dx.doi.org/10.5599/admet.841 http://creativecommons.org/licenses/by/3.0/