Mathematical modelling of variable porosity coatings for controlled drug release

TL;DR

This study models how variable porosity coatings with multiple layers can be engineered to optimize controlled drug release, showing that microstructural differences significantly influence release profiles.

Contribution

It introduces a mathematical model analyzing multi-layer variable porosity coatings, highlighting how layer properties can be tuned for desired drug release profiles.

Findings

Contrast in layer properties enhances control over drug release

Drug distribution across layers modulates early-stage release

Microstructural differences can be used to tailor release profiles

Abstract

In this paper we investigate the extent to which variable porosity drug-eluting coatings can provide better control over drug release than coatings where the porosity is constant throughout. In particular, we aim to establish the potential benefits of replacing a single-layer with a two-layer coating of identical total thickness and initial drug mass. In our study, what distinguishes the layers (other than their individual thickness and initial drug loading) is the underlying microstructure, and in particular the effective porosity and the tortuosity of the material. We consider the effect on the drug release profile of varying the initial distribution of drug, the relative thickness of the layers and the relative resistance to diffusion offered by each layer's composition. Our results indicate that the contrast in properties of the two layers can be used as a means of better controlling the release, and that the quantity of drug delivered in the early stages can be modulated by varying the distribution of drug across the layers. We conclude that microstructural and loading differences between multi-layer variable porosity coatings can be used to tune the properties of the coating materials to obtain the desired drug release profile for a given application.

Tables

Table 1. Table 1: Reference dimensional parameter values used in the baseline simulations. The two layers have the same physical parameters and unhindered transport between the layers, and so they are equivalent to one layer. ( ∗ ) (^{*}) In reality, we wish to impose K 2 , P → ∞ → subscript 𝐾 2 𝑃 K_{2},P\rightarrow\infty , however, for the purposes of the numerical simulations it was found that the value 10 10 superscript 10 10 10^{10} was sufficient to represent this case.

Parameter	Value (layer 1  $-$  layer 2)
$D_i^e\left(m^2{s}^{-1}\right)$	$5\cdot {10}^{-11}-5\cdot {10}^{-11}$
$L_i\left(m\right)$	$5\cdot {10}^{-5}-5\cdot {10}^{-5}$
${\varphi }_i$	$0.6-0.6$
$K_i$	$0-{10}^{10}{(}^{\ast })$
$P\left({\mathrm{ms}}^{-1}\right)$	${10}^{10}{(}^{\ast })$

Table 2. Table 2: Range of non-dimensional parameters simulated. In each Study, three values of δ ​ ( 0.2 , 0.5 , 0.8 ) 𝛿 0.2 0.5 0.8 \delta(0.2,0.5,0.8) were used. λ 1 m ​ i ​ n superscript subscript 𝜆 1 𝑚 𝑖 𝑛 \lambda_{1}^{min} is an indicator of the release time.

Study	$\chi$	$C^0$	$\varphi$	${\lambda }_1^{min}\left(\delta =0.5\right)$
Baseline	1	1	1	1.57
	0.5	1	1	1.16
1	2	1	1	2.03
	1	0	1	1.57
2	1	5	1	1.57
	1	1	2/3	1.61
3	1	1	3/2	1.50

Equations

$${\varphi }_1\frac{\partial c_1}{\partial t}$$

$${\varphi }_2\frac{\partial c_2}{\partial t}$$

$$-D_1^e\frac{\partial c_1}{\partial x}$$

$$-D_2^e\frac{\partial c_2}{\partial x}$$

$$-D_1^e\frac{\partial c_1}{\partial x}$$

$$c_1=c_1^0,-L_1\le x\le 0,t=0,$$

$$c_2=c_2^0,0<x\le L_2,t=0.$$

$$x^{\prime }=x/L,t^{\prime }=D_1^{e}t/{\varphi }_1{L}^2,c_1^{\prime }=c_1/c_1^0,c_2^{\prime }=c_2/c_1^0.$$

$$\frac{\partial c_1}{\partial t}$$

$$\frac{\partial c_2}{\partial t}$$

$$-\frac{\partial c_1}{\partial x}$$

$$\frac{\partial c_1}{\partial x}$$

$$-\chi \frac{\partial c_2}{\partial x}$$

$$c_1=1,-\delta \le x\le 0,t=0.$$

$$c_2=C^0,0<x\le 1-\delta ,t=0.$$

$$\delta =\frac{L_1}{L},\chi =\frac{D_2^e}{D_1^e},C^0=\frac{c_2^0}{c_1^0},\varphi =\frac{{\varphi }_2}{{\varphi }_1},$$

$$c_1(x,t)=X_1(x)G_1(t),c_2(x,t)=X_2(x)G_2(t).$$

$$\frac{G_1^{\prime }}{G_1}=-{\lambda }_1^2,\frac{\varphi }{\chi }\frac{G_2^{\prime }}{G_2}=-{\lambda }_2^2,$$

$$G_1(t)=\exp (-{\lambda }_1^{2}t),G_2(t)=\exp \left(-\frac{\chi }{\varphi }{\lambda }_2^{2}t\right),$$

$$X_1^{\prime \prime }$$

$$-X_1^{\prime }$$

$$X_1^{\prime }$$

$$X_2^{\prime \prime }$$

$$-\chi X_2^{\prime }$$

$$-X_1^{\prime }$$

$${\lambda }_1=\sqrt{\frac{\chi }{\varphi }}{\lambda }_2.$$

$$X_1(x)=a_1\cos ({\lambda }_{1}x)+b_1\sin ({\lambda }_{1}x),X_2(x)=a_2\cos ({\lambda }_{2}x)+b_2\sin ({\lambda }_{2}x),$$

$$a_1({\lambda }_1\sin ({\lambda }_1\delta )+{\Gamma }_1\cos ({\lambda }_1\delta ))+b_1({\lambda }_1\cos ({\lambda }_1\delta )-{\Gamma }_1\sin ({\lambda }_1\delta ))$$

$$a_2[-\chi {\lambda }_2\sin ({\lambda }_2(1-\delta ))+{\Gamma }_2\cos ({\lambda }_2(1-\delta ))]+b_2[\chi {\lambda }_2\cos ({\lambda }_2(1-\delta ))+{\Gamma }_2\sin ({\lambda }_2(1-\delta ))]$$

$$b_1{\lambda }_1$$

$$-b_1{\lambda }_1$$

$$\displaystyle\varphi(\lambda_{1},\lambda_{2})=\bigg{(}\lambda_{1}\sin(\lambda_{1}\delta)+\Gamma_{1}\cos(\lambda_{1}\delta)\bigg{)}\,\bigg{[}\Pi\bigg{(}\chi\lambda_{2}\cos(\lambda_{2}(1-\delta))+$$

$$\displaystyle\Gamma_{2}\sin(\lambda_{2}(1-\delta))\bigg{)}+\chi\lambda_{2}\bigg{(}-\chi\lambda_{2}\sin(\lambda_{2}(1-\delta))+\Gamma_{2}\cos(\lambda_{2}(1-\delta))\bigg{)}\bigg{]}-$$

$$\displaystyle\Pi\sqrt{\chi\phi}\bigg{(}\lambda_{1}\cos(\lambda_{1}\delta)-\Gamma_{1}\sin(\lambda_{1}\delta)\bigg{)}\bigg{(}-\chi\lambda_{2}\sin(\lambda_{2}(1-\delta))+\Gamma_{2}\cos(\lambda_{2}(1-\delta))\bigg{)}=0$$

$$a_2$$

$$a_1$$

$$b_1$$

$$X_{1m}$$

Full text

1 Introduction

The topic of drug delivery is a truly multi-disciplinary research area and has been attracting the interest of engineers, mathematicians, chemists and life scientists for decades. In particular, controlled drug delivery has received much attention, particularly concerning the design of tablets [1, 2, 3] and local drug delivery devices such as stents [4], transdermal patches [5], contact lenses [6] and orthopaedic implants [7] (Figure 1). Controlled release of drug from each of these vehicles can in principle be obtained by varying system design parameters. Some of the most common include the device geometry and materials; the physico-chemical properties of the drug and; the drug loading configuration. In the case of experimental studies, it is often demonstrated that different drug release profiles can be obtained by either varying the experimental conditions (e.g. in-vitro versus in-vivo) or physical delivery system properties, whilst in the case of mathematical and computational modelling, it is usual for a sensitivity analysis of the underlying model parameters to be conducted, and release profiles subsequently simulated. Both approaches are useful and indeed can be complementary in the quest for device design optimisation. In the case of tablets, there is a body of literature concerning multi-layer systems (see e.g. [1, 2, 3]), where the individual layers contain either different drugs or chemicals, or contrasting material properties from which the same drug or chemical is released in a bi- or multi-modal fashion.

However, the literature concerning multi-layer drug release from the aforementioned drug delivery devices is lacking somewhat, particularly in relation to mathematical modelling (see [12] as a rare exception). This will be the focus of the current manuscript.

Much of the present research concerned with such drug-eluting medical devices is focussed on developing sophisticated computational models which accurately simulate drug release and the subsequent distribution in the biological environment. The complexity of these models is increasing, with more and more realistic features being accounted for, including accurate 3D geometrical representations of the device and anatomical features; anisotropic and spatially-varying drug transport properties within the body and; complex features such as nonlinear binding reactions. If, on the one hand, these models are indeed necessary to accurately simulate drug transport within the device and in the biological environment, on the other hand it is clear that device manufacturers cannot intervene on the underlying biology. What they can control, however, are the properties of the device platform. Therefore, in this paper, we take a step back from the fully coupled computational models (see e.g. [13]) and focus instead solely on the properties of the drug-containing coating.

The drug is typically contained within some durable/biodegradable polymeric coating attached to the device platform or embedded within a nanoporous structure. The drug release profile depends on a number of factors including the porosity of the coating or bulk structure; the drug loading and initial distribution; the physico-chemical properties of the drug (e.g. molecule size, solubility, etc.) and; the release medium. A certain level of control is required: an excessive amount of drug delivered too quickly can result in toxicity whilst too little drug will have no therapeutic effect. However, the most desirable release profile is not always known and may in fact be patient-specific and therapy-dependent.

Motivated by today’s advances in material fabrication and by the increased capabilities of the miniaturization of structures offered by micro and nanotechnology, we propose variable porosity multi-layer coatings as an additional means of controlling the drug delivery and tailoring the release profile to the desired application. Our initial goal is to gain a better understanding of the potential benefits of replacing a single-layer with a two-layer drug-eluting coating of identical total thickness and initial drug mass. In our study, what distinguishes the layers (other than their thickness and initial drug loading) is the underlying microstructure, and in particular the effective porosity and tortuosity of the material (Figure 2). The structure of the paper is as follows. In Section 2 we provide the mathematical formulation of the problem and define a suitable non-dimensionalisation. We then propose, in Section 3, a semi-analytical solution method which makes use of separation of variables and expresses the solution as a Fourier series. A special case which admits an analytical solution is also presented. In the penultimate section we provide our results and investigate the sensitivity of the release profile to variations in the model parameters. Finally, in Section 5, we provide the conclusions of our study.

2 Mathematical formulation

A drug delivery device typically includes a polymeric matrix coating containing drug which is in contact with some release medium. The particular geometry of the device varies between applications, but the drug-eluting coating can usually reasonably be idealised as a slab (layer) of some thickness $L$. In Figure 2 we display an example of the situation we wish to model in the present work: two adjacent coating layers with different microstructural properties. Since the total thickness of drug-eluting coatings is typically small relative to the lateral coating dimensions, and the net drug transport is along a single direction, we restrict our attention to a one-dimensional model (Figure 3). We consider layers 1 and 2 to have thickness $L_{1}$ and $L_{2}$, respectively, with $L=L_{1}+L_{2}$ the total coating thickness, which we keep fixed in the following.

We represent each layer of the porous coating as a homogeneous material and define some representative elementary volume (r.e.v.) of size larger than the pore scale, but smaller than the typical length scale of the phenomenon. Within the r.e.v., we have solid and void parts. We choose to define all concentrations as intrinsically averaged variables, that is, averaged with respect to the void volume, rather than the total r.e.v.

Let $c_{1}$ and $c_{2}$ then denote the intrinsic concentrations of drug in layer $1$ of constant porosity $\phi_{1}$ and layer $2$ of constant porosity $\phi_{2}$, respectively. We further define $\phi_{i}^{e},\;i=1,2$ ($0<\phi_{i}^{e}\leq\phi_{i}$) as the effective transport-through porosities, which may be smaller than the overall porosity of each layer if, for example, there are small inaccessible pores or dead-end pores [15]. Additionally, we directly account for the fact that the molecules may have to travel through an increased path length due to the circuitous nature of the pores by introducing a tortuosity parameter $\tau_{i},\;i=1,2$. Assuming that the coating is rapidly wetted and that the drugs are readily soluble, it can be shown that drug transport satisfies the following diffusion equations

[TABLE]

where $D_{1}^{e}=\displaystyle\frac{\phi_{1}^{e}D^{w}}{\tau_{1}}$ and $D_{2}^{e}=\displaystyle\frac{\phi_{2}^{e}D^{w}}{\tau_{2}}$ are the effective diffusion coefficients in each layer and $D^{w}$ is the corresponding free diffusion of drug in water [16] . We emphasize that $D^{w}$ is independent of the microstructure and that we consider only the case of the same drug in each layer. In this work we envisage medical implant coatings which release drug through fluid-filled pores only. As a consequence, we do not consider diffusion in the solid phase, which can be several orders of magnitude slower than in the liquid phase.

For the sake of generality, we impose a mixed-type condition at both ends:

[TABLE]

where we may, in principle, choose $K_{1}$ and $K_{2}$ to match experimentally measured flux. The above boundary conditions allow us to explore the two extreme cases of zero flux and infinite sink conditions (see Section 4). If, for example, the coating is attached to an impermeable device (e.g. a stent) and drug release is measured under infinite sink conditions, we can let $K_{1}=0$ and $K_{2}\rightarrow\infty$.

At the interface between the two layers we impose continuity of flux. To keep the problem general, this flux accounts for a possible drug partitioning or a non-perfect contact, modelled through a mass transfer coefficient $P(m/s)$:

[TABLE]

We assume that initially the drug is loaded at uniform concentrations $c_{1}^{0}$ and $c_{2}^{0}$ in layers 1 and 2, respectively:

[TABLE]

The case of a single layer can easily be recovered, as will be demonstrated in Section 4.

2.1 Non-dimensionalisation

We now proceed to non-dimensionalise equations (2.1)-(2.7). We choose

[TABLE]

The non-dimensionalised equations (after dropping primes) are then:

[TABLE]

[TABLE]

where

[TABLE]

We note that the non-dimensional parameter $\chi$ contains all the important microstructural parameters which influence drug release.

3 Solution procedure

3.1 Solution by Separation of Variables

The model given by (2.8)-(2.14) is amenable to solution by separation of variables, an approach we have adopted in previous work considering two-layer and multi-layer problems [18, 17]. We let

[TABLE]

Equations (2.8)–(2.9) give rise to the ordinary differential equations (ODEs):

[TABLE]

which yield the solution:

[TABLE]

and the Sturm-Liouville eigenvalue system:

[TABLE]

[TABLE]

obtained by setting $G_{1}=G_{2}$, which implies

[TABLE]

The general solution of the ODEs (3.4) and (3.7) is:

[TABLE]

where the eigenvalues $\lambda_{i}$ ($i=1,2$) and the unknown coefficients $a_{i}$ and $b_{i}$ may be computed by imposing the boundary and interface conditions as follows. From equations (3.5) and (3.8), we have:

[TABLE]

From the interface conditions (3.6) and (3.9), it follows:

[TABLE]

Equations (LABEL:sec41)–(3.15) form a system of four homogeneous linear algebraic equations in the four unknowns $a_{1},b_{1},a_{2}$ and $b_{2}$ . To obtain a solution different from the trivial one $(0,0,0,0)$, it is necessary that the determinant of the coefficient matrix associated with the above system is equal to zero, that is:

[TABLE]

By replacing $\lambda_{1}$ with $\lambda_{2}$ through the relation (3.10), if the above transcendental equation (eigen condition) in $\lambda_{2}$ is satisfied, the coefficients may be taken as:

[TABLE]

where the multiplicative constant $b_{2}$ is arbitrary and its value depends on the initial condition (see below). We note that $\varphi$ depends on the parameters $\Pi,\delta,\chi,\phi,\Gamma_{1},\Gamma_{2}$ (but not $C^{0}$) and has infinitely many roots (eigenvalues), which are real and distinct.

For each eigenvalue couple $(\lambda_{1m},\lambda_{2m}),\>m=0,1,2,...,$ satisfying (3.16), the constants $a_{1m}$, $b_{1m}$ and $a_{2m}$ are obtained from (3.18), (3.19) and (3.17) respectively, and thus the corresponding eigenfunctions $X_{1m}$ and $X_{2m}$ defined in (3.11) are computed as:

[TABLE]

where the tilde indicates a variable which has been scaled by $b_{2m}$.

Furthermore, the corresponding time-variable functions $G_{1m}$ and $G_{2m}$ defined by equations (3.3) are computed as:

[TABLE]

($G_{1m}=G_{2m}$). Finally, the complete solution of the problem is given by a linear superposition of the fundamental solutions (3.1) in the form:

[TABLE]

where the arbitrary constants $A_{m}(=b_{2m})$ are determined through the initial conditions (2.14). The damping factors $\exp(-\lambda_{1m}^{2}t)$ and $\exp\left(-\displaystyle{\chi\over\phi}\lambda_{2m}^{2}t\right)$, $m=1,2,...,$ measure the attenuation of the various terms in summations (3.23). Because of the fast exponential convergence, the series (3.23) will be truncated at a finite number of terms, in accordance with the accuracy desired at the time of interest. Since $\max_{x}|A_{m}\tilde{X}_{im}(x)|<1$ for any $i=1,2$, $m>1$, to reach an accuracy of $10^{-r}$, it is sufficient to consider a finite series summation up to the index $j>1$ such that

[TABLE]

and the series is truncated at the first $j$ terms. A value of $j=30$ is considered for all times in the simulations.

Application of the initial condition

By evaluating (3.23) at $t=0$ and multiplying it by $\tilde{X}_{1n},\tilde{X}_{2n}$, after integration we obtain:

[TABLE]

and

[TABLE]

By combining equations (3.24) and (3.25) and by using the orthogonality property of $(X_{1m},X_{2m})$ [18] :

[TABLE]

we have:

[TABLE]

3.2 Computing Mass

The total mass of drug at any time can be evaluated by integrating the drug concentrations in each layer over their respective spatial domain. If we normalize the total mass by its initial value, then the non-dimensional total mass of drug in the coating is given by

[TABLE]

Letting $\theta_{i}$ represent the non-dimensional mass of drug in each layer as a fraction of the total mass, we then have:

[TABLE]

where $c_{i}$ ($i=1,2$) are given by (3.23). It is then straightforward to evaluate the total non-dimensional mass of drug in the coating as $M=\theta_{1}+\theta_{2}$ and the cumulative fraction of drug release, $M_{frac}$ (the release profile) as $M_{frac}=1-\left(\theta_{1}+\theta_{2}\right)$.

The depletion of the drug in coating as a result of the release process is governed by an exponential decay as in the above equations. The analytical solution indicates that a complete release is reached only asymptotically and equations (3.29)–(3.30) allow one to estimate the release time $T_{r}$, within a given tolerance $\epsilon$, through:

[TABLE]

The smallest eigenvalue $\lambda_{1}^{min}=\displaystyle\sqrt{\chi\over\phi}\lambda_{2}^{min}$ relates to the dominant damping factor in the series (3.23), (3.29)-(3.30). Comparing $\lambda_{1}^{min}$ between different parameter regimes provides an indication of the relative rate of release (see Table 2). For the particular parameter regime of interest and a given initial mass per cross-sectional area $M^{0}$, the initial loading concentrations are calculated through:

[TABLE]

3.3 Special Case

We note that in the special case where the microstructural properties of the layers are identical ($\chi=\phi=1$) then we can obtain an analytical solution. In this case, the eigenvalues $\lambda_{1}=\lambda_{2}=\lambda$, say, and are obtained by solving

[TABLE]

The difference between the solutions $c_{1}$ and $c_{2}$ then arises only through $A_{m}$, which is calculated using the initial condition. It can be shown that the solution in this case is

[TABLE]

and the release profile is readily computed as:

[TABLE]

4 Results and Discussion

In all simulations, we consider the typical boundary conditions $\Gamma_{1}=0$ and $\Gamma_{2}\rightarrow\infty$. This situation is representative of the most common case where the coating is in contact with an impermeable material (e.g. metal structure of the device) on one side and is exposed to an infinite sink at the other side, where the drug is washed away instantaneously. It is usual for infinite sink conditions to be maintained during in-vitro drug release experiments. At the interface between the two layers, we choose $\Pi\rightarrow\infty$ to reflect ‘unhindered’ transport.

The free diffusion coefficient of molecules in liquids $D^{w}$ is typically of the order of $10^{-9}\;m^{2}\;s^{-1}$ [16]. By definition, $0<\phi_{i}^{e}<1$. However, extremely low ($\phi_{i}^{e}<0.1$) and extremely high ($\phi_{i}^{e}>0.9$) porosities would likely result in drug loading and mechanical constraints, respectively. A typical range of tortuosity values is $1<\tau_{i}<6$, although values as high as 10 have been reported [16]. Taken together, we expect that the effect of the microstructure in each layer is to result in an effective diffusion coefficient at most two orders of magnitude smaller than the free diffusion coefficient in water. We note that drug diffusion coefficients in some polymers have been reported to be as low as $10^{-17}\;m^{2}\;s^{-1}.$ However, it should be noted that these are usually apparent diffusion coefficients which likely incorporate other effects such as absorption and desorption (possibly in addition to the microstructure effects that we consider here). The effect of such processes can be to reduce the overall diffusion coefficient by several orders of magnitude. In all simulations we fix $D_{1}^{e}=5\cdot 10^{-11}$ and $\phi_{1}=0.6$ and consider the effects of varying the microstructure of each layer by varying $\phi$ and $\chi$.

Since the purpose of this study is to establish the benefits of replacing a single layer with a two-layer coating of identical total thickness $L$ and initial drug mass, we fix $L=10^{-4}m$ in all simulations. In reality, of course, the values of $D_{1}^{e}$, $\phi_{1}$ and $L$ will vary depending on the particular application. Since our focus is to investigate the effect of the results on varying the ratio between the parameters of each layer, we have decided to choose broadly typical values, whilst acknowledging that this will not cover all cases. In all of our simulations the initial non-dimensional mass is 1. We choose not to mathematically implement a fixed dimensional mass. As a consequence, the initial dimensional loading concentrations $c_{1}^{0}$ and $c_{2}^{0}$ are to be back-calculated using (3.31) such that the desired initial dimensional mass is achieved.

4.1 Baseline model

To assess the effect on drug release of variations in system parameters, we preliminarily assume that layer 1 and layer 2 have identical microstructural parameters ($\chi=\phi=1$) and equal initial drug concentrations ($C^{0}=1$): in this case we can use the analytical solutions given by (3.3-3.3). The result is that our baseline model (see Table 1) essentially reduces to a single layer system (the solution is independent of the choice of $\delta$). The resulting non-dimensional parameters are $\chi=\phi=C^{0}=1$.

In Figure 4 we display the results of the baseline case, where each layer has identical initial drug loading and microstructure, so that we effectively have a single layer. As a result of the infinite sink boundary condition at the release medium, drug is rapidly released from layer 2 in the early stages, whilst there is a small delay before drug concentrations in layer 1 drop from their initial value. Drug release from layer 1 proceeds at a slower rate than in layer 2, and therefore there is a difference in both the shape and the duration of release in each layer. All of the drug has been released from the system by approximately $t=3$ (non-dimensional time).

4.2 Sensitivity analysis

We are ultimately interested in quantifying the effect on drug release of having two separate layers (as opposed to a single layer) with different microstructure and drug loading parameters. Therefore, it is of interest to vary the parameters, one at a time, around the baseline values and to compare the resulting drug release profiles. We consider three cases (see Table 2): in Study 1 we assess the effect of varying $\chi$, whilst in Study 2 and Study 3 we vary $C^{0}$ and $\phi$, respectively. In each case we consider three values of $\delta$.

Study 1: effect of varying microstructure ratio $\chi$

We now assess the effect of varying the relative microstructural parameters between the two layers. In Figure 5, left column, we choose $\chi$ such that the effective diffusion coefficient in layer 2 is half that of layer 1, whilst in the right column the effective diffusion coefficient is 2 times greater. In the first case we observe that drug release from layer 1 is hindered by the lower effective diffusion coefficient in layer 2 and as a result there is a delay in drug being released from layer 1. Despite the lower effective diffusion coefficient in layer 2, there is still a burst release as a result of the infinite sink conditions, but this effect is smaller than the baseline case. In the second case, the faster effective diffusivity in layer 2 results not only in significantly faster drug release from layer 2, but also from layer 1 (Figure 5, right). In Figure 6 we plot the overall release profiles in these two cases and compare with the baseline (dashed red line). We display only the case of $\delta=0.5$. From Figure 6 we conclude that the parameter $\chi$ has a strong influence on both the shape (rate of release) and the duration of release. This is perhaps unsurprising since $\chi$ appears prominently in the exponential damping factor (see (3.23) and (3.30)). For $\chi=0.5$ and $\chi=2$ we observe faster release and slower release, respectively, as we increase $\delta$ (not shown).

The implication here is that, simply by varying the microstructure of the two layers, not only it is possible to alter the shape of the release profile, but it is also possible to ensure that drug is delivered over some defined period of time. We note that although $\chi$ contains parameters relating to both the porosity and tortuosity of each layer, it is the combination of these values (rather than their individual size) which defines the release profile. For example, a value of $\chi=2$ could be obtained by doubling the effective porosity of layer 2 (in comparison with layer 1) or, by doubling the tortuosity of layer 1 (in comparison with layer 2). Therefore, this parameter is highly important as it offers much flexibility from the manufacturing point of view.

Study 2: effect of varying ratio of initial concentrations $C^{0}$

We now elucidate the effect of varying the initial drug concentration between the two layers. In the first case (Figure 7, left) we choose the initial drug concentration in layer 2 to be zero, whilst in the second case (Figure 7, right) we choose the concentration in layer 2 to be five times that of layer 1. In the first case we observe that layer 2 is initially rapidly infiltrated with drug, before drug is subsequently released after it has traversed the thickness of the second layer. In the second case we observe that whilst layer 2 is depleted rapidly as a result of the infinite sink condition, at early times an increase in drug concentration (and consequently drug mass) is observed in layer 1 due to the concentration gradient between the two layers (we are assuming that no drug can diffuse between the layers prior to the coating being placed in the release medium). As layer 2 continues to be depleted of drug, eventually the concentration gradient at the interface changes direction and drug then diffuses from layer 1 into layer 2 before being released. In each case, the value of $\delta$ has a significant impact on the concentration profile in each layer.

In Figure 8 we plot the overall release profiles in these two cases with $\delta=0.5$ and compare with the baseline (dashed red line). For both $C^{0}=0$ and $C^{0}=5$ we observe faster release as we increase $\delta$ (not shown). We conclude that having a drug-free second layer can delay the start of the drug-release process, which may be desirable in certain applications. In contrast, choosing a higher initial drug concentration in the second layer can result in a larger burst of drug which also may be advantageous in other circumstances. However, in all cases the overall duration of release results the same. Therefore, the non-dimensional parameter $C^{0}$ can be used as a tuning parameter to vary the proportion of drug delivered in the initial stages. The inflection point at $t=0$ (Figure 8, $C^{0}=0$ - black curve) indicates a retardation time due to the filling of the second layer which is initially empty.

Study 3: effect of varying porosity ratio $\phi$

By varying $\phi$ from $2/3$ to $3/2$, no significant differences are found between the concentration profiles (not shown). The release profiles are vitually indistinguishable for high values of $\delta$, although minor differences in the profiles are observed as $\delta$ is reduced (not shown). These results reinforce the idea that it is the effective porosity in each layer $\phi_{i}^{e}$ that drives the drug release, rather than the overall porosity $\phi_{i}$.

Having studied separately the influence of the individual parameters, we note that a combination of the above cases should be considered in order to meet the precise manufacturing requirements or with the aim of optimising some quantity. For example, if the objective is to slow down the release, then it appears that the simultaneous occurrence of the two cases $\chi<1,\;C^{0}=0$ will boost this property: in particular, a configuration with a lower effective porosity in layer 2 faced with one of higher effective porosity in layer 1 acts more favourably to achieve this goal: the time scale for release from layer 2 is increased, and layer 1 acts as a reservoir that continuously supplies drug during elution.

5 Conclusions

In this paper we have presented a mathematical model of drug diffusion through two adjacent porous layers and we have carried out a systematic study of the effect on drug release of changes to system parameters. Our results indicate that the contrast in properties of the two layers can be used as a means of better controlling the release, and that the quantity of drug delivered in the early stages can be modulated by varying the distribution of drug across the layers. We conclude that microstructural and loading differences between variable porosity coating layers can be utilized to tune the properties of the coating materials to obtain the desired drug release profile for a given application. We expect that our results will generalise to the multi-layer case, with increasing numbers of layers exhibiting contrasting properties potentially providing additional flexibility for targeting a specific release profile. Finally, as we reduce the thickness of each layer, in the limit we can obtain a continuously changing porosity. Whilst we acknowledge that this may provide even more flexibility in terms of controlling the release, the model we consider is a useful starting point to assess the effect of variable porosity. Additionally, the manufacturing of such a system is likely to be very challenging given the typical coating thicknesses.

We would like to emphasise that we have made a number of simplifications in this work. Perhaps the most significant is the assumption that drug is transported via a diffusive mechanism only. Depending on the particular coating material under consideration, it may be more appropriate to account for: polymer-drug interactions; diffusion through the solid phase; erosion; swelling and/or degradation. Additionally, in cases where fluid penetration into the coating is slow and/or the drug in question is poorly soluble, then the model may need to account for the drug dissolution process. Nevertheless, the model we have provided here will be relevant in a number of drug delivery cases. Now that we have established that variable porosity coatings for drug-eluting devices are worth further consideration, we will seek to study some of these additional features in future work.

Acknowledgments

Sean McGinty would like to acknowledge the funding provided by EPSRC under grant nos. EP/J007242/1 and EP/J007579/1. David King would like to acknowledge funding provided by EPSRC under grant number EP/M506539/1.