Similarity solution of a Stefan drug-release subdiffusion problem

The advection-diffusion transport in random velocity fields often exhibits a non-standard behavior [1] which is referred to as anomalous diffusion. The essence of the phenomenon of the anomalous diffusion is the violation of the relation $\langle | {\bf{x}}(t)-{\bf{x}}(0){| }^{2}\rangle \sim t$ for a random motion x = x(t) of a molecule. Instead, the mean square displacement of the particle scales as $\langle | {\bf{x}}(t)-{\bf{x}}(0){| }^{2}\rangle \sim {t}^{\beta }$, where the exponent β is either larger than 1 (superdiffusion) or smaller than 1 (subdiffusion). Numerous examples of both kinds of anomalous diffusion in flows can be found in review paper [2]. In addition to turbulent convection [3], plasma turbulence [4, 5], and atmospheric turbulence [6], anomalous diffusion has been detected in contaminant transport [7], transport in heterogeneous porous media [8] and gels [9–12].

Anomalous diffusion in gels has potential applications in the problem of controlled drug delivery. Mathematical modeling of controlled drug delivery has been attracting much attention of scientists [13–19]. Analysis of the drug delivery process is often based on traditional models that typically involve a standard Fickian diffusion equation. Due to the scaling properties of this equation, the flux of a drug released from a pill decreases with time as t^−1/2. In applications, it is desirable to release the drug with a nearly constant rate. Anomalous diffusion in gels has been suggested as a mechanism that would provide a slower decrease of the drug flux than in the case of normal diffusion [20–24].

In the present paper we consider drug release from a polymer matrix that undergoes a gel/glass transition. Drug molecules are assumed to perform a subdiffusive motion in the gel and be quiescent in the glass region. We construct an exact self-similar solution for the drug concentration and determine the amount of released drug as a function of time.

Before formulating the problem, we note that subdiffusive transport of the drug in the gel can be described by a fractional diffusion equation for the drug concentration C(x, t)

$$\begin{eqnarray}&&\displaystyle \frac{\partial C(x,t)}{\partial t}=\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}[K{}_{0}{D}_{t}^{1-\beta }C(x,t)].\end{eqnarray} \tag{ 1 }$$

The equation can be derived by means of probabilistic methods [25]. Here

$$\begin{eqnarray}&&{}_{0}{D}_{t}^{1-\beta }C(x,t)=\displaystyle \frac{1}{{\rm{\Gamma }}(\beta )}\displaystyle \frac{\partial }{\partial t}{\displaystyle \int }_{0}^{t}{(t-\tau )}^{\beta -1}C(x,\tau ){\rm{d}}\tau \end{eqnarray} \tag{ 2 }$$

is the Riemann–Liouville fractional derivative and K is a constant fractional diffusion coefficient. One can rewrite equation (1) as

$$\begin{eqnarray}&&{}_{0}^{c}{D}_{t}^{\beta }C(x,t)=\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}[K\ C(x,t)],\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&{}_{0}^{c}{D}_{t}^{\beta }C(x,t)=\displaystyle \frac{1}{{\rm{\Gamma }}(1-\beta )}{\displaystyle \int }_{0}^{t}{(t-\tau )}^{-\beta }\displaystyle \frac{\partial C(x,\tau )}{\partial \tau }{\rm{d}}\tau \end{eqnarray} \tag{ 4 }$$

is the Caputo fractional derivative. Equation (1) can be also written in the form of a conservation law

$$\begin{eqnarray}&&\displaystyle \frac{\partial C(x,t)}{\partial t}+\displaystyle \frac{\partial j(x,t)}{\partial x}=0,\end{eqnarray} \tag{ 5 }$$

where the subdiffusive flux

$$\begin{eqnarray}&&j(x,t)=-K{}_{0}{D}_{t}^{1-\beta }\displaystyle \frac{\partial C(x,t)}{\partial x}.\end{eqnarray} \tag{ 6 }$$

Let us assume now that the polymer in the state of gel, where the subdiffusion of a drug takes place, occupies the region 0 < x < s(t), while in the region $s(t)\lt x\lt \infty$ the polymer is in the state of glass with a vanishing diffusion coefficient. Assume also that the glass/gel transition front moves to the right, and it passes any point x at a certain time instant t_f (x), so that s(t_f (x)) = x. Then the subdiffusive flux in the gel is

$$\begin{eqnarray}&&j(x,t)=-K{}_{{t}_{f}(x)}{D}_{t}^{1-\beta }\displaystyle \frac{\partial C(x,t)}{\partial x},\ 0\lt x\lt s(t),\end{eqnarray} \tag{ 7 }$$

which leads to the subdiffusion equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial C(x,t)}{\partial t}=K\displaystyle \frac{\partial }{\partial x}\left[{}_{{t}_{f}(x)}{D}_{t}^{1-\beta }\displaystyle \frac{\partial C(x,t)}{\partial x}\right],\ 0\lt x\lt s(t).\end{eqnarray} \tag{ 8 }$$

We note the difference between the generic expression for the subdiffusive flux (6) and the expression (7) that we are using. In (6), the history of a drug molecule motion (i.e., the lower limit of integration in the Riemann–Liouville derivative) starts at t = 0. In (7), the history starts at the time t = t_f (x) when the molecule enters the gel region since no motion of the molecule occurs in the glass region. This difference is discussed in greater detail in [26].

At the left boundary of the gel region, x = 0, we assume full consumption of the drug

$$\begin{eqnarray}&&C(0,t)=0.\end{eqnarray} \tag{ 9 }$$

At the front of the glass/gel transition, we expect a jump of the drug concentration due to a restricted solubility of the drug in the solvent:

$$\begin{eqnarray}&&C(s{(t)}^{-},t)={C}_{s},\ C(s{(t)}^{+},t)={C}_{0},\ {C}_{s}\lt {C}_{0}.\end{eqnarray} \tag{ 10 }$$

The conservation of the drug mass leads to the following Stefan boundary condition:

$$\begin{eqnarray}&&({C}_{0}-{C}_{s})\displaystyle \frac{{\rm{d}}s(t)}{{\rm{d}}t}=-j(s{(t)}^{-},t),\ x=s{(t)}^{-}.\end{eqnarray} \tag{ 11 }$$

In the glass region

$$\begin{eqnarray}&&C(x,t)={C}_{0},\ x\gt s(t).\end{eqnarray} \tag{ 12 }$$

The initial conditions are:

$$\begin{eqnarray}&&s(0)=0;\ C(0,0)=0;\ C(x,0)={C}_{0},\ x\gt 0.\end{eqnarray} \tag{ 13 }$$

We remark that while in the case of the normal diffusion, β = 1, the system (8)–(13) has a self-similar solution which describes propagation of the transition front with s(t) ∼ t^1/2 and a finite ∂C/∂x at the front, in the case of subdiffusion, where the flux is determined by the memory integral (7), j(s(t)⁻, t) = 0, unless ∂C/∂x is singular as x → s(t)⁻. In order to determine the limiting value of j(x, t) as x → s(t)⁻, we derive equations that determine the self-similar subdiffusive solution.

Taking into account the scaling properties of the system, we assume

$$\begin{eqnarray}&&s(t)={{at}}^{\beta /2}\end{eqnarray} \tag{ 14 }$$

and introduce the self-similar variable

$$\begin{eqnarray}&&u=t{\left(\displaystyle \frac{x}{a}\right)}^{-2/\beta }.\end{eqnarray} \tag{ 15 }$$

For the self-similar solution C(x, t) = f (u), we find:

$$\begin{eqnarray}&&\displaystyle \frac{\partial C}{\partial x}={f}^{\prime }(u)\displaystyle \frac{\partial u}{\partial x}=-\displaystyle \frac{2}{\beta }{x}^{-1}{{uf}}^{\prime }(u),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial C}{\partial t}={f}^{\prime }(u)\displaystyle \frac{\partial u}{\partial t}={\left(\displaystyle \frac{x}{a}\right)}^{-2/\beta }{f}^{\prime }(u),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&j(x,t)=-K\displaystyle \frac{1}{{\rm{\Gamma }}(\beta )}\displaystyle \frac{\partial }{\partial t}{\displaystyle \int }_{{t}_{f}(x)}^{t}{(t-\tau )}^{\beta -1}\displaystyle \frac{\partial C(x,\tau )}{\partial x}\ {\rm{d}}\tau \\ &&=\displaystyle \frac{2K}{\beta }{\left(\displaystyle \frac{x}{a}\right)}^{1-2/\beta }{a}^{-1}\displaystyle \frac{1}{{\rm{\Gamma }}(\beta )}\displaystyle \frac{{\rm{d}}}{{\rm{d}}u}{\displaystyle \int }_{1}^{u}{(u-v)}^{\beta -1}{{vf}}^{\prime }(v)\mathrm{dv}\\ &&=\displaystyle \frac{2K}{\beta }{\left(\displaystyle \frac{x}{a}\right)}^{1-2/\beta }{a}^{-1}{}_{1}{D}_{u}^{1-\beta }[{{uf}}^{\prime }(u)].\end{eqnarray} \tag{ 18 }$$

Substituting (17), (18) into (5), we obtain the following equation for the self-similar solution f(u):

$$\begin{eqnarray}f^{\prime} (u) & = & \displaystyle \frac{4K}{{\beta }^{2}{a}^{2}}\left\{\left(1-\displaystyle \frac{\beta }{2}\right){}_{1}{D}_{u}^{1-\beta }[{{uf}}^{\prime }(u)]\right.\\ & & +\left.u\displaystyle \frac{{\rm{d}}}{{\rm{d}}u}({}_{1}{D}_{u}^{1-\beta }[{{uf}}^{\prime }(u)])\right\},\quad u\gt 1.\end{eqnarray} \tag{ 19 }$$

It is convenient to rewrite this equation in terms of the function $g(u)={uf}^{\prime} (u)$ as

$$\begin{eqnarray}{Ag}(u) & = & \left(1-\displaystyle \frac{\beta }{2}\right)\ u{}_{1}{D}_{u}^{1-\beta }[g(u)]\\ & & +{u}^{2}\displaystyle \frac{{\rm{d}}}{{\rm{d}}u}({}_{1}{D}_{u}^{1-\beta }[g(u)]),\quad u\gt 1,\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray*}&&A=\displaystyle \frac{{\beta }^{2}{a}^{2}}{4K}.\end{eqnarray*}$$

Boundary conditions (9) and (10) give

$$\begin{eqnarray*}&&f(\infty )=0,\quad f(1)={C}_{s},\end{eqnarray*}$$

respectively, from which

$$\begin{eqnarray}&&{\displaystyle \int }_{1}^{\infty }{f}^{\prime }(v)\ {\rm{d}}v=-{C}_{s}.\end{eqnarray} \tag{ 21 }$$

Finally, the boundary condition (11) yields

$$\begin{eqnarray}&&{a}^{2}=-K{\left(\displaystyle \frac{2}{\beta }\right)}^{2}\displaystyle \frac{1}{({C}_{0}-{C}_{s})}{{}_{1}{D}_{u}^{1-\beta }[{{uf}}^{\prime }(u)]| }_{u=1},\end{eqnarray} \tag{ 22 }$$

which can be written as

$$\begin{eqnarray}&&A=-\displaystyle \frac{1}{({C}_{0}-{C}_{s})}{{}_{1}{D}_{u}^{1-\beta }[g(u)]| }_{u=1}.\end{eqnarray} \tag{ 23 }$$

The solution g(u) of equation (20) can be written in the form

$$\begin{eqnarray*}g(u) & = & {(u-1)}^{1-\beta }{r}_{1}(u)+{(u-1)}^{3-2\beta }{r}_{2}(u)\\ & & +{(u-1)}^{5-3\beta }{r}_{3}(u)+\cdots ,\end{eqnarray*}$$

where r_j(u) are analytic functions that can be expanded in powers of u − 1, with ${r}_{1}(1)=k\ne 0$. All the coefficients of the expansions can be determined in terms of k, while k is arbitrary as it should be since equation (20) is a linear homogeneous equation. We reiterate that for u close to 1

$$\begin{eqnarray*}&&g(u)\sim k{(u-1)}^{1-\beta },\end{eqnarray*}$$

so that

$$\begin{eqnarray*}{\displaystyle \int }_{1}^{u}{(u-v)}^{\beta -1}g(v)\ {\rm{d}}v & \sim & k(u-1){\displaystyle \int }_{0}^{1}{(1-\tau )}^{\beta -1}{\tau }^{1-\beta }\ {\rm{d}}\tau \\ & = & k(u-1)B(2-\beta ,\beta )\\ & = & k(u-1){\rm{\Gamma }}(2-\beta ){\rm{\Gamma }}(\beta )\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}u}{\displaystyle \int }_{1}^{u}{(u-v)}^{\beta -1}g(v)\ {\rm{d}}v{| }_{u=1}=k{\rm{\Gamma }}(2-\beta ){\rm{\Gamma }}(\beta ).\end{eqnarray*}$$

We can rewrite (23) as

$$\begin{eqnarray}&&A=-\displaystyle \frac{k{\rm{\Gamma }}(2-\beta )}{({C}_{0}-{C}_{s})}.\end{eqnarray} \tag{ 24 }$$

Rather than attempting to solve (20) by the series expansion, we solve it by using the Laplace transform. We rewrite (20) for the function y(ξ) = g(u), ξ = u − 1 to obtain

$$\begin{eqnarray}&&{Ay}(\xi ){\rm{\Gamma }}(\beta )=\left(1-\displaystyle \frac{\beta }{2}\right)(\xi +1)\displaystyle \frac{{\rm{d}}}{{\rm{d}}\xi }{\displaystyle \int }_{0}^{\xi }{(\xi -\tau )}^{\beta -1}y(\tau )\ {\rm{d}}\tau \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray*}&&+{(\xi +1)}^{2}\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{\xi }^{2}}{\displaystyle \int }_{0}^{\xi }{(\xi -\tau )}^{\beta -1}y(\tau )\ {\rm{d}}\tau ,\quad \xi \gt 0.\end{eqnarray*}$$

Applying the Laplace transform to (25) and using

$$\begin{eqnarray*}&&{ \mathcal L }\left\{{\displaystyle \int }_{0}^{\xi }{(\xi -\tau )}^{\beta -1}y(\tau )\ {\rm{d}}\tau \right\}=\displaystyle \frac{{\rm{\Gamma }}(\beta )}{{s}^{\beta }}\hat{y}(s),\end{eqnarray*}$$

$$\begin{eqnarray*}&&{ \mathcal L }\left\{\displaystyle \frac{{\rm{d}}}{{\rm{d}}\xi }{\displaystyle \int }_{0}^{\xi }{(\xi -\tau )}^{\beta -1}y(\tau )\ {\rm{d}}\tau \right\}={\rm{\Gamma }}(\beta ){s}^{1-\beta }\hat{y}(s),\end{eqnarray*}$$

$$\begin{eqnarray*}{ \mathcal L }\left\{\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{\xi }^{2}}{\displaystyle \int }_{0}^{\xi }{(\xi -\tau )}^{\beta -1}y(\tau )\ {\rm{d}}\tau \right\} & = & {\rm{\Gamma }}(\beta ){s}^{2-\beta }\hat{y}(s)\\ & & -k{\rm{\Gamma }}(\beta ){\rm{\Gamma }}(2-\beta ),\end{eqnarray*}$$

$$\begin{eqnarray*}&&{ \mathcal L }\left\{(\xi +1)\displaystyle \frac{d}{{\rm{d}}\xi }{\displaystyle \int }_{0}^{\xi }{(\xi -\tau )}^{\beta -1}y(\tau )\ {\rm{d}}\tau \right\}\\ &&\quad =\left(-\displaystyle \frac{{\rm{d}}}{{\rm{d}}s}+1\right){\rm{\Gamma }}(\beta ){s}^{1-\beta }\hat{y}(s),\end{eqnarray*}$$

$$\begin{eqnarray*}&&{ \mathcal L }\left\{{(\xi +1)}^{2}\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{\xi }^{2}}{\displaystyle \int }_{0}^{\xi }{(\xi -\tau )}^{\beta -1}y(\tau )\ {\rm{d}}\tau \right\}\\ &&\quad =\left(\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{s}^{2}}-2\displaystyle \frac{{\rm{d}}}{{\rm{d}}s}+1\right){\rm{\Gamma }}(\beta ){s}^{2-\beta }\hat{y}(s)-k{\rm{\Gamma }}(\beta ){\rm{\Gamma }}(2-\beta ),\end{eqnarray*}$$

we reduce the equation to

$$\begin{eqnarray*}&&s({s}^{1-\beta }\hat{y})^{\prime\prime} +\left(1+\displaystyle \frac{\beta }{2}-2\;s\right){({s}^{1-\beta }\hat{y})}^{\prime }\end{eqnarray*}$$

$$\begin{eqnarray*}&&+\;\left[-\left(1+\displaystyle \frac{\beta }{2}\right)+s-{{As}}^{\beta -1}\right]{s}^{1-\beta }\hat{y}=k{\rm{\Gamma }}(2-\beta ).\end{eqnarray*}$$

Two linearly independent solutions of the homogeneous equation are

$$\begin{eqnarray*}&&{\hat{y}}_{1}(s)={s}^{\beta /2-1}\mathrm{exp}\left(s+\displaystyle \frac{2\sqrt{A}}{\beta }{s}^{\beta /2}\right),\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\hat{y}}_{2}(s)={s}^{\beta /2-1}\mathrm{exp}\left(s-\displaystyle \frac{2\sqrt{A}}{\beta }{s}^{\beta /2}\right),\end{eqnarray*}$$

so that using variation of parameters we find the general solution as

$$\begin{eqnarray*}\hat{y}(s) & = & \displaystyle \frac{k{\rm{\Gamma }}(2-\beta )}{2\sqrt{A}}\left[{\hat{y}}_{1}(s)\displaystyle \int \mathrm{exp}\left(-s-\displaystyle \frac{2\sqrt{A}}{\beta }{s}^{\frac{\beta }{2}}\right)\ {\rm{d}}s\right.\\ & & \left.-{\hat{y}}_{2}(s)\displaystyle \int \mathrm{exp}\left(-s+\displaystyle \frac{2\sqrt{A}}{\beta }{s}^{\frac{\beta }{2}}\right)\ {\rm{d}}s\right]\\ & = & \displaystyle \frac{k{\rm{\Gamma }}(2-\beta )}{2\sqrt{A}}\left[{\hat{y}}_{1}(s)\displaystyle \int \displaystyle \frac{{s}^{\frac{\beta }{2}-1}}{{\hat{y}}_{1}(s)}\ {\rm{d}}s-{\hat{y}}_{2}(s)\displaystyle \int \displaystyle \frac{{s}^{\frac{\beta }{2}-1}}{{\hat{y}}_{2}(s)}\ {\rm{d}}s\right].\end{eqnarray*}$$

We are interested in the solution that decays as $s\to \infty$. This solution can be written in the form

$$\begin{eqnarray}&&\hat{y}(s)=\displaystyle \frac{k{\rm{\Gamma }}(2-\beta )}{2\sqrt{A}}{\hat{y}}_{0}(s),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray*}&&{\hat{y}}_{0}(s)=-{\hat{y}}_{1}(s){\displaystyle \int }_{s}^{\infty }\displaystyle \frac{{\tau }^{\displaystyle \frac{\beta }{2}-1}}{{\hat{y}}_{1}(\tau )}\ {\rm{d}}\tau +{\hat{y}}_{2}(s){\displaystyle \int }_{s}^{\infty }\displaystyle \frac{{\tau }^{\displaystyle \frac{\beta }{2}-1}}{{\hat{y}}_{2}(\tau )}\ {\rm{d}}\tau .\end{eqnarray*}$$

Next, we determine the coefficient k by using the condition (21). We observe that

$$\begin{eqnarray*}\hat{y}(s) & = & {\displaystyle \int }_{0}^{\infty }y(\xi ){{\rm{e}}}^{-s\xi }\ {\rm{d}}\xi ={\displaystyle \int }_{1}^{\infty }y(v-1){{\rm{e}}}^{-s(v-1)}\ {\rm{d}}v\\ & = & {{\rm{e}}}^{s}{\displaystyle \int }_{1}^{\infty }{{vf}}^{\prime }(v){{\rm{e}}}^{-{sv}}\ {\rm{d}}v=-{{\rm{e}}}^{s}\displaystyle \frac{{\rm{d}}}{{\rm{d}}s}{\displaystyle \int }_{1}^{\infty }{f}^{\prime }(v){{\rm{e}}}^{-{sv}}\ {\rm{d}}v.\end{eqnarray*}$$

Thus

$$\begin{eqnarray*}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}s}{\displaystyle \int }_{1}^{\infty }{f}^{\prime }(v){{\rm{e}}}^{-{sv}}\ {\rm{d}}v=-{{\rm{e}}}^{-s}\hat{y}(s)\end{eqnarray*}$$

and integrating with respect to s from s to infinity yields

$$\begin{eqnarray*}&&{\displaystyle \int }_{1}^{\infty }{f}^{\prime }(v){{\rm{e}}}^{-{sv}}\ {\rm{d}}v={\displaystyle \int }_{s}^{\infty }{{\rm{e}}}^{-p}\hat{y}(p)\ {\rm{d}}p.\end{eqnarray*}$$

Evaluating at s = 0, we get

$$\begin{eqnarray*}&&{\displaystyle \int }_{1}^{\infty }{f}^{\prime }(v)\ {\rm{d}}v={\displaystyle \int }_{0}^{\infty }{{\rm{e}}}^{-s}\hat{y}(s)\ {\rm{d}}s.\end{eqnarray*}$$

Using the condition (21) and the expression (26) for $\hat{y}(s)$ in the above equation results in the formula for k

$$\begin{eqnarray*}&&k{\rm{\Gamma }}(2-\beta )=-2{C}_{s}\sqrt{A}{\left[{\displaystyle \int }_{0}^{\infty }{{\rm{e}}}^{-s}{\hat{y}}_{0}(s){\rm{d}}s\right]}^{-1}.\end{eqnarray*}$$

Using this result in (24), we obtain an equation for A

$$\begin{eqnarray}&&\sqrt{A}{\displaystyle \int }_{0}^{\infty }{{\rm{e}}}^{-s}{\hat{y}}_{0}(s)\ {\rm{d}}s=2\displaystyle \frac{{C}_{s}}{{C}_{0}-{C}_{s}}.\end{eqnarray} \tag{ 27 }$$

The equation can be simplified. Indeed

$$\begin{eqnarray*}&&\sqrt{A}{\displaystyle \int }_{0}^{\infty }{{\rm{e}}}^{-s}{\hat{y}}_{0}(s)\ {\rm{d}}s\\ &&\quad =\sqrt{A}{\displaystyle \int }_{0}^{\infty }\left[-{{\rm{e}}}^{-s}{\hat{y}}_{1}(s){\displaystyle \int }_{s}^{\infty }\displaystyle \frac{{\tau }^{\frac{\beta }{2}-1}}{{\hat{y}}_{1}(\tau )}\ {\rm{d}}\tau \right.\\ &&\quad \quad +\left.{{\rm{e}}}^{-s}{\hat{y}}_{2}(s){\displaystyle \int }_{s}^{\infty }\displaystyle \frac{{\tau }^{\frac{\beta }{2}-1}}{{\hat{y}}_{2}(\tau )}\ {\rm{d}}\tau \right]\ {\rm{d}}s\\ &&\quad =\sqrt{A}{\displaystyle \int }_{0}^{\infty }\left[-{s}^{\frac{\beta }{2}-1}\mathrm{exp}\left(\displaystyle \frac{2\sqrt{A}}{\beta }{s}^{\frac{\beta }{2}}\right){\displaystyle \int }_{s}^{\infty }\displaystyle \frac{{\tau }^{\frac{\beta }{2}-1}}{{\hat{y}}_{1}(\tau )}\ {\rm{d}}\tau \right.\\ &&\quad \quad \left.+{s}^{\frac{\beta }{2}-1}\mathrm{exp}\left(\displaystyle \frac{-2\sqrt{A}}{\beta }{s}^{\frac{\beta }{2}}\right){\displaystyle \int }_{s}^{\infty }\displaystyle \frac{{\tau }^{\frac{\beta }{2}-1}}{{\hat{y}}_{2}(\tau )}\ {\rm{d}}\tau \right]\ {\rm{d}}s\\ &&\quad =-{\displaystyle \int }_{0}^{\infty }\left[\displaystyle \frac{{\rm{d}}}{{\rm{d}}s}\left\{\mathrm{exp}\left(\displaystyle \frac{2\sqrt{A}}{\beta }{s}^{\frac{\beta }{2}}\right)\right\}{\displaystyle \int }_{s}^{\infty }\displaystyle \frac{{\tau }^{\frac{\beta }{2}-1}}{{\hat{y}}_{1}(\tau )}\ {\rm{d}}\tau \right.\\ &&\quad \quad \left.+\displaystyle \frac{{\rm{d}}}{{\rm{d}}s}\left\{\mathrm{exp}\left(\displaystyle \frac{-2\sqrt{A}}{\beta }{s}^{\frac{\beta }{2}}\right)\right\}{\displaystyle \int }_{s}^{\infty }\displaystyle \frac{{\tau }^{\frac{\beta }{2}-1}}{{\hat{y}}_{2}(\tau )}\ {\rm{d}}\tau \right]\ {\rm{d}}s\\ &&\quad ={\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{\tau }^{\frac{\beta }{2}-1}}{{\hat{y}}_{1}(\tau )}\ {\rm{d}}\tau +{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{\tau }^{\frac{\beta }{2}-1}}{{\hat{y}}_{2}(\tau )}\ {\rm{d}}\tau -{\displaystyle \int }_{0}^{\infty }[\mathrm{exp}\\ &&\quad \quad \left.\times \left(\displaystyle \frac{2\sqrt{A}}{\beta }{s}^{\frac{\beta }{2}}\right)\displaystyle \frac{{s}^{\frac{\beta }{2}-1}}{{\hat{y}}_{1}(s)}+\mathrm{exp}\left(\displaystyle \frac{-2\sqrt{A}}{\beta }{s}^{\frac{\beta }{2}}\right)\displaystyle \frac{{s}^{\frac{\beta }{2}-1}}{{\hat{y}}_{2}(s)}\right]\ {\rm{d}}s\\ &&\quad =2{\displaystyle \int }_{0}^{\infty }{{\rm{e}}}^{-s}\mathrm{cosh}\left(\displaystyle \frac{2\sqrt{A}}{\beta }{s}^{\frac{\beta }{2}}\right)\ {\rm{d}}s-2.\end{eqnarray*}$$

Thus, equation (27) simplifies to

$$\begin{eqnarray*}&&{\displaystyle \int }_{0}^{\infty }{{\rm{e}}}^{-s}\mathrm{cosh}\left(\displaystyle \frac{2\sqrt{A}}{\beta }{s}^{\displaystyle \frac{\beta }{2}}\right)\ {\rm{d}}s=\displaystyle \frac{{C}_{0}}{{C}_{0}-{C}_{s}}.\end{eqnarray*}$$

Recalling the definition of A, we get

$$\begin{eqnarray}&&F(a,\beta )\equiv {\displaystyle \int }_{0}^{\infty }{{\rm{e}}}^{-s}\mathrm{cosh}\left(\displaystyle \frac{a}{\sqrt{K}}{s}^{\displaystyle \frac{\beta }{2}}\right)\ {\rm{d}}s=\displaystyle \frac{{C}_{0}}{{C}_{0}-{C}_{s}}.\end{eqnarray} \tag{ 28 }$$

Since the right-hand side of (28) is greater than 1 and the left-hand side monotonically increases from 1 to infinity as a goes from zero to infinity, there is a unique solution of the equation for any β, $0\leqslant \beta \leqslant 1$. Note that in the case of normal diffusion, β = 1, the integral on the left-hand side can be computed yielding the well known result

$$\begin{eqnarray*}&&\sqrt{\pi }\displaystyle \frac{a}{2\sqrt{K}}\mathrm{exp}\left(\displaystyle \frac{{a}^{2}}{4K}\right){\rm{erf}}\left(\displaystyle \frac{a}{2\sqrt{K}}\right)=\displaystyle \frac{{C}_{s}}{{C}_{0}-{C}_{s}}.\end{eqnarray*}$$

In the limiting case β → 0, we obtain from (28) that

$$\begin{eqnarray*}&&\displaystyle \frac{a}{\sqrt{K}}=\text{arc cosh}\left(\displaystyle \frac{{C}_{0}}{{C}_{0}-{C}_{s}}\right).\end{eqnarray*}$$

We remark that the function F(a, β) monotonically increases with β, 0 ≤ β ≤ 1, so that a is a decreasing function of β.

The total amount of drug Q(t) released from the polymer by the time t, which is an important characteristic of the process, can be found using (18) as

$$\begin{eqnarray}&&Q(t)=-{\displaystyle \int }_{0}^{t}j(0,\tau )\ {\rm{d}}\tau ={{Mt}}^{\beta /2},\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray*}&&M=\displaystyle \frac{4K}{a{\beta }^{2}{\rm{\Gamma }}(\beta )}\underset{u\to \infty }{\mathrm{lim}}H\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&H=-{u}^{1-\beta /2}\displaystyle \frac{{\rm{d}}}{{\rm{d}}u}{\displaystyle \int }_{1}^{u}{(u-v)}^{\beta -1}{{vf}}^{\prime }(v)\ {\rm{d}}v.\end{eqnarray*}$$

To calculate M, we observe that H can be written as

$$\begin{eqnarray}H & = & -{(\xi +1)}^{1-\beta /2}\displaystyle \frac{{\rm{d}}}{{\rm{d}}\xi }{\displaystyle \int }_{0}^{\xi }{(\xi -\eta )}^{\beta -1}y(\eta )\ {\rm{d}}\eta \\ & & =-{(\xi +1)}^{1-\beta /2}\displaystyle \frac{{\rm{d}}}{{\rm{d}}\xi }\left[{\xi }^{\beta }{\displaystyle \int }_{0}^{1}{(1-\tau )}^{\beta -1}y(\xi \tau )\ {\rm{d}}\tau \right].\end{eqnarray} \tag{ 30 }$$

We have to compute the limit of H as $\xi \to \infty$.To determine the limit, we note that asymptotic behavior of y(η) for large η can be found from the behavior of the Laplace transform $\hat{y}(s)$ for small s. The expression (26) for $\hat{y}(s)$ gives for $s\ll 1$

$$\begin{eqnarray*}&&\hat{y}(s)\sim \displaystyle \frac{k{\rm{\Gamma }}(2-\beta )}{\sqrt{A}}{s}^{\beta /2-1}{I}_{s},\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&{I}_{s}={\displaystyle \int }_{0}^{\infty }{{\rm{e}}}^{-\tau }\mathrm{sinh}\left(\displaystyle \frac{a}{\sqrt{K}}{\tau }^{\displaystyle \frac{\beta }{2}}\right)\ {\rm{d}}\tau ,\end{eqnarray*}$$

so that for η ≫ 1

$$\begin{eqnarray}&&y(\eta )\sim \displaystyle \frac{k{\rm{\Gamma }}(2-\beta )}{\sqrt{A}}{\eta }^{-\beta /2}\displaystyle \frac{{I}_{s}}{{\rm{\Gamma }}\left(1-\displaystyle \frac{\beta }{2}\right)}.\end{eqnarray} \tag{ 31 }$$

Using (31) in (30), we obtain

$$\begin{eqnarray*}&&\underset{u\to \infty }{\mathrm{lim}}H=-\displaystyle \frac{\beta k{\rm{\Gamma }}(\beta ){\rm{\Gamma }}(2-\beta )}{}2{\rm{\Gamma }}\left(1+\displaystyle \frac{\beta }{2}\right)\sqrt{A}{I}_{s}\end{eqnarray*}$$

and accounting for the definition of A and relation (24), we simplify the expression for M yielding

$$\begin{eqnarray}&&Q(t)=({C}_{0}-{C}_{s})\displaystyle \frac{\sqrt{K}{I}_{s}}{{\rm{\Gamma }}\left(1+\displaystyle \frac{\beta }{2}\right)}{t}^{\beta /2}.\end{eqnarray} \tag{ 32 }$$

In the case of normal diffusion, β = 1, the integral I_s in (32) can be evaluated and (32) reduces to the square root dependence

$$\begin{eqnarray*}&&Q(t)=({C}_{0}-{C}_{s})a\mathrm{exp}\left(\displaystyle \frac{{a}^{2}}{4K}\right){t}^{1/2}.\end{eqnarray*}$$

We studied a Stefan problem which describes propagation of a gel/glass transition front in a polymer in the context of drug-release devices. In our model, there is no transport of drug molecules through the glass region, while in the gel region the evolution of drug concentration is governed by a subdiffusion equation. We constructed an exact self-similar solution for the drug concentration and determined the amount of drug released from the glass region with time.

There are several other works in which anomalous diffusion in the gel is studied for drug-release devices [20–24]. The models considered in these papers differ from our model. For example, in [20] the authors obtain a self-similar solution of equation (3), which is equivalent to (1), with the fractional derivative ${}_{0}{D}_{t}^{1-\beta }$ rather than ${}_{{t}_{f}(x)}{D}_{t}^{1-\beta }$. In [21, 22], in addition, the boundary condition (11) is replaced with

$$\begin{eqnarray*}&&({C}_{0}-{C}_{s}){\ }_{0}^{c}{D}_{t}^{\beta }s(t)=K\displaystyle \frac{\partial C(x,t)}{\partial x},\end{eqnarray*}$$

which is not equivalent to (11). In [24], the amount of released drug Q(t) is calculated as proportional to t^1 − β/2 (see (32)) as a result of using the normal diffusion flux through the boundary x = 0 rather than a subdiffusive flux,which should be used in the subdiffusion case.

According to our predictions, the drug flux from a subdiffusive gel will decrease with time as t^β/2 − 1, i.e., faster than in the case of normal diffusion. Therefore, gel subdiffusion cannot be considered as a way to provide a nearly constant drug release, which is desirable in applications. Actually, there are other explanations of the deviation from the Fickian scaling observed in experiments (see, e.g., [27] and the references therein) due to polymer swelling, dissolution, erosion, etc. Analysis of these mechanisms of a non-Fickian drug release is beyond the scope of the present paper.

The work has been partially supported by the US-Israel Binational Science Foundation (grant No. 2008122). VAV gratefully acknowledges the support of NSF grant DMS 1108624.