Drug Release from a Spherical Matrix: Theoretical Analysis for a Finite Dissolution Rate Affected by Geometric Shape of Dispersed Drugs

Amending the neglect of finite dissolution in traditional release models, this study proposed a more generalized drug release model considering the simultaneous dissolution and diffusion procedure from a drug-loaded spherical matrix. How the shape factor (n = 0, 1/2, and 2/3 for the planar, cylindrical, and spherical geometry, respectively) of dispersed drug particles affected the release from the matrix was examined for the first time. Numerical solutions of this generalized model were validated by consensus with a short-time analytical solution for planar drugs and by the approach of the diffusion-controlled limits with Higuchi’s model. The drug release rate increases with the ratio of dissolution/diffusion rate (G) and the ratio of solubility/drug loading (K) but decreases with the shape factor of drug particles. A zero-order release profile is identified for planar drugs before starting the surface depletion layer, and also found for cylindrical and spherical dispersed drugs when K and G are small, i.e. the loaded drug is mainly un-dissolved and the drug release rate is dissolution-controlled. It is also shown that for the case of a small G value, the variation of drug release profile, due to the drug particle geometry, becomes prominent. Detailed comparison with the results of the traditional Higuchi’s model indicates that Higuchi’s model can be applied only when G is large because of the assumption of an instantaneous dissolution. For K = 1/101–1/2, the present analysis suggests an error of 33–85% for drug release predicted by Higuchi’s model for G = 100, 14–44% error for G = 101, while a less than 5% error for G ≧ 103.


Introduction
A reservoir-type drug carrier with un-dissolved drug particles loaded uniformly in a matrix is one of the most common formulations for a controlled release system. However, various mathematical models have been proposed to describe the drug release from this system over the years. Higuchi's model proposed in the early 1960s remain the best known and most popular one [1][2][3]. In his model, a moving front separates the inner region coexisting dissolved and un-dissolved drugs from the outer region of depletion layer, where the un-dissolved drug is exhausted and the concentration of dissolved drug is analyzed by a pseudo-steady state approximation [4,5]. The assumption of instantaneous dissolution in Higuchi's model simplifies the release system to become a purely diffusion-controlled release of dispersed drugs.
After Higuchi's model, many works directed towards various improvements to extend its applications [6][7][8]. In spite of these further refinements, the effect of a finite dissolution rate is still of the surface area A of dispersed drugs and the concentration difference between the drug solubility C s and C. Therefore, S(r,t) = kA(C s − C) (1) where k is the dissolution rate constant with the unit of area −1 time −1 . By the proportion expression, the surface area can be described by the power law of C a and expressed as where A o , C ao , and n are the initial drug surface area per unit volume of sphere, the initial dispersed drug concentration, and the geometry-dependent factor of dispersed drugs, respectively. Assume only the thickness or radius of the drug particles decreases during the release process. It can be derived that the shape factors for planar, cylindrical, and spherical drugs are 0, 1/2, and 2/3, respectively. Applying Fick's second law of diffusion and mass balances for dissolved and dispersed drugs, one obtains the following coupled governing equations of the dissolution-diffusion release model for a spherical matrix: where D is the drug diffusion coefficient and u(r − r * ) is a unit step function. r * in the unit step function indicates the location of the moving boundary of the depletion zone, where the un-dissolved drugs are exhausted.
Before solving the problem, the following non-dimensional variables are introduced: ϕ = C/C t and φ = C a /C t (4) τ = Dt/r o 2 and η = r/r o (5) where r o is the radius of the spherical matrix and C t , defined as the sum of C s and C ao , is the initial drug loading. By substituting these dimensionless parameters into the coupled governing equations for the drug diffusion and dissolution processes, Equation (3) can be transformed into where h * = r * /r o , K = C s /C t , and G = kA o r o 2 /D. K is the ratio of solubility/drug loading and G is a dimensionless number representing the ratio of dissolution rate/diffusion rate. For a diffusion-controlled system, the dissolution rate is much larger than the diffusion rate and therefore G >> 1. It is assumed that saturated dissolved drugs appear originally, and, therefore, initial conditions can be expressed as Applying symmetric and the perfect sink conditions at the center and surface of the drug carrier, respectively, one obtains the following boundary conditions:

Critical Time for the Formation of a Surface Depletion Zone
Let τ c be the critical time for the formation of a surface depletion zone. It is the time at which un-dissolved drugs at the sphere's surface are exhausted. Before the formation of a surface depletion zone, by introducing the boundary condition of ϕ (τ,1) = 0 in Equation (8b) and the initial condition of φ (0,η) = 1 − K in (7b) into the governing equation of Equation (6b), the concentration of un-dissolved drugs at sphere's surface, φ (τ,1), can be obtained, i.e., For τ = τ c , φ(τ,1) = 0, one thus obtains Therefore, τ c = (1 − K)/GK, 2(1 − K)/GK and 3(1 − K)/GK for planar, cylindrical, and spherical drug particles, respectively. For τ < τ c , the drug release process remains a fixed boundary problem. However, for τ > τ c , the boundary of the depletion zone proceeds inward with time and the drug release process turns into a moving boundary problem. Figure 1 shows a schematic diagram of the concentration profile described in the model for τ > τ c . Here, the moving boundary h * , defined by the dimensionless radius of the un-depleted zone, separates a depletion zone containing only dissolved drugs from a region containing both dissolved and un-dissolved drugs.

Critical Time for the Formation of a Surface Depletion Zone
Let  be the critical time for the formation of a surface depletion zone. It is the time at which un-dissolved drugs at the sphere's surface are exhausted. Before the formation of a surface depletion zone, by introducing the boundary condition of ϕ(τ,1) = 0 in Equation (8b) and the initial condition of φ(0,η) = 1-K in (7b) into the governing equation of Equation (6b), the concentration of un-dissolved drugs at sphere's surface, φ (τ,1), can be obtained, i.e., For τ =  , φ(τ,1) = 0, one thus obtains Therefore,  = (1-K)/GK, 2(1-K)/GK and 3(1-K)/GK for planar, cylindrical, and spherical drug particles, respectively. For τ <  , the drug release process remains a fixed boundary problem. However, for τ >  , the boundary of the depletion zone proceeds inward with time and the drug release process turns into a moving boundary problem. Figure 1 shows a schematic diagram of the concentration profile described in the model for τ >  . Here, the moving boundary η * , defined by the dimensionless radius of the un-depleted zone, separates a depletion zone containing only dissolved drugs from a region containing both dissolved and un-dissolved drugs.

Special Case of Planar Drugs
For planar drug particles, the shape factor n = 0 and prior to the formation of the depletion layer, η * = 1 and the unit step function u(η -η * ) in Equation (6) equals zero. Therefore, the drug release problem described in Equation (6) becomes linear, i.e.,

Special Case of Planar Drugs
For planar drug particles, the shape factor n = 0 and prior to the formation of the depletion layer, h * = 1 and the unit step function u(η − η * ) in Equation (6) equals zero. Therefore, the drug release problem described in Equation (6) becomes linear, i.e., The analytic solution of this special case of the boundary value problem was first solved by Harland et al. [17], which gives the expressions of the final solutions of ϕ(τ, η) and φ(τ, η) as

Numerical Solutions
All calculations were performed using MATLAB programs (MathWorks Inc., USA). The PDEPE routine was applied to obtain the numerical solution of Equations (6). In this routine, a spatial discretization technique developed by Skeel and Berzins [30] was applied to transform the partial differential equation into an ordinary differential equation. Time integration using routine ODE15S was then employed to solve the resulting ordinary differential equation. The number of mesh points required for convergence of the numerical solution depended upon the G value. Convergence tests indicated that the larger the G value, the greater the need for more spatial mesh points. To have an acceptable accuracy, a fixed number of spatial mesh points, 2000, was applied for all G values. The relative error tolerance was set to 1 × 10 −3 in PDEPE. The fraction of released drugs was obtained by evaluating the integral of the concentration profile, 1 − ϕ − φ, using the QUAD routine which controlled the relative error tolerance fixed at 1 × 10 −6 . The accuracy of the numerical evaluation of the infinite sums in Equation (12) was controlled at 1 × 10 −4 for all time.

Higuchi's Model
On the basis of an instantaneous dissolution assumption and the presence of a pseudo-steady state drug release at the surface depletion zone, Higuchi's equation gave the following expression for drug release from a spherical matrix immersed in a perfect sink [5]: Here, the concentration of drugs beyond the surface depletion zone was assumed to be constant, and the drugs remaining in the depleted zone could be obtained by integrating the dissolved drug concentration from η * to 1. Therefore,

Validation of the Model
The drug release model in this study was verified for the special case of n = 0. Drug release profiles of numerical and analytical solutions in this work were compared with Higuchi's model for (a) a fixed G with a different K, and (b) a fixed K with a different G in Figure 2. The markers represent analytical solutions describing the release before the formation of a depletion zone, and the solid lines and dash lines, respectively, represent the solution by Higuchi's model and numerical solutions obtained herein. Because the closed form analytical solution can only be obtained before the formation of the depletion zone, it is noticed that analytic solutions are shown only before τ c , i.e., 10 −1 , 10 −2 , and 10 −3 for K = 1/101, 1/11, and 1/2, respectively, in Figure 2a; and 10, 1, and 10 −1 for G = 10, 10 2 , and 10 3 , respectively, in Figure 2b. The exact correspondence of markers to the lines indicates excellent agreement between the analytical and numerical solutions and suggests the accuracy of the numerical solutions. Besides, it is instructive to show the analytical solutions covering the major part of the release process for a small G value because a small G value represents relatively slow dissolution, and the depletion layer will occur in the late release process. The results also show that the effect of K values on applicable ranges of analytical solutions in the release process were slight.
Pharmaceutics 2020, 12, x 6 of 16 numerical solutions. Besides, it is instructive to show the analytical solutions covering the major part of the release process for a small G value because a small G value represents relatively slow dissolution, and the depletion layer will occur in the late release process. The results also show that the effect of K values on applicable ranges of analytical solutions in the release process were slight. The release profiles of different G values in Figure 2b differed from each other. With increasing G values, the drug release rate increased and the release profile approached that of Higuchi's model. In other words, the drug release curves in Figure 2b will be insensitive to the dissolution effect when the G value is large. Higuchi's solution stands for the limiting case of the G value approaching infinity and can be regarded as an asymptotic behavior of G values. Furthermore, numerical solutions in this work showed an initial delay compared with Higuchi's model. The delay was due to the finite The release profiles of different G values in Figure 2b differed from each other. With increasing G values, the drug release rate increased and the release profile approached that of Higuchi's model. In other words, the drug release curves in Figure 2b will be insensitive to the dissolution effect when the G value is large. Higuchi's solution stands for the limiting case of the G value approaching infinity and can be regarded as an asymptotic behavior of G values. Furthermore, numerical solutions in this work showed an initial delay compared with Higuchi's model. The delay was due to the finite dissolution rate, and it was more conspicuous for slowly dissolving drugs in Figure 2b and enhanced by a high initial drug loading in Figure 2a. Figure 2a shows that the deviation between the release profiles in Higuchi's model and this work does not apparently decrease with the decrease in K values. The explanation for this is that the G value of 10 3 used in Figure 2a is still not large enough for Higuchi's model to work. On the other hand, it is somewhat surprising that Higuchi's solution for a K value as large as 1/2 is not very different from the exact numerical solution as shown in Figure 2a. This suggests that Higuchi's formula might also be applicable to drug release systems without a high drug loading ratio, and the accuracy of Higuchi's model is more sensitive to the G value and less sensitive to the K value.

Effects of the Drug Shape Factor
After verifying the exactitude of the numerical solutions, Figure 3 explores the effects of the shape of un-dissolved drug particles, i.e. n values, on the release profile. Calculations were performed for n = 2/3, 1/2, and 0, representing the spherical, cylinder, and planar shaped un-dissolved drug particles dispersed in the matrix, respectively for G = 10 −1 , 10 0 , and 10 1 in the cases of K = 1/2 (Figure 3a-c) and K = 1/101 (Figure 3d-f). As a whole, the release profile displayed a sinusoidal shape of increase with the logarithm of time. The effects of n values became conspicuous in the late stage of the release profile.
Considering the surface area of solid drugs for dissolution in Equation (2), there is a fast decrease in the surface area with decreasing concentration of solid drugs in the release process for a large n. As shown in Figure 3, the case of n = 2/3 had the slowest release rate among the three n values due to the smallest surface area provided for dissolution during the release process. Therefore, the slower the dissolution rate was, the slower the release rate became. The effect of n values was conspicuous for a small G value, because for a large G value the drug release process was diffusion-controlled.
For systems with the same diffusivity, a larger G value indicates a higher dissolution rate and results in a faster exhausting rate as observed in Figure 3a-c for K = 1/2 and Figure 3d-f for K = 1/101. Among these three G values, G = 10 1 (Figure 3c,f) indicated a quick increase in the fraction of released drugs, and the time required to exhaust the entirety of the loaded drugs was shorter than in the case where G = 10 −1 (Figure 3a,d). Furthermore, the effects of K show that with the decrease in the K value, the high un-dissolved drug ratio decreased the drug release rate, and it took a longer time to exhaust drugs within the matrix. Figure 4 further investigates the release rate in Figure 3. Results indicated that a more intense release rate was found for the early release profile, and the release rate decreased with time. For K = 1/2 (Figure 4a-c), half the loaded drugs were dissolved; therefore, the result of a fast release rate for a large G only occurred in late release time. For K = 1/101 (Figure 4d-f), the high ratio of un-dissolved drugs brought about an obvious dissolution effect that the release rate increased with the G value. On the other hand, the difference among three cases of n was small in early release time, but the deviation enlarged in late release time. This phenomenon was enhanced for a small K.
Results discussing the effects of the shape factor of the solid drugs were given by the release curve ( Figure 3) and its corresponding release rate (Figure 4). These figures give data for G values of 10 −1 to 10 1 only. Results for G > 10 1 were not shown because the shape factors of the solid drug particles no longer had a sensitive parameter for drug release in a diffusion-controlled system. From the release rate data, an interesting constant release rate region can be found for planar drug particles, and the effects of the shape factor on the release rate for G values of 10 −1 to 10 5 and K values of 1/2 to 1/101 were further analyzed quantitatively in the following section.
The profile of planar drugs in Figure 4 shows that release profiles in τ < τ c exhibited an initial burst in a short lag-time and a constant release rate thereafter. A zero-order release profile was identified for planar drugs before starting the surface depletion layer. It is sensible to define the τ c as the end point of the constant release rate region. Furthermore, since the release rate is nearly constant within the constant rate region, the starting point of the constant rate region, τ s , is, therefore, defined as the time in a release process first reaching the release rate at τ c with a deviation of less than 1%. The coverage of this zero-order release profile increased for a small K and G. A nearly zero-order release profile was also observed for cylindrical and spherical dispersed drugs when K and G were small. The linear relationship between the release fraction and time resulted from an interesting phenomenon in the profiles of drug concentration before the formation of the depletion layer. Due to the sustained exhausting dissolved drugs, un-dissolved drugs continued dissolving to supply the loss of dissolved drugs and brought about a constantly decreasing concentration of un-dissolved drugs with time [31]. However, the concentration profile of dissolved drugs remained nearly unchanged after the lag time required for the system to reach a quasi-steady state. This feature could cause the embedded drugs to be released at a constant rate [18]. Table 1 summarizes the characteristics of the constant release rate for different K and G values in Figure 4. Results showed that the constant release rate was evident for planar drugs with a high drug loading. The constant drug release rate increased with G, but the duration of constant release rate decreased with G. Although the constant drug release rate decreased with the n value, the constant drug release rate did not vary much for three n cases. The coverage of constant release rate in the whole release profile decreased obviously from n = 0 to 1/2 and 2/3. For K = 101, the coverage of constant release rate in the whole release profile for planar drugs reached 98.28%, 92.60%, and 66.88% for G = 10 −1 , 10 0 , and 10 1 , respectively. Nevertheless, the coverage dropped for cylindrical and spherical drugs.

Evaluation of Possible Deviation by Higuchi Approximation
Higuchi's model is a milestone to estimate the release of a drug-loaded particle, and nowadays is used as a correlation function. However, there is no available data for people to estimate the possible error that might be introduced due to the effect of the finite dissolution rate of the dispersed solid drug particle. In order to provide a reference tool for people to check up the possible error, Figure 5 examines the deviations between Higuchi's model by Equations (13) and (14) and this present model for spherical drugs in different K (1/2, 1/11, and 1/101) and G (10 0 , 10 1 , 10 3 , and 10 5 ) values. The deviation of released drugs was defined by subtracting the release fraction in this present model from Higuchi's model. Positive deviation meant that the release fraction in Higuchi's model was greater than in the present model. Therefore, the positive deviation revealed the overestimation of the release rate by Higuchi's model. were all smaller than 5% for G ≧ 10 3 , regardless of K, Higuchi's model, therefore, afforded an acceptable description of the drug release profile when G ≧ 10 3 .
The practical usage of this proposed model will be limited if one has to perform the calculations by oneself for cases with specific K and G values. As a remedy to this disadvantage, we thus analyzed the deviations between Higuchi's model and this present model and came up with results given in Figure 5 and Table 2. These results cover a wide range of K and G values and provide direct information for the possible maximum error one might have if the release curve is evaluated by the simple Higuchi's model. For a case with a given tolerance on the deviation of drug release, one is thus able to justify the application of the simple Higuchi's model for certain ranges of K and G values. To explain the possible negative deviation, the dynamic radial concentration profiles of dissolved drugs, un-dissolved drugs, and total drugs are demonstrated in Figure 6a-c, respectively, for G = 10 3 and K = 1/2 at τ = 0.003, 0.01, 0.1, and 0.2. The fraction of remaining drugs can be obtained by integrating the sum of the radial concentration profiles of dissolved and un-dissolved drugs. In Figure 6c, it is clearly shown that the total drug concentration evaluated by Higuchi's model is higher than that evaluated by the present model for τ = 0.1 and 0.2, suggesting that the remaining drugs in the matrix predicted by Higuchi's model are higher than the present model and, therefore, a negative deviation occurs at these two time points.

Concentration of un-dissolved drugs
Concentration of total drugs These results can be explained from two aspects: one is the fast development of the surface depletion zone and the possible lower transient surface concentration gradient of Higuchi's model. Firstly, although the initial release rate predicted by Higuchi's model is always faster, the release rate might drop quickly due to the quick development of the surface depletion layer. Since the diffusional resistance is proportional to the square of travel distance, the increase of the surface depletion layer thickness should greatly reduce the drug diffusion rate. As shown in Figure 6b, the concentration profiles of un-dissolved drugs predicted by Higuchi's model give a depletion layer thickness of 0.07 and 0.12 at τ = 0.003 and 0.01, respectively, while a depletion layer thickness of~0 and 0.07 is observed by the present model. This result suggests that the quickly built surface depletion layer might temporarily cause a negative deviation of Higuchi's model. Besides, the negative deviation can also be explained by the concentration profiles of dissolved drugs as shown in Figure 6a. By closely examining the concentration profiles at the surface of the drug carrier, it is observed that the surface concentration gradient of dissolved drugs by the present model is higher than that by Higuchi's model for τ = 0.003 and 0.01. Since the drug release rate is proportional to the surface concentration gradient of dissolved drugs, this result again suggests that Higuchi's model might give a negative deviation during a dynamic release process. Conceptually, the negative deviation of Higuchi's model observed for the rare cases at low drug loading can be explained by the feature of its surface erosion type of drug release. In a surface erosion process, the remaining drugs are mainly allocated in the core region of the carrier. Meanwhile, in the present model, the undissolved drugs are more evenly consumed along radial direction. As a result, the un-dissolved drugs in a carrier with a finite dissolution rate can remain closer to the surface of the carrier matrix and continue to serve as a reservoir to maintain a higher surface concentration gradient of the carrier. For the case of low drug loading and with a saturated initial dissolved drug concentration, the effects of the distribution of remaining drugs in the carrier might be greater than the effect of a finite dissolution rate and, therefore, cause a transient negative deviation.
In Figure 5, this deviation was found to be prominent in the case of a small G value. In the early stage, the deviation was insensitive to the K value, but became prominent at the late stage. In general, the absolute deviation was small for a large G and large K value, i.e., a high dissolution rate and low drug loading such as G = 10 3 and K = 1/2 in this study, which approached the basic assumption in Higuchi's model. Table 2 summarizes the maximum deviations for the release conditions in Figure 5. As shown in the table, Higuchi's model can approximately describe the release curves of G = 10 3 -10 5 under K = 1/101-1/2 within a 5% error. Since the absolute deviations of Higuchi's approximations were all smaller than 5% for G 10 3 , regardless of K, Higuchi's model, therefore, afforded an acceptable description of the drug release profile when G 10 3 . The practical usage of this proposed model will be limited if one has to perform the calculations by oneself for cases with specific K and G values. As a remedy to this disadvantage, we thus analyzed the deviations between Higuchi's model and this present model and came up with results given in Figure 5 and Table 2. These results cover a wide range of K and G values and provide direct information for the possible maximum error one might have if the release curve is evaluated by the simple Higuchi's model. For a case with a given tolerance on the deviation of drug release, one is thus able to justify the application of the simple Higuchi's model for certain ranges of K and G values.

Conclusions
This study derived a more accurate model for drug release that considers a finite dissolution rate and different drug particle shapes from a spherical matrix. The effects of the dispersed drug shape on release profiles were investigated for the first time. Numerical solutions of planar drugs were compared with a short-time analytical solution and with the corresponding diffusion-controlled limits in Higuchi's model to validate this release model developed herein. The results indicated that Higuchi's model is suitable for a large ratio of dissolution/diffusion rate (G) and holds for a ratio of drug solubility/initial loading (K) as large as 1/2. A constant rate of release was found before the formation of the depletion layer for planar drugs, and a large portion of the release process could be described by the analytical solutions in the case of a small G value. As for the effects of drug shape, the rate of release decreased with the shape factor (n) and the variation diminished gradually with the increase in G. The deviation analysis revealed that Higuchi's model could well predict the actual release within a 5% error for the case of G = 10 3 -10 5 under K = 1/101-1/2.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Nomenclature
A drug surface area per volume unit sphere, 1/cm A o initial drug surface area per volume unit sphere, 1/cm C(r,t) dissolved drug concentration, g/mL C a (r,t) un-dissolved drug concentration, g/mL C ao initial un-dissolved drug concentration, g/mL C s solubility of drug in matrix, g/mL C t initial total drug concentration C t = C ao + C s , g/mL D diffusivity of drug in matrix, cm 2 /s G ratio of diffusion rate to dissolution rate kA o r o 2 /D K ratio of drug solubility to initial total drug concentration C s /C t k dissolution rate constant, 1/(cm 2 ·s) n shape factor of a drug particle r radial coordinate from the left of a sphere, cm r o radius of a sphere, cm r * moving front, cm t time, s u unit step function Greek symbols ϕ(η,τ) dimensionless dissolved concentration C/C t φ(η,τ) dimensionless un-dissolved concentration C a /C t η dimensionless radius r/r o η * dimensionless moving front η/r o τ dimensionless time Dt/r 2 o τ c dimensionless critical time τ s dimensionless starting time of constant release rate region