Accuracy of the Single-mode Model of Controlled Release from Hollow Porous Nanospheres

Brian J. Edwards,^1,*  Aili Wang,²  Carl N. Edwards,¹  and Hengbo Yin²

¹Materials Research and Innovation Laboratory (MRAIL) Department of Chemical and Biomolecular Engineering, University of Tennessee, Knoxville, TN 37996, USA

²Faculty of Chemistry and Chemical Engineering, Jiangsu University, Zhenjiang, 212013, China

* Corresponding Author(E-mail): bje@utk.edu (B.J. Edwards)

---

Table of Content

The kinetics of chemical release from hollow porous nanospheres are examined in terms of a single-mode transport model and the general time-dependent multimode solution to the diffusion equation.

 

 

---

ABSTRACT

The kinetics of diffusion-controlled release from hollow nanoporous spheres of varying thickness were recently investigated using classical mass transport theory. A new model, expressed in terms of a single diffusive mode, was developed to describe the time-dependent mass transport of a loaded chemical species through the nanosphere shell based on physically reasonable boundary and initial conditions, as well as an assumption of monotonicity during the approach to steady-state diffusive behavior. The purpose of this communication is to demonstrate the validity of these assumptions in the long-time limit and to estimate the error associated with the short-time transport. To this end, the kinetics of the single-mode model are compared with a full solution to the time-dependent diffusion equation obtained using numerical techniques, which demonstrates the utility of the single-mode model for all but the initial moments of the transport process. As a result, it is evident that the single-mode model, which expresses the diffusion of the solute from the nanospheres in terms of a single, universal time constant, is a valid approximate model for describing diffusion in these nanoporous systems for controlled release applications.

---

Keywords: Controlled release; Drug delivery; Diffusion; Fick’s law; Hollow nanospheres

---

Introduction

Controlled-released chemical delivery is currently a very active area of research in the fields of medicine, agriculture, chemistry, and materials science.^1-6 The primary purpose of controlled-release technology is to enhance effectiveness by delivering a prescribed dosage of a chemical over a prolonged period of time, rather than all at once.^{7, 8} The principal objectives of controlled-release systems are (1) to maintain a constant and long-term release rate, (2) to achieve a predetermined concentration level over a specified period of time, (3) to deliver a chemical to the target site of application, (4) to assist the chemical in crossing physiochemical barriers, and (5) to protect the chemical from premature elimination by the host. Hence controlled-release technology is currently studied in the agricultural, pharmaceutical, biomedical, and food industries as a potentially beneficial technology for herbicides, pharmaceuticals, fertilizers, and other targeted chemical applications.^9-14

Chemically loaded nanosphere materials are actively being researched as controlled release carriers because of their high specific surface area and nanoporous structure;^15-24 hollow nanosphere materials are particularly interesting because they exhibit a very high chemical loading capacity and allow for efficient transport into and out of their hollow interior cores.^25-30 Modelling diffusion through these hollow nanoporous materials has thus become a critical subject in advancing this new technology so that a complete understanding can be gained of the physiochemical mechanisms that control the chemical transport through the nanosphere walls, based on relevant parameters such as shell thickness and porosity. Although a large body of theoretical methodology for modelling controlled release kinetics in general systems has been developed over the previous decades,^31-37 this technology has not yet been fully applied to hollow nanosphere materials, primarily because of the relative newness of the application.

In a recent publication,³⁸ Wang and Edwards developed a model for the controlled release of the agrichemical glyphosate from hollow silica nanospheres of various diameter and porosity. This model was expressed in terms of a single diffusion mode, which represented the longest length and time scale mode of an infinite series expansion of the full solution to the partial differential equation (PDE) for diffusion from hollow spherical media. The basis of this single-mode solution was that only the longest mode satisfied the requisite physical constraints in the approach to an eventual steady state. In terms of a single diffusive mode, a natural (and exact) solution of the diffusion PDE revealed a universal scaling and a single time constant for the diffusive process, which matched very well the long-time experimental data for release of glyphosate from nanosphere systems possessing a range of diameters and porosities. The issue that remains, which this communication intends to address, is whether this single-mode diffusion model is an accurate enough representation of the short length and time scale diffusion behavior of these systems as compared to a full multimode solution comprised of an infinite series of diffusional modes. To reach a definitive conclusion, a series of numerical calculations were performed on the diffusion PDE to obtain multimodal solutions, as far as known boundary and initial conditions can be ascribed, which were then compared with the predictions of the single-mode model of Ref. 38.

Experiments were conducted on four different silica nanosphere systems, each loaded with 0.237 g of glyphosate via an impregnation method.^38-42 Table 1 provides details of the dimensions and physical parameters of these four nanoporous materials made using standard characterization techniques. Fig. 1 (left panel) shows a typical transmission electron microscopy (TEM) image of one of the nanosphere systems used in the experiments. Deviations from the average values reported in Table 1 were found to be less than 7%.^{38, 42} Fig. 1 (right panel) also displays the released mass fraction of glyphosate from the four samples as a function of time after the samples were dispersed in 1000 ml of distilled water. After a short initial induction period where the release rate of glyphosate was relatively rapid, the released mass fraction decelerated and eventually plateaued as the diffusional driving force between the nanospheres and the solvent vanished.

Table 1 Specific surface area, specific pore volume, average inner shell radius (a), outer shell radius (b), and shell thickness (Δ_r) for the four nanosphere samples.

Sample	Surface area (m2/g)	Pore volume (cm3/g)	a (nm)	b (nm)	Δr (nm)
A	176.8	0.265	161.8	183.8	22.0
B	123.3	0.218	164.3	188.5	24.2
C	84.49	0.189	163.9	190.9	27.0
D	76.11	0.144	170.6	203.4	32.8

The physical system of controlled release from a dispersion of hollow nanospheres can be modelled based on the diffusion from a single particle expressed in terms of a partial differential equation of mass transport through the shell of the sphere. Assuming that diffusion of the solute through the shell is the rate-limiting process (i.e., the release rate from the outer surface of the shell to the surrounding solvent is relatively fast as compared to diffusion through the porous shell, corresponding to a small mass transfer resistance at the outer shell surface), the diffusion equation for this system is unambiguously expressed as⁴³

where C(r, t) is the concentration of solute within the shell, which is a function of the radial coordinate, r ∈ [a, b], and the time, t ∈ [0, ∞), and D is the (constant) pseudo-binary diffusion coefficient for transport of the solute between the inner surface (located at coordinate r = a) and outer surface (r = b) of the nanosphere shell. In terms of a dimensionless concentration,

which assumes values lying within the range of 0 ≤ θ ≤ 1 over r ∈ [a, b], Eq. (1) becomes

where C0 ≡ C(a, 0) is the concentration of solute at the inner shell surface at the initial instant and C∞ ≡ C(b,∞) is the concentration at the outer surface in the infinite time limit; i.e., once the diffusive process has completely decayed away. If no mass transfer barrier  were to exist at the interface of the outer shell and solvent, and assuming an infinite bath of solvent, C∞ ≡ 0 and furthermore C∞ ≡ C(b, t) at all instants in time. However, given the nanoporous nature of the silica spheres, an osmotic pressure differential might exist such that C∞ ≠ 0.

Nevertheless, when expressed in terms of dimensionless concentration, several boundary and initial conditions can be deduced as a bare minimum to begin the solution of Eq. (3). These are

The concentration of solute within the core (hollow center) of the spheroid is assumed to be homogenously distributed, and hence no radial variation exists for r < a. Unfortunately, the conditions of (4) are insufficient to obtain a fully general solution to PDE (3). In order  to derive such a solution, two more conditions are necessary.

The first necessary condition is a mathematical expression for the dimensionless concentration at the inner surface as a function of time; i.e., θ(a, t) = ϕ(t). Another key point is that the total amount of glyphosate within the nanospheres is not concentrated exclusively within the cores at t = 0, but is also contained within the pores of the nanosphere shells. Therefore, one must also have knowledge of the concentration profile throughout the shell at the initial instant, θ(r, 0) = ξ(r). (Note that one could also possibly formulate a problem in terms of Neumann or Robin boundary/initial conditions, but the essential arguments here are unchanged.) With these additional conditions complementing the conditions of (4), a completely general yet unique solution to PDE (3) can be derived, as discussed (for the analogous case of heat conduction from a hollow sphere) by Carslaw and Jaeger.⁴⁴ Unfortunately, definitive mathematical expressions for neither of these two additional conditions can be determined from experimental measurements or physical arguments. 

Given the ambiguity noted above in the functions ϕ(t) and ξ(r), it is unsurprising that a number of solutions to PDE (3) have been proposed over the years, each corresponding to simple choices for these functions. These solutions typically turn out to be infinite series  expansions expressed in terms of diffusive modes, such as 

which corresponds to the solution of PDE (1) when one assumes that C(r, 0) = ξ(r) = C1 and C(a, t) = C(b, t) = ϕ(t) = C2, where C1 and C2 are constants.⁴³ Obviously, this choice of the two functions ϕ(t) and ξ(r) is not physically relevant to the current system, but regardless, this generic form of the solution carries through to other particular solutions of the diffusion PDE based on alternative choices for ϕ(t) and ξ(r), all of which involve an infinite series expansion in terms of individual diffusive modes.

Expressions such as Eq. (5) are only as accurate as the boundary and initial conditions that were used to derive them; however, it is very interesting to observe the consistency between the various particular solutions corresponding to different sets of boundary and initial conditions. These solutions tend invariably to contain an infinite series expansion in terms of both an exponential time decay function and a periodic trigonometric function of the spatial coordinate. The reason for this is clear: PDE (3) is separable in terms of two single variable functions, θ(r, t) = F(t)G(r), such that (3) can be split into two independent contributions,

each of which must be equal to the same undetermined constant, herein set to −λ2D for subsequent convenience. The first equation yields an exponential function and the second results in a trigonometric one. For the exponential function, the time constant of the  series scales with the square of the mode number, (n2), which indicates that the individual modes decay much faster as the mode number increases. This implies that most of the modes die out during the initial moments of a typical experiment, leaving only the longest mode playing a role at the approach to steady-state behavior. As for the trigonometric function, for mode numbers larger than unity (n > 1), the concentration profile, even at steady state, will involve multiple periods in which the spatial variation of concentration can be both positive and negative, and have a slope which can also be both positive and negative. This implies that within certain spatial regions, mass is diffusing contrarily to the applied global concentration gradient, in contradiction of Fick’s Law. Although such an occurrence could be possible based on the imposed initial and boundary conditions in the initial stages of the solution to PDE (1), this sort of behavior should be considered unphysical as the system approaches a steady state in the long time limit.

In response to these issues, Wang and Edwards³⁸ derived a new kinetic model based on PDE (1) involving only a single diffusive mode under the premise that the longest relaxation (diffusion) mode governed the diffusive process over most of the experimental timescale. Under this constraint, a solution to PDE (3) was derived as (see the Supplementary Material of Ref. 38 for the explicit derivation) 

which provides a quantitative profile of the concentration throughout the spherical shell as a function of space and time. In this expression, λ is a purely geometric scaling parameter defined as

The solution (7) was obtained using the boundary and initial conditions of Eq. (4), effectively supplemented by the expressions

When solving PDE (3) according to the boundary and initial conditions of (9), a fully general solution would yield an infinite series composed of only the odd integer modes. In this situation, the first (the longest) exponential in the expansion is nine times (n2 = 9) larger than the second exponential, twenty-five times larger than the third, and so on. Consequently, other than during the initial moments of the process, the first term of the expansion should dominate the overall diffusive process. 

Integration of the concentration profile resulting from Eqs. (2) and (6) from the center of the core to the outer shell surface allowed calculation of the cumulative retained and released fractions of solute as functions of time according to

and

in which MC(t) is the mass remaining in the sphere at any instant in time,  is the initial mass within the sphere at t = 0, and  is the steady-state mass of solute retained within the core/shell; i.e., the long-term release limit (as t → ∞). [Note that in the Supplementary Material of Ref. 38, a typo appeared in Eq. (S.27); however, this typo did not affect any of the subsequent equations or analysis. For completion, a corrected version of that equation is presented in the Appendix of this communication.]

A key observation from Eq. (10) is that a scaled retained mass fraction can be derived as

which is expressed in terms of single characteristic timescale, τ ≡ (λ2D)−1, that dictates the kinetics of the time-release process. Such an expression is a direct result of the fact that the solution to PDE (3) was expressed in terms of only a single diffusive mode and the fact that the solution is separable with respect to space and time, hence implying that the product λ2D is a constant of integration. The utility of this scaled expression is that it removes the effects of the initial and long-term concentrations from the retained mass fraction profile, thus allowing the various nanosphere systems to be compared solely on the basis of their underlying kinetics. Furthermore, Eq. (12) explicitly shows that, in a process governed by only the longest relaxation mode, the product λ2D should be a universal constant independent of the geometry (λ) or porosity (and hence the diffusion coefficient, D) of the nanosphere system. This represents a very important conclusion, which was evidently borne out by the available experimental data.³⁸

The experimental data of the retained mass fraction is plotted in Fig. 2 according to the scaled mass fraction of Eq. (12). As evident in the figure, time-release data of all four systems mostly collapse onto a common curve when expressed in scaled form with a universal time constant of 33.3 min or 2,000 s.³⁸ This universal fit describes the data reasonably well quantitatively over the entire time period, although the fit at short times over predicts the actual data, which are likely multimodal, and causes an associated under prediction at later times in the unimodal exponential profile. The universal timescale can be decomposed into the diffusion coefficient and geometric parameter for each nanosphere system; these values are collected in Table 2. 

Equations (10) and (11) can be used to rescale the data of Fig. 2; a graph of the unscaled retained mass fraction is presented in Fig. 3. In this plot, there remains evidence of either a multimodal solution or a second mass transport mechanism at small time values;  however, in general, the single time constant assumption appears to serve adequately to represent the available data. The remaining questions, therefore, are to what extent a multimodal solution can correct for the short time deviations from unimodal behavior and what is the error associated with neglecting the higher-order diffusive modes over the entire time duration of the experiments? These are the primary questions to be investigated in this communication. 

PDE (3) was solved subject to boundary and initial conditions (4) within the MATLAB numerical computing environment using the pdepe function, which solves initial/boundary value problems involving parabolic/elliptic PDEs. PDE (3) was solved for constants a, b, and D using two 1-dimensional arrays, one for the spatial coordinate (r ∈ [a, b]) and one for the temporal coordinate (t ∈ [0, 7200 s]). In addition to conditions (4), two independent cases of 

Table 2 The product λ²D, the effective diffusive time constant, τ, the geometric parameter, λ, and the diffusion coefficient, D, of the four nanosphere samples.

Sample	λ2D (min-1)	τ (s)	λ (nm-1)	D (cm2/s)
A	0.03	2,000	0.0714	9.80 × 10−16
B	0.03	2,000	0.0649	1.19 × 10−15
C	0.03	2,000	0.0582	1.48 × 10−15
D	0.03	2,000	0.0479	2.18 × 10−15

auxiliary conditions were examined. In the first (Case I), auxiliary conditions (9) were employed, which are consistent with the singlemode diffusion model of Wang and Edwards.³⁸ In the second (Case II), in lieu of any physical rationale or theoretical expectation, the initial concentration profile within the nanosphere shell was arbitrarily set to the constant value C0, implying a uniform concentration profile θ(r, 0) = ξ(r) = 1 across the shell diameter (a ≤ r < b) at t = 0. At the initial instant, at the outer surface θ(b, 0) = 0, with the  boundary condition θ(b, t) = 0 subsequently active at all times. In this way, the pdepe function could determine the most appropriate fit to the boundary and initial conditions compatible with a multimode solution to PDE (3) by integrating inward to the inner surface at r = a. (Explicit details of the pdepe function can be found in the online MATLAB documentation.⁴⁵) Any solution thus obtained was mandated from the separability of PDE (3), as displayed in Eq. (6), to possess an exponential time character and a periodic trigonometric  patial dependence at r = a. Accordingly, time-dependent concentration profiles in terms of either C(r, t) or θ(r, t) were computed over the relevant parameter space [a, b,D]. The retained mass at any instant of time, MC(t), was calculated from the computed concentration profiles using the theoretical expression

Accordingly, the scaled retained mass fraction can be expressed in terms of θ(r, t) as

which can be evaluated numerically using the trapezoidal method.

The dimensionless concentration surface, θ(r, t), is plotted in 3- dimensional format versus the spatial (r) and temporal (t) coordinates for Sample A (a = 161.8 nm, b = 183.8 nm, D = 9.80 × 10−16 cm2/s) in the upper panel of Fig. 4. These data were obtained by solving  PDE (3) using the MATLAB pdepe function, as described above, subject to conditions (4) supplemented with the uniform initial concentration profile θ(r, 0) = ξ(r) = 1, corresponding to Case II. At the initial instant, the concentration profile is constant throughout the shell, as required, dropping precipitously to θ(b,0) = 0 at the outer shell surface; i.e., a Heaviside function. For later times, the profile gradually decays to the uniform steady-state condition θ(r,∞) = 0 at long times.

A more quantitative analysis can be achieved by plotting the same data in 2-dimensional format (the lower panel of Fig. 4) using spatial profiles of θ(r, t) at specific instants of t (Case II: solid lines). Also displayed in the lower panel of Fig. 4 are the predictions of the single-mode model (Eq. (7): dashed lines) at the same instants in time used for Case II. [The Case I solution of the full PDE (3) yields essentially equivalent predictions to those of the single-mode model, demonstrating that the numerical solution to PDE (3) defaults to the single-mode solution under application of the appropriate boundary/initial conditions, as it should.] The general Case II solution decays rapidly from its initial Heaviside function at t = 0 to a monotonically decreasing profile, which assumes a sine function

behavior after approximately 7 minutes have elapsed. In fact, at t ≈ 7 min, the fully general solution matches quite closely the single-mode model prediction at t = 0. Furthermore, for all t > 7 min, the general Case II solution tracks quite closely the single-mode model with a lag time of approximately 7 minutes. As time ultimately becomes long, both solutions essentially decay away to the long time solution θ(r,∞) = 0 with absolute error between the two cases approaching zero with increasing time. Given that the two cases start with vastly different initial concentration profiles across the shell (i.e., one is constant whereas the other is given by the sinusoidal function of Eq. (9)), this behavior is remarkable. What it implies is that the higher-order temporal diffusion modes in the multimode solution rapidly decay, leaving the kinetics controlled primarily by only the longest diffusional mode, which is the only mode present in the single-mode model. Consequently, regardless of the initial concentration profile used as an initial condition on the multimode solution to PDE (3), this solution rapidly decays to the solution provided by the single-mode model. 

The scaled mass fraction as a function of time of the multimode solution is calculated by integrating the concentration profiles of Fig. 4 over the 3-dimensional spatial coordinates according to Eq. (14), whereas this quantity can be determined for the single-mode model  directly from Eq. (12). These two cases are plotted alongside the experimental data for all four samples in Fig. 5. As evident from the figure, for long times (t > 30 min), the two cases essentially coincide, indicating that the approach to steady-state behavior is exactly the same. Furthermore, even for t < 30 min, there is very little relative error, less than 10% at its maximum value. The multimode solution has the higher value since there is more area under the concentration curve of Fig. 4 at t = 0 for the multimode solution than the single-mode model, which implies a larger integral value. Therefore, it is no surprise that for small time values the single-mode model underpredicts the multimode solution. However, the small relative error between the two cases once again vindicates the assumption of  single-mode behavior due to the rapid damping of the higher order temporal modes of the fully general solution. 

Interestingly, the multimode solution does not explain the issue with overprediction of the short time behavior of the experimental data; indeed, it actually enhances the overprediction by shifting the integral upward instead of down. Consequently, this overprediction  cannot be caused by the neglect of higher-order diffusional modes, and therefore it must be due to a secondary diffusional process (i.e., a second mass transport mechanism) that is not captured by PDE (3). Wang and Edwards³⁸ hypothesized that this rapid decay of  the scaled retained mass fraction at small time values was caused by a relatively rapid  (with respect to diffusive transport) release of excess solute initially absorbed on the surface of the nanospheres upon exposure to the solvent. In other words, they conjectured that not all of the glyphosate was absorbed inside the nanospheres during the chemical loading phase, and that this excess material rapidly solvated upon dilution into the solvent. This explanation is still reasonable, based on the conclusions of the present analysis.

Having determined the limitations on the accuracy of the singlemode model with respect to the fully multimode solution, it can now be more confidently employed as a modeling tool for controlled release studies of agrichemicals and pesticides from hollow porous nanosphere systems. As demonstrated above, the primary limitation of the single-mode model is the description of the short time behavior; however, the multiple-mode behavior quickly decays to the longest diffusional mode over the initial stages of a controlled release application. Use of the model was illustrated for several interesting nanosphere systems in prior work.³⁸

Contributions of Authors

BJE performed theoretical/numerical calculations and preparation of the manuscript. AW and HY were responsible for the experimental data/imaging and manuscript editing. CNE developed the numerical methodology, performed numerical calculations, and prepared images for the manuscript. The scaled mass fraction as a function of time of the multimode solution is calculated by integrating the concentration profiles of Fig. 4 over the 3-dimensional spatial coordinates according to Eq. (14), whereas this quantity can be determined for the single-mode model directly from Eq. (12). These two cases are plotted alongside the experimental data for all four samples in Fig. 5. As evident from the figure, for long times (t > 30 min), the two cases essentially coincide, indicating that the approach to steady-state behavior is exactly the same. Furthermore, even for t < 30 min, there is very little relative error, less than 10% at its maximum value. The multimode solution has the higher value since there is more area under the concentration curve of Fig. 4 at t = 0 for the  multimode solution than the single-mode model, which implies a larger integral value. Therefore, it is no surprise that for small time values the single-mode model underpredicts the multimode solution. However, the small relative error between the two cases once again  vindicates the assumption of single-mode behavior due to the rapid damping of the higher order temporal modes of the fully general solution.

Conflicts of interest

None declared.

Acknowledgments

This work was financially supported by research funding from the National Natural Science Foundation of China (201506078) and the Chinese Postdoctoral Science Foundation (2016M601739).

References

1. Y. Zhu, J. Shi, Y. Li, H. Chen, W. Shen and X. Dong, Micropor. Mesopor. Mater., 2005, 85, 75–81.[CrossRef][Scopus]
2. F. Chang, X. Huang, H. Wei, K. Chen, C. Shan and X. Tang, Mater. Lett., 2014, 125, 128–131.[CrossRef][Scopus]
3. B. Chen, G. Quan, Z. Wang, J. Chen, L. Wu, Y. Xu, G. Li and C. Wu, Powder Technol., 2013, 240, 48–53.[CrossRef][Scopus]
4. Z. Xu, P. Ma, C. Li, Z. Hou, X. Zhai, S. Huang and J. Lin, Biomaterials, 2011, 32, 4161–4173.[CrossRef][Scopus]
5. A. Chevillard, H. Angellier-Coussy, V, Guillard, N. Gontard and E. Gastaldi, J. Hazard. Mater., 2012, 205, 32–39.[CrossRef][Scopus]
6. A. Roy, S.K. Singh, J. Bajpai and A.K. Bajpai, Cent. Eur. J. Chem., 2014, 12, 453–469. [CrossRef][Scopus]
7. A. Metaxa, E. Efthimiadou, N. Boukos, E. Fragogeorgi, G. Loudos and G. Kordas, J. Colloid Interf. Sci., 2014, 435, 171–181.[CrossRef][Scopus]
8. L. Yao, Y. Shi, S.R. Jin, M.J. Li and L.M. Zhang, Mater. Res. Bull., 2010, 45, 1351–1356. [CrossRef][Scopus]
9. S. Tashima, S. Shimada, I. Ando, K. Matsumoto, R. Takeda and T. Shiraishi, J. Pesticide Sci., 2000, 25, 128–132.[CrossRef][Scopus]
10. M. Davies, Trends Biotechnol., 2001, 19, 489. [CrossRef][Scopus]
11. E. Kenawy and M.A. Sakran, J. Appl. Polymer Sci., 2001, 80, 415-421.[CrossRef][Scopus]
12. N.A. Peppas and W. Leobandung, J. Biomater. Sci. Polymer Edn., 2004, 15, 125-144.[CrossRef][Scopus]
13. J. Liu, S. Lin, L. Li and E. Liu, Int. J. Pharmaceut., 2005, 298, 117–125. [CrossRef][Scopus]
14. P.P. Nadiminti, J.E. Rookes, Y.D. Dong, C. Sayer, B.J. Boyd and D.M. Cahill, Crop Prot., 2016, 82, 17–29. [CrossRef][Scopus]
15. C.T. Kresge, M.E. Leonowicz, W.J. Roth, J.C. Vartuli and J.S. Beck, Nature, 1992, 359, 710–712. [CrossRef][Scopus]
16. M. Vallet-Regi, A. Rámila, R.P. del Real and J. Pérez-Pariente, Chem. Mater., 2001, 13, 308–311. [CrossRef][Scopus]
17. W. Zhao, J. Gu, L. Zhang, H. Chen and J. Shi, J. Am. Chem. Soc., 2005, 127, 8916–8917. [CrossRef][Scopus]
18. I. Slowing, J. Vivero-Wscoto, C.W. Wu and V.S.Y. Lin, Adv. Drug Delivery Rev., 2008, 60, 1278–1288. [CrossRef][Scopus]
19. C. Lee, L. Lo, C. Mou and C. Yang, Adv. Funct. Mater., 2008, 18, 3283– 3292.[CrossRef][Scopus]
20. M. Manzao, M. Colilla, M. Vallet-Regí and Expert Opin. Drug Delivery, 2009, 6, 1383–1400.[CrossRef][Scopus]
21. P. Yang, Z. Quan, Z. Hou, C. Li, X. Kang, Z. Cheng and J. Lin, Biomaterials, 2009, 30, 4786–4795.[CrossRef][Scopus]
22. J. Lee, N. Lee, H. Kim, J. Kim, S.H. Choi, J.H. Kim, T. Kim, I.C. Song, W.K. Moon and T. Hyeon, J. Am. Chem. Soc., 132 (2) (2010) 552–557.[CrossRef][Scopus]
23. A. Szegedi, M. Popova, I. Goshev and J. Mihály, J. Solid State Chem., 2011, 184, 1201–1207.[CrossRef][Scopus]
24. D. Carmona, F. Balas and J. Santamaria, Mater. Res. Bull., 2014, 59, 311–322.[CrossRef][Scopus]
25. Y. Zhu, J. Shi, W. Shen, H. Chen, X. Dong and M. Ruan, Nanotechnology, 2005, 16, 2633–2638.[CrossRef][Scopus]
26. C. Lin, S. Qiao, C.Z. Yu, S. Ismadji and G. Lu, Micropor. Mesopor. Mater., 2009, 117, 213–219.[CrossRef][Scopus]
27. L. Hao, X. Gong, S. Xuan, H. Zhang, X. Gong, W. Jiang and Z. Chen, Appl. Surf. Sci., 2006, 252, 8724–8733.[CrossRef][Scopus]
28. Y. Zhu, Y. Fang, L. Borchardt and S. Kaskel, Micropor. Mesopor. Mater., 2011, 141, 199–206.[CrossRef][Scopus]
29. Y. Zhu and J. Shi, Micropor. Mesopor. Mater., 2007, 103, 243–249.[CrossRef][Scopus]
30. A. Pei, Z. Shen and G. Yang, Mater. Lett., 2007, 61, 2757–2760.[CrossRef][Scopus]
31. G. Flynn, S. Yalkowsky and T. Roseman, J. Pharm. Sci., 1974, 63, 479-510.[CrossRef][Scopus]
32. J. Siepmann, F. Lecomte and R. Bodmeier, J. Control. Rel., 1999, 60, 379-389.[CrossRef][Scopus]
33. J. Siepmann and F. Siepmann, Internat. J. Pharm., 2008, 364, 328-343.[CrossRef][Scopus]
34. F. Siepmann, B. Leclercq and B. Carlin, J. Siepmann, J. Control. Rel., 2009, 135, 71-79. [CrossRef][Scopus]
35. J. Siepmann and N.A. Peppas, Internat. J. Pharm., 2011, 418, 6-12.[CrossRef][Scopus]
36. J. Siepmann and F. Siepmann, J. Control. Rel., 2012, 161, 351-362.[CrossRef][Scopus]
37. D. Caccavo, S. Cascone, G. Lamberti and A.A. Barba, Mol. Pharmaceutics, 2015, 12, 474-483.[CrossRef][Scopus]
38. A. Wang and B.J. Edwards, Int. J. Heat Mass Trans., 2016, 103, 997-1007.[CrossRef][Scopus]
39. C. Ge, D. Zhang, A. Wang, H. Yin, M. Ren, Y. Liu, T. Jiang and L. Yu, J. Phys. Chem. Solids, 2009, 70, 1432–1437. [CrossRef][Scopus]
40. C. Liu, C. Ge, A. Wang, H. Yin, M. Ren, Y. Zhang, L. Yu, T. Jiang and Korean J. Chem. Eng., 2011, 28, 1458-1463. [CrossRef][Scopus]
41. C. Liu, H. Yin, A. Wang, Z. Wu, G. Wu, T. Jiang, Y. Shen and T. Jiang, Trans. Nonferrous Met. Soc. China, 2012, 22, 1161-1168.[CrossRef][Scopus]
42. C. Ge, D. Zhang, A. Wang, H. Yin, M. Ren, Y. Liu, T. Jiang and L. Yu, J. Phys. Chem. Solids, 2009, 70, 1432–1437.[CrossRef][Scopus]
43. J. Crank, The Mathematics of Diffusion, 2nd ed., Oxford University Press: Oxford, 1975.[CrossRef][Scopus]
44. H.S. Carslaw and J.C. Jaeger, Conduction of Heat in Solids, 2nd ed., Oxford University Press: Oxford, 1959.
45. MATLAB Product Documentation, MathWorks Publishing; Natick, MA (2018):