Drug delivery from microcapsules: how can we estimate the release time?

TL;DR

This paper introduces two methods based on temporal moments to estimate drug release time from multi-layer spherical microcapsules, enabling quick comparison of different capsule designs without full simulations.

Contribution

It presents novel moment-based approaches for predicting drug release times from multi-layer microcapsules, avoiding complex transient solution calculations.

Findings

Both methods accurately estimate release times.

Approaches are computationally inexpensive.

Methods enable comparison of capsule configurations.

Abstract

Predicting the release performance of a drug delivery device is an important challenge in pharmaceutics and biomedical science. In this paper, we consider a multi-layer diffusion model of drug release from a composite spherical microcapsule into an external surrounding medium. Based on this model, we present two approaches that provide useful indicators of the release time, i.e. the time required for the drug-filled capsule to be depleted. Both approaches make use of temporal moments of the drug concentration versus time curve at the centre of the capsule, which provide useful insight into the timescale of the process and can be computed exactly without explicit calculation of the full transient solution of the multi-layer diffusion model. The first approach, which uses the zeroth and first temporal moments only, provides simple algebraic expressions involving the various parameters in the model (e.g. layer diffusivities, mass transfer coefficients, partition coefficients) to characterize the release time while the second approach yields an asymptotic estimate of the release time that depends on consecutive higher moments. Through several test cases, we show that both approaches provide a computationally-cheap and useful measure to compare \textit{a priori} the release time of different composite microcapsule configurations.

Tables

Table 1. Table 1: Convergence of the asymptotic estimate of the release time t r ∗ superscript subscript 𝑡 𝑟 ∗ t_{r}^{\ast} ( 4.9 ) for increasing values of the index k 𝑘 k applied to the drug diffusion model ( 2.1 )–( 2.10 ) with a single hydrogel layer ( n = 1 𝑛 1 n=1 ) and the physical parameters ( 5.1 )–( 5.2 ). The value of t r ∗ superscript subscript 𝑡 𝑟 ∗ t_{r}^{\ast} is calculated using three different choices for ε 𝜀 \varepsilon and is displayed in the format days:hours:minutes:seconds.

	$\epsilon ={10}^{-4}$
$k$	$t_r^{\ast }$	${\tilde{c}}_0\left(t_r^{\ast }\right)$
1	02:06:22:59	$7.2152\cdot {10}^{-4}$
2	07:23:03:31	$8.6792\cdot {10}^{-5}$
3	07:21:39:42	$8.8261\cdot {10}^{-5}$
4	07:14:19:34	$9.6411\cdot {10}^{-5}$
5	07:11:49:55	$9.9360\cdot {10}^{-5}$
6	07:11:02:47	$1.0031\cdot {10}^{-4}$
7	07:10:48:25	$1.0060\cdot {10}^{-4}$
8	07:10:44:09	$1.0069\cdot {10}^{-4}$
9	07:10:42:54	$1.0071\cdot {10}^{-4}$
10	07:10:42:33	$1.0072\cdot {10}^{-4}$
11	07:10:42:27	$1.0072\cdot {10}^{-4}$
12	07:10:42:25	$1.0072\cdot {10}^{-4}$
13	07:10:42:25	$1.0072\cdot {10}^{-4}$
14	07:10:42:25	$1.0072\cdot {10}^{-4}$

Table 2. Table 2: Sensitivity of the characteristic timescales, t c ( 1 ) superscript subscript 𝑡 𝑐 1 t_{c}^{(1)} and t c ( 2 ) superscript subscript 𝑡 𝑐 2 t_{c}^{(2)} , and release time, t r ∗ superscript subscript 𝑡 𝑟 ∗ t_{r}^{\ast} with ε = 10 − 4 𝜀 superscript 10 4 \varepsilon=10^{-4} , when varying the outer radius of the capsule ( R 1 subscript 𝑅 1 R_{1} ); the mass transfer coefficient at the external coating ( P 𝑃 P ); the diffusivity in the core ( D 0 subscript 𝐷 0 D_{0} ); and the diffusivity in the hydrogel layer ( D 1 subscript 𝐷 1 D_{1} ). When varying a parameter all other parameters are held fixed at the values given in Eqs ( 5.1 )–( 5.2 ). All results are based on the drug diffusion model ( 2.1 )–( 2.10 ) with a single hydrogel layer ( n = 1 𝑛 1 n=1 ). Values of R 1 subscript 𝑅 1 R_{1} , P 𝑃 P , D 0 subscript 𝐷 0 D_{0} , D 1 subscript 𝐷 1 D_{1} are displayed to three significant figures while values of t c ( 1 ) superscript subscript 𝑡 𝑐 1 t_{c}^{(1)} , t c ( 2 ) superscript subscript 𝑡 𝑐 2 t_{c}^{(2)} and t r ∗ superscript subscript 𝑡 𝑟 ∗ t_{r}^{\ast} are reported to two decimal places (cf. Figure 3 ).

Varying $R_1$
Case	1	2	3	4	5	6
$R_1$ [m]	$1.52\cdot {10}^{-3}$	$1.74\cdot {10}^{-3}$	$1.99\cdot {10}^{-3}$	$2.28\cdot {10}^{-3}$	$2.62\cdot {10}^{-3}$	$3.00\cdot {10}^{-3}$
$t_c^{\left(1\right)}$ [h]	$1.04$	$1.48$	$1.86$	$2.20$	$2.49$	$2.74$
$t_c^{\left(2\right)}$ [h]	$5.32$	$6.03$	$6.80$	$7.61$	$8.45$	$9.32$
$t_r^{\ast }$ [h]	$177.78$	$178.91$	$180.40$	$182.39$	$185.03$	$188.55$
Varying $P$
Case	7	8	9	10	11	12
$P$ [$m{s}^{-1}$]	$1.00\cdot {10}^{-8}$	$2.51\cdot {10}^{-8}$	$6.31\cdot {10}^{-8}$	$1.58\cdot {10}^{-7}$	$3.98\cdot {10}^{-7}$	$1.00\cdot {10}^{-6}$
$t_c^{\left(1\right)}$ [h]	$12.23$	$5.72$	$3.13$	$2.10$	$1.69$	$1.52$
$t_c^{\left(2\right)}$ [h]	$21.86$	$12.73$	$8.77$	$7.08$	$6.39$	$6.11$
$t_r^{\ast }$ [h]	$214.19$	$191.96$	$183.86$	$180.74$	$179.51$	$179.03$
Varying $D_0$
Case	13	14	15	16	17	18
$D_0$ [$m^2{s}^{-1}$]	$3.00\cdot {10}^{-11}$	$7.54\cdot {10}^{-11}$	$1.89\cdot {10}^{-10}$	$4.75\cdot {10}^{-10}$	$1.19\cdot {10}^{-9}$	$3.00\cdot {10}^{-9}$
$t_c^{\left(1\right)}$ [h]	$4.54$	$2.45$	$1.62$	$1.29$	$1.16$	$1.10$
$t_c^{\left(2\right)}$ [h]	$11.35$	$7.54$	$6.22$	$5.74$	$5.55$	$5.48$
$t_r^{\ast }$ [h]	$183.16$	$180.17$	$178.99$	$178.53$	$178.34$	$178.27$
Varying $D_1$
Case	19	20	21	22	23	24
$D_1$ [$m^2{s}^{-1}$]	$5.00\cdot {10}^{-12}$	$1.26\cdot {10}^{-11}$	$3.15\cdot {10}^{-11}$	$7.92\cdot {10}^{-11}$	$1.99\cdot {10}^{-10}$	$5.00\cdot {10}^{-10}$
$t_c^{\left(1\right)}$ [h]	$5.83$	$2.88$	$1.70$	$1.23$	$1.05$	$0.97$
$t_c^{\left(2\right)}$ [h]	$15.36$	$8.70$	$6.41$	$5.63$	$5.34$	$5.23$
$t_r^{\ast }$ [h]	$189.86$	$182.33$	$179.41$	$178.26$	$177.81$	$177.63$

Equations

$$c(x,t)\simeq 0\qquad\text{for {\color[rgb]{0,0,0}{all $x\geq R_{\infty}$ and $t>0$.}}}$$

$$\frac{\partial c_0}{\partial t}=\frac{D_0}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial c_0}{\partial r}\right),$$

$$\frac{\partial c_i}{\partial t}=\frac{D_i}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial c_i}{\partial r}\right),$$

$$\frac{\partial c_e}{\partial t}=\frac{D_e}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial c_e}{\partial r}\right),$$

$$\frac{\partial c_0}{\partial r}=0,$$

$$c_i={\sigma }_i{c}_{i+1},D_i\frac{\partial c_i}{\partial r}=D_{i+1}\frac{\partial c_{i+1}}{\partial r},$$

$$D_n\frac{\partial c_n}{\partial r}=D_e\frac{\partial c_e}{\partial r},D_e\frac{\partial c_e}{\partial r}=P({\sigma }_n{c}_e-c_n),$$

$$c_e=0,$$

$$c_0(r,0)=C_0,$$

$$c_i(r,0)=0,$$

$$c_e(r,0)=0,$$

$$c_i^{\infty }(r)=\underset{t\to \infty }{\lim }{c}_i(r)=0,\text{for}i=0,1,\dots ,n,e.$$

$$t_c^{(1)}={\int }_0^{\infty }{\tilde{c}}_0(t)\mathrm{d}t,$$

$$u_i(r)=\frac{1}{C_0}{\int }_0^{\infty }{c}_i(r,t)\mathrm{d}t,$$

$$t_c^{(1)}=u_0(0).$$

$$\mathcal{L}\phi :=\frac{D_i}{r^2}\frac{\text{d}}{\text{d}r}\left(r^2\frac{\text{d}\phi }{\text{d}r}\right),$$

$$\frac{D_i}{r^2}\frac{\text{d}}{\text{d}r}\left(r^2\frac{\text{d}{u}_i}{\text{d}r}\right)$$

$$\frac{D_0}{r^2}\frac{\text{d}}{\text{d}r}\left(r^2\frac{\text{d}{u}_0}{\text{d}r}\right)=-1,$$

$$\frac{D_i}{r^2}\frac{\text{d}}{\text{d}r}\left(r^2\frac{\text{d}{u}_i}{\text{d}r}\right)=0,$$

$$\frac{D_e}{r^2}\frac{\text{d}}{\text{d}r}\left(r^2\frac{\text{d}{u}_e}{\text{d}r}\right)=0,$$

$$\frac{\text{d}{u}_0}{\text{d}r}=0,$$

$$u_i={\sigma }_i{u}_{i+1},D_i\frac{\text{d}{u}_i}{\text{d}r}=D_{i+1}\frac{\text{d}{u}_{i+1}}{\text{d}r},$$

$$D_n\frac{\text{d}{u}_n}{\text{d}r}=D_e\frac{\text{d}{u}_e}{\text{d}r},D_e\frac{\text{d}{u}_e}{\text{d}r}=P({\sigma }_n{u}_e-u_n),$$

$$u_e=0,$$

$$u_0(r)=\frac{{\alpha }_0}{r}+{\alpha }_1-\frac{r^2}{6D_0},$$

$$u_1(r)=\frac{{\alpha }_2}{r}+{\alpha }_3,$$

$$u_e(r)=\frac{{\alpha }_4}{r}+{\alpha }_5,$$

$${\alpha }_1-\frac{R_0^2}{6D_0}={\sigma }_0\left(\frac{{\alpha }_2}{R_0}+{\alpha }_3\right),$$

$$\frac{R_0}{3}=\frac{D_1{\alpha }_2}{R_0^2},$$

$$\frac{D_1{\alpha }_2}{R_1^2}=\frac{D_e{\alpha }_4}{R_1^2},$$

$$-\frac{D_e{\alpha }_4}{R_1^{2}P}={\sigma }_1\left(\frac{{\alpha }_4}{R_1}+{\alpha }_5\right)-\frac{{\alpha }_2}{R_1}-{\alpha }_3,$$

$$\frac{{\alpha }_4}{R_{\infty }}+{\alpha }_5=0.$$

$$u_0(r)=\frac{R_0^2-r^2}{6D_0}+\frac{{\sigma }_0{R}_0^2(R_1-R_0)}{3D_1{R}_1}+\frac{{\sigma }_0{\sigma }_1{R}_0^3(R_{\infty }-R_1)}{3D_e{R}_1{R}_{\infty }}+\frac{{\sigma }_0{R}_0^3}{3R_1^{2}P},$$

$$u_1(r)=\frac{R_0^3(r-R_1)}{3D_1{R}_{1}r}+\frac{{\sigma }_1{R}_0^3(R_{\infty }-R_1)}{3D_e{R}_1{R}_{\infty }}+\frac{R_0^3}{3R_1^{2}P},$$

$$u_e(r)=\frac{R_0^3(R_{\infty }-r)}{3D_e{R}_{\infty }r},$$

$$t_c^{(1)}=\frac{R_0^2}{6D_0}+\frac{{\sigma }_0{R}_0^2\left(R_1-R_0\right)}{3D_1{R}_1}+\frac{{\sigma }_0{\sigma }_1{R}_0^3\left(R_{\infty }-R_1\right)}{3D_e{R}_1{R}_{\infty }}+\frac{{\sigma }_0{R}_0^3}{3R_1^{2}P}.$$

$$t_c^{(2)}={\int }_0^{\infty }{\tilde{c}}_0(t)\mathrm{d}t+\sqrt{{\int }_0^{\infty }t{\tilde{c}}_0(t)\mathrm{d}t}=t_c^{(1)}+\sqrt{{\int }_0^{\infty }t{\tilde{c}}_0(t)\mathrm{d}t},$$

$$t_c^{(2)}=\frac{R_0^2}{6D_0}+\frac{{\sigma }_0{R}_0^2\left(R_1-R_0\right)}{3D_1{R}_1}+\frac{{\sigma }_0{\sigma }_1{R}_0^3\left(R_{\infty }-R_1\right)}{3D_e{R}_1{R}_{\infty }}+\frac{{\sigma }_0{R}_0^3}{3R_1^{2}P}$$

$$+\{\frac{7R_0^4}{180D_0^2}+\frac{7{\sigma }_0{R}_0^4(R_1-R_0)}{45D_0{D}_1{R}_1}+\frac{7{\sigma }_0{\sigma }_1{R}_0^5(R_{\infty }-R_1)}{45D_0{D}_e{R}_1{R}_{\infty }}+\frac{7{\sigma }_0{R}_0^5}{45D_{0}P{R}_1^2}$$

$$+\frac{2{\sigma }_0{R}_0^3(R_1-R_0{)}^2({\sigma }_0{R}_0+R_1-R_0)}{9D_1^2{R}_1^2}+\frac{4{\sigma }_0{\sigma }_1{R}_0^3(R_{\infty }-R_1)({\sigma }_0{R}_0^3-R_0^3+R_1^3)}{9D_{e}P{R}_1^3{R}_{\infty }}$$

$$+\frac{2{\sigma }_0{R}_0^3(R_1-R_0)(2{\sigma }_0{R}_0^2-2R_0^2+R_0{R}_1+R_1^2)}{9D_1{R}_1^2}\left[\frac{{\sigma }_1(R_{\infty }-R_1)}{D_e{R}_{\infty }}+\frac{1}{P{R}_1}\right]$$

$${+\frac{2{\sigma }_0{\sigma }_1{R}_0^3(R_{\infty }-R_1{)}^2({\sigma }_0{\sigma }_1{R}_0^3-{\sigma }_1{R}_0^3+{\sigma }_1{R}_1^3-R_1^3+R_1^2{R}_{\infty })}{9D_e^2{R}_1^2{R}_{\infty }^2}\}}^{1/2}.$$

$${\tilde{c}}_0(t_r)=\epsilon ,$$

Full text

Keywords: mass diffusion, drug release, release time, composite capsule, asymptotic estimates.

1 Introduction

Polymeric microcapsules are commonly used in pharmaceutical or medical processes as drug carriers or for the encapsulation of organic cells [1, 2]. The main functions of these micro- or nano-sized vesicles are the efficient transport and controlled release of the therapeutic agent into the external environment. Typically, drug-filled microcapsules exhibit a multi-layered structure consisting of a spherical core surrounded by a thin, protective semi-permeable polymeric shell. The release properties depend crucially on the nature of this coating structure [3]. For this reason, in recent years, materials innovation and nanotechnology have stimulated novel research and progress in biodegradable, biocompatible, environment-responsive, and targeted delivery systems [4, 5]. Mathematical modeling plays an important role in elucidating the drug release mechanisms, thus facilitating the development of advanced materials and the assessment of smart products by a systematic, rather than trial-and-error, approach. Within this area, special attention is given to semi-empirical and mechanistic models, which often provide a good fit with experimental data [6, 7].

Analysis and solution of mechanistic models taking into account the multi-layer structure of the microcapsule has been recently addressed by Kaoui et al. [8] and Carr and Pontrelli [9]. However, even in situations where the transport phenomena are well represented mathematically, the need exists to develop analytical tools or simple indicators to extract meaning from the model and the data. Such tools make it possible to answer practical questions, such as how long it takes to attain a therapeutic flux, or what processing conditions need to be adjusted and by how much, in order to reach a desired delivery rate without solving the complete mechanistic model. Furthermore, when a full mathematical representation is not available, or is too complicated to be of any practical use, simple performance indicators are required to elucidate the main transport mechanisms and identify the most critical components of the process. In most circumstances, for example, rather than the full transient solution, the time required to reach a steady state plasma drug concentration solely determines the effectiveness of the delivery system.

Among several mathematical techniques available for control systems analysis, Laplace transform and linearization techniques are frequently applied to describe a process’ dynamics through an effective time constant [10]: this represents a useful indicator of the time elapsed before reaching a steady state. The concept of a time constant is crucial to pharmaceutical developers who want to create controlled-release devices able to deliver drugs at a desired rate. Estimation of the release time, i.e. the time required for the drug-filled capsule to be depleted, is of great importance in the design of drug delivery systems, because it allows product manufacturers to tune specific properties to ensure a precise release within a determined time [11]. The idea of using a single timescale to characterize how fast the drug concentration reaches the equilibrium value has been explored by several authors. Pontrelli and Simon [12] define a timescale equal to the mean of a normalized probability density function representing the transition of the concentration profile from initial to equilibrium state. Calculating this mean requires either computing the full transient solution or evaluating the Laplace transform of the transient solution, and ultimately for multi-layered problems, produces quite complex expressions for the timescale in terms of the relevant parameters in the model. Another approach is the concept of mean action time, which has been used to characterize how long a diffusion process takes to reach steady state [13, 14, 15, 16]. Here, the transition from initial to equilibrium state is represented as a cumulative distribution function, with the mean of the corresponding probability density function defining the characteristic timescale. The attraction of the mean action time is that it completely avoids any calculation of the transient solution and produces simple explicit algebraic formulas for the timescale in terms of the parameters in the model [14, 15]. Finally, we note that the notion of release time is different but related to the concept of penetration time, which that has been studied in a number of configurations for thermal disturbances in heat conduction problems [17].

In this paper, inspired by recent literature on diffusion processes [18, 14, 19, 20, 15], we present two approaches that provide useful indicators of the release time for a drug-filled capsule. Both approaches utilize temporal moments of the concentration versus time curve at the centre of the capsule. The first approach, which uses the zeroth and first moments only, provides a simple and cheap way to characterize the release time in terms of the various parameters in the model (e.g. layered diffusivities, mass transfer coefficients, partition coefficients). The second approach yields an asymptotic estimate of the release time, defined as the time when the drug concentration at the centre of the capsule is a small prescribed distance away from its steady-state value. Attractively, both approaches can be used to avoid explicit (analytical or numerical) solution of the mathematical model and are useful to compare a priori different capsule configurations.

The remaining parts of this article are ordered in the following way: (i) the diffusion model for drug release from a spherical multi-layer capsule is presented (ii) both approaches for providing indicators of the release time are described (iii) a specific case study on a core-shell capsule is considered and numerical results presented and discussed.

2 Drug diffusion model for a multi-layer sphere

Multi-layer capsules consist of a drug-loaded (fluid or solid) spherical core surrounded by one or more polymeric layers [1, 2]. Such layer-by-layer assembly enhances a selective diffusion and allows for better control of the transfer rate. In this framework, and in the most general case of $n$ layers, the composite system can be treated as a sequence of enveloping concentric shells, constituted by spheres of increasing radii satisfying $0<R_{0}<R_{1}<\ldots<R_{n}$ (see Figure 1 for the case $n=1$). To prevent fast delivery, the capsule’s outmost layer is protected with a thin semi-permeable shell (coating). This coating shields and preserves the encapsulated drug from degradation and chemical aggression, protects the inner structure, and guarantees a more controlled and sustained release [21]. The coating is in contact with the targeted external release medium (either a bulk fluid or a tissue). Strictly speaking, for a pure diffusion problem into a homogeneous medium, the concentration field undergoes an exponential decay and vanishes asymptotically at infinite distance. In other words, at a given time, the concentration is gradually damped, reaching zero at infinite distance [22]. Nevertheless, for computational purposes, we can confine the diffusion process within an enveloping spherical layer of finite extent, at a distance “far enough” from the capsule surface. In our model, this additional layer is defined by setting a cut-off length $R_{\infty}\gg R_{n}$ (sometimes named release distance or diffusion length), beyond which the concentration reduces to zero, within a prescribed tolerance, at all times (Figure 1) [23]:

[TABLE]

Due to homogeneity and isotropy, we can assume that net drug diffusion occurs along the radial direction only, and thus we restrict our study to a one-dimensional model (Figure 1) that reflects a perfectly radially symmetric system. We adopt the model formulated in [9, 8], where the evolution of concentrations in the layers is governed by a set of one-dimensional linear diffusion equations:

[TABLE]

paired with interlayer and boundary conditions

[TABLE]

In the above equations, the parameters $D_{i}$ are the diffusion coefficients of drug in each layer, $\sigma_{i}$ are the partition coefficients, and $P$ is the mass transfer coefficient at the external coating [9, 8]. For a releasing drug-loaded core, the initial conditions are:

[TABLE]

where $C_{0}>0$ is a constant. The steady state solution of the drug diffusion model (2.1)–(2.10) is trivially given by

[TABLE]

In our previous studies the transient solution of Eqs (2.1)–(2.10) was obtained through a separation of variables [8], or a Laplace transform approach [9]. In some circumstances, however, rather than working with the complicated expressions of the full solution, simple and cheap measures of the performance of the delivery system are desired. In the following sections, we propose two ways to quantify the release time of the microcapsule, that is, the time taken for the capsule to be depleted of the drug.

3 Characterizing the release time

In this section, we define two indicators that provide a useful characterisation of the release time scale. Figure 2 shows the typical spatio-temporal behaviour of the dimensionless concentration, $c_{i}(r,t)/C_{0}$, arising from the solution of the diffusion model (2.1)–(2.10). An important observation from this plot is that the centre of the capsule, $r=0$, takes the longest time to reach steady state. In Figure 2b we plot the dimensionless concentration at the centre of the capsule, $\widetilde{c}_{0}(t):=c_{0}(0,t)/C_{0}$, versus time as it progresses towards the steady state solution of zero concentration (2.11). The area underneath the curve $\widetilde{c}_{0}(t)$, or the zeroth temporal moment of $\widetilde{c}_{0}(t)$, is defined as

[TABLE]

and shaded in Figure 2b. The value of $t_{c}^{(1)}$ will tend to be small for a fast release and large for a slow release. Moreover, when comparing two different capsule configurations, if the corresponding $\widetilde{c}_{0}(t)$ curves for the two configurations do not intersect then the configuration with the larger release time will have a larger area (larger value of $t_{c}^{(1)}$). Therefore, it is reasonable to conclude that $t_{c}^{(1)}$ provides a useful characterization of the timescale of release111We remark that Eq (3.1) differs from the standard definition of mean action time (see e.g. [14]) commonly used to characterize the time taken for a diffusion process to reach steady state. The reason for this difference is that the standard definition is not defined when $c_{i}(r,0)=c_{i}^{\infty}(r)$, as is the case in Eqs (2.8)–(2.11)..

A major attraction of working with Eq (3.1) is that a simple closed-form expression can be derived for $t_{c}^{(1)}$ involving the model parameters without requiring the full expression for $c_{0}(0,t)$. This is achieved by extending and modifying similar ideas presented elsewhere (see, e.g., [14, 15]). First, we define:

[TABLE]

which allow us to write down an equivalent form of Eq (3.1)

[TABLE]

Next, applying the linear operator $\mathcal{L}$, defined as

[TABLE]

to both sides of Eq (3.2) and making use of Eqs (2.1)–(2.3) yields the following differential equation:

[TABLE]

Supplementary boundary and interlayer conditions corresponding to Eqs (2.4)–(2.7) are derived using the definition of $u_{i}(r)$ in Eq (3.2) (see, e.g. [18, 14]). In summary, recalling the initial conditions (2.8)–(2.10) and the steady-state solution (2.11), we have the following boundary value problem satisfied by $u_{i}(r)$ ($i=0,1,\ldots,n,e$):

[TABLE]

The above boundary value problem admits a closed-form analytical solution. By way of example, we consider the simplest case of the core-shell model, a drug-filled core surrounded by one hydrogel shell ($n=1$). In this case, the differential equations (3.6)–(3.8) possess the general solution:

[TABLE]

where $\alpha_{0},\alpha_{1},\ldots,\alpha_{5}$ are arbitrary constants. Immediately, the boundary condition (3.9) requires $\alpha_{0}=0$. The remaining constants satisfy the following algebraic system generated by substituting Eqs (3.13)–(3.15) into the four interface conditions (3.10)–(3.11) with $n=1$:

[TABLE]

Eqs (3.16)–(3.20) can then be solved sequentially: Eq (3.17) for $\alpha_{2}$, Eq (3.18) for $\alpha_{4}$, Eq (3.20) for $\alpha_{5}$, Eq (3.19) for $\alpha_{3}$ and finally Eq (3.16) for $\alpha_{1}$. Combining these results with Eqs (3.13)–(3.15) yields the solution to the boundary value problem (3.6)–(3.12)222Note that evaluating either $u_{0}(r)$, $u_{1}(r)$ or $u_{e}(r)$ allows the area under the $\widetilde{c}_{0}(t)$ curve between $t=0$ and $t\rightarrow\infty$ to be calculated for any value of $r$, which may also be of practical interest.:

[TABLE]

Evaluating Eq (3.21) at $r=0$ gives the following expression for the characteristic timescale (3.3):

[TABLE]

In summary, Eq (3.24) provides a simple formula for characterizing the timescale of release: for small values of $t_{c}^{(1)}$ we expect a rapid release while for large values of $t_{c}^{(1)}$ we expect a slow release [14].

An alternative to $t_{c}^{(1)}$ is to incorporate the first temporal moment, $\int_{0}^{\infty}t\widetilde{c}_{0}(t)\,\mathrm{d}t$, into the calculation [15]:

[TABLE]

where the square root ensures that the second term has units of time. Note that $t_{c}^{(2)}>t_{c}^{(1)}$ with the inclusion of this additional term in Eq (3.25) introducing a penalty to $\widetilde{c}_{0}(t)$ curves with heavy tails that exhibit a slower decay to zero, behaviour that may not be captured by the zeroth moment alone. By deriving a similar boundary value problem to Eqs (3.6)–(3.12) for the first moment, as described later in Section 4, the following closed-form expression can be derived:

[TABLE]

The advantage of incorporating information about the first moment is that the expression (3.26), although more complex and costly than (3.24), provides a better indication of the release time (see section 5 for additional comments on the results). Other time scale indicators can be similarly defined by the natural extension of (3.25) to higher moments, but because of their increasingly complicated form, they likely have limited practical use.

4 Estimating the release time using high order moments

The characteristic timescales, $t_{c}^{(1)}$ (3.24) or $t_{c}^{(2)}$ (3.26), provide an indicator for, not an estimate of, the release time of the capsule. In this section, we present an asymptotic estimate of the release time defined as the time $t_{r}>0$ satisfying:

[TABLE]

where $\varepsilon$ is a small specified tolerance (see Figure 2b). In other words, $t_{r}$ measures the time taken for the capsule to be depleted within a small tolerance. Since $\widetilde{c}_{0}(t)$ monotonically decreases from one to zero, $t_{r}$ is unique for a given choice of $\varepsilon$. Under this definition, the release time is not an absolute concept, but dependent on the desired level of accuracy required to measure a complete depletion.

To derive the asymptotic estimate of the release time, we extend to multi-layer diffusion in spherical coordinates previous work by Carr [18], which focussed on monolayer diffusion in Cartesian coordinates. At $r=0$, the analytical solution to the drug diffusion model (2.1)–(2.10) possesses the following functional form [8]:

[TABLE]

where $\gamma_{j}$ and $\beta_{j}$ are constants. It follows then that the long time behaviour of the concentration is exponentially decreasing:

[TABLE]

Combining Eqs (4.1) and (4.3) gives the following asymptotic estimate of the release time333To guarantee $t_{r}>0$, the condition $\varepsilon<\gamma_{0}$ is required.:

[TABLE]

The formula (4.4) requires knowledge of $\gamma_{0}$ and $\beta_{0}$, which are related to the dominant eigenvalue and eigenfunction pair of the underlying Sturm-Liouville problem [8]. Alternatively, $\gamma_{0}$ and $\beta_{0}$ can be calculated by using appropriate temporal moments of the concentration as we now describe. Define the $k$th temporal moment444Under this definition the functions $u_{i}(r)$ defined in Eq (3.2) can be thought as the zeroth temporal moment, i.e., $u_{i,k}(r)$ with $k=0$.:

[TABLE]

where the first subscript indicates the layer ($i=0,1,...,e$) and the second subscript denotes the $k$th moment ($k=0,1,\ldots$). Combining Eqs (4.2) and (4.5) and carrying out the integration yields:

[TABLE]

Since $\beta_{0}<\beta_{j}$ for all $j=1,2,\ldots$, we have the following asymptotic relation for the higher moments [18]:

[TABLE]

Evaluating Eq (4.7) at $k$ and $k-1$ and solving the resulting algebraic system for $\gamma_{0}$ and $\beta_{0}$ gives:

[TABLE]

Finally, substituting Eq (4.8) into Eq (4.4) yields the following asymptotic estimate of the release time:

[TABLE]

where the superscript ($\ast$) is used to signify that $t_{r}^{\ast}$ is an asymptotic estimate of $t_{r}$. Due to Eq (4.7), the expectation is that $t_{r}^{\ast}$ becomes more accurate as $k$ increases [18]. The attraction of the formula for $t_{r}^{\ast}$ is that the moment expressions appearing in Eq (4.9) can be calculated without computing the full transient solution $c_{i}(r,t)$, which appears in the definition (4.5). This is achieved by deriving a similar boundary value problem to the one satisfied by $u_{i}(r)$ (or, equivalently $u_{i,0}(r)$) given in Eqs (3.6)–(3.12). Applying the linear operator $\mathcal{L}$ in Eq (3.4) to Eq (4.5) and utilising Eqs (2.1)–(2.3) yields:

[TABLE]

Integrating by parts and noting that $c_{i}(r,t)\rightarrow 0$ exponentially as $t\rightarrow\infty$ produces the following differential equation for $u_{i,k-1}(r)$:

[TABLE]

involving the $(k-1)$th moment, $u_{i,k-1}(r)$. Similar boundary and interlayer conditions to those in Eqs (3.6)–(3.12) apply [14] giving the following boundary value problem for the $k$th moment:

[TABLE]

The release time is computed iteratively by solving the sequence of boundary value problems (4.12)–(4.18) for increasing values of $k$ [18, 20]. First, an initial estimate of $t_{r}^{\ast}$ is calculated using $k=1$ in Eq (4.9) by solving Eqs (4.12)–(4.18) with $k=1$ and by using the zeroth order moment $u_{i,0}(r)$ computed previously (Section 3) in the right-hand side of Eqs (4.12)–(4.14). Next, an improved estimate of $t_{r}^{\ast}$ is computed using $k=2$ in Eq (4.9) by solving Eqs (4.12)–(4.18) with $k=2$ and using the previously computed moment $u_{i,1}(r)$ in the right-hand side of Eqs (4.12)–(4.14). The process repeats over $k$ until the value of $t_{r}^{\ast}$ converges sufficiently to an accurate estimate of $t_{r}$ as defined in Eq (4.1). We remark that numerically solving the system (4.12)–(4.18) is computationally inexpensive, since the coefficient matrix that arises from spatial discretisation of the boundary value problem remains unchanged throughout the iterations and only needs to be factorized once, with only the right-hand side of Eqs (4.12)–(4.14) changing with $k$.

5 Results and discussion

In this section, we apply the characteristic timescale of Section 3 and the release time estimate of Section 4 to several test cases. For all our numerical experiments we consider the drug diffusion model (2.1)–(2.10) with $n=1$ (core-shell model) and the following base values for the geometrical and physical parameters [9, 8, 21]:

[TABLE]

Using this set of parameter values, we first investigate the convergence behaviour of the asymptotic estimate of the release time, $t_{r}^{\ast}$ (4.9), for increasing values of $k$. In Table 1, we report the value of $t_{r}^{\ast}$ and $\widetilde{c}_{0}(t_{r}^{\ast})$ for $k=1,\ldots,14$ and $\varepsilon=10^{-4},10^{-5},10^{-6}$. These results show that $t_{r}^{\ast}$ is converging to a highly accurate estimate of the release time, with $\lim_{k\rightarrow\infty}\widetilde{c}_{0}(t_{r}^{\ast})$ close to the prescribed value of $\varepsilon$ in all three cases. For the two smaller values of $\varepsilon$, $\lim_{k\rightarrow\infty}\widetilde{c}_{0}(t_{r}^{\ast})$ is closer to the specified value of $\varepsilon$ because the one-term approximation (4.3) becomes more accurate for larger values of $t$. For this parameter set, the release time estimate converges to the second by $k=12$ for all three values of $\varepsilon$.

We now investigate the sensitivity of $t_{r}^{\ast}$ when varying the outer radius of the capsule, $R_{1}$, the mass transfer coefficient at the external coating, $P$, the diffusivity in the core, $D_{0}$, and the diffusivity in the hydrogel layer, $D_{1}$, about the base values given in Eqs (5.1)–(5.2). Each of the parameters, $R_{1}$, $P$, $D_{0}$ and $D_{1}$, are varied one at a time holding the other three parameters fixed at the values given in Eqs (5.1)–(5.2). We consider six logarithmically-spaced values of $R_{1}$ between $1.515\cdot 10^{-3}$ and $3\cdot 10^{-3}$, $P$ between $10^{-8}$ and $10^{-6}$, $D_{0}$ between $3\cdot 10^{-11}$ and $3\cdot 10^{-9}$ and $D_{1}$ between $5\cdot 10^{-12}$ and $5\cdot 10^{-10}$. In total, there are 24 test cases as labelled in Table 2. For all simulations, the release time estimate $t_{r}^{\ast}$ (4.9) is calculated with $k=15$ and $\varepsilon=10^{-4}$. Results in Figures 3(a)–(d) show the effect that varying each parameter has on the temporal profile of the dimensionless concentration at the centre of the capsule, $\widetilde{c}_{0}(t)$. These curves demonstrate that, over the specified ranges of each parameter, the time required for $\widetilde{c}_{0}(t)$ to approach zero increases for increasing values of $R_{1}$ and decreasing values of $P$, $D_{0}$ and $D_{1}$. This translates to larger values of $t_{r}^{\ast}$ for increasing values of $R_{1}$ and decreasing values of $P$, $D_{0}$ and $D_{1}$ (holding all other parameters fixed) as confirmed in Figures 3(e)–(h) and Table 2.

An interesting exercise is to determine how well the characteristic timescales, $t_{c}^{(1)}$ and $t_{c}^{(2)}$, indicate the value of the release time. To explore this, we plot $t_{c}^{(1)}$ and $t_{c}^{(2)}$ versus $R_{1}$, $P$, $D_{0}$ and $D_{1}$ in Figure 3(e)–(h) together with the corresponding values of the release time estimate $t_{r}^{\ast}$. The values of $t_{c}^{(1)}$ and $t_{c}^{(2)}$ are quite different from each other and differ considerably from $t_{r}^{\ast}$. However, this is not surprising because, as already mentioned, the goal of $t_{c}^{(1)}$ and $t_{c}^{(2)}$ is not to estimate the release time but rather provide a cheap way to rank different capsule configurations in terms of release time. Evident from these results is that both $t_{c}^{(1)}$ and $t_{c}^{(2)}$ follows the same general trend as $t_{r}^{\ast}$, that is, both increase with $R_{1}$ and decrease with $P$, $D_{0}$ and $D_{1}$. As the $\widetilde{c}_{0}(t)$ curves in Figures 3(a)–(d) do not intersect, we are guaranteed that a smaller value of $t_{c}^{(1)}$ implies a smaller value of $t_{r}^{\ast}$, as confirmed in Figures 3(e)–(h). An equivalent guarantee is true for $t_{c}^{(2)}$. Comparing across all 24 test cases listed in Table 2, we see that Case 24 gives rise to both the smallest characteristic timescales and release time while Case 7 yields both the largest ones.

The above results seem to indicate that the characteristic timescale always provides a useful indication of the release time. However, this is not necessarily true and caution is required when using either $t_{c}^{(1)}$ (3.24) or $t_{c}^{(2)}$ (3.26) as an indicator of the release time of different capsule configurations. For example, consider two configurations $A$ and $B$, where configuration $A$ has physical parameters (5.1)–(5.2) except with $P=5\cdot 10^{-8}\,\mathrm{m}\,\mathrm{s}^{-1}$ while configuration $B$ has physical parameters (5.1)–(5.2) except with $P=5\cdot 10^{-8}\,\mathrm{m}\,\mathrm{s}^{-1}$ and $R_{1}=2.2\cdot 10^{-3}\,\mathrm{m}$. In this case, the value of $t_{c}^{(1)}$ is larger for configuration $A$ despite the value of $t_{r}^{\ast}$ (calculated using $\varepsilon=10^{-4}$ and $k=15$) being smaller (see caption of Figure 4 for the numerical values). In Figure 4, we also see that the dimensionless concentration at the centre of the capsule $\widetilde{c}_{0}(t)$ for configurations $A$ and $B$ intersect, meaning that a smaller value of $t_{c}^{(1)}$ is not guaranteed to produce a smaller release time (as pointed out in Section 3). Conversely, the value of $t_{c}^{(2)}$ correctly correlates with $t_{r}^{\ast}$. Here, the second term in Eq (3.25) is larger for configuration $B$ due to its heavier tail, which causes the value of $t_{c}^{(2)}$ to be larger for configuration $B$ than for configuration $A$, consistent with $t_{r}^{\ast}$ (see caption of Figure 4 for the numerical values).

6 Conclusions

Polymer capsules, pellets, tablets, layer-by-layer vehicles and other micro-engineered drug releasing implants are attracting a great deal of attention for their potential use for therapeutic applications. The design of these novel drug delivery systems poses major challenges, such as the unknown significance of process parameters. Specifically, modelling and computational tools are provided to assist in the development of drug-release devices that can control the time to establish a steady-state flux or to deliver at a desired rate. The performance of a composite microcapsule can be sensibly enhanced if the release mechanism is understood and an appropriate mathematical model is used to characterize the releasing ability of the system. In this study, novel approaches to characterize and estimate the release time of diffusion problems from spherical multi-layer capsules are proposed under a limited number of physical assumptions. The method is based on the linearity of a pure diffusive system, that holds or dominates in most circumstances. No explicit solution of the diffusion problem is required; instead temporal moments of the drug concentration versus time curve at the centre of the capsule are used to derive analytical expressions that provide a priori quantitative indication of the drug release time. The proposed methodology, which combines in just one quantity the relevant geometrical and physical parameters, provides a simple tool to measure microcapsule dynamic performance.

Acknowledgments

G.P. acknowledges funding from the European Research Council under the European Unions Horizon 2020 Framework Programme (No. FP/2014-2020)/ERC Grant Agreement No. 739964 (COPMAT).