Total fraction of drug released from diffusion-controlled delivery systems with binding reactions

In diffusion-controlled drug delivery, it is possible for drug molecules to bind to the carrier material and never be released. A common way to incorporate this phenomenon into the governing mechanistic model is to include an irreversible first-order reaction term, where drug molecules become permanently immobilised once bound. For diffusion-only models, all the drug initially loaded into the device is released, while for reaction-diffusion models only a fraction of the drug is ultimately released. In this short paper, we show how to calculate this fraction for several common diffusion-controlled delivery systems. Easy-to-evaluate analytical expressions for the fraction of drug released are developed for monolithic and core-shell systems of slab, cylinder or sphere geometry. The developed formulas provide analytical insight into the effect that system parameters (e.g. diffusivity, binding rate, core radius) have on the total fraction of drug released, which may be helpful for practitioners designing drug delivery systems.


Introduction
Increasingly, mechanistic mathematical models of drug delivery (based on physical conservation laws) are being developed to improve understanding of the transport mechanisms that control the release rate, explore the effect of varying design parameters on the release profile and avoid costly and time consuming experiments [1].In this field of research, mechanistic mathematical models of diffusion-controlled drug delivery are typically based on Fick's second law, where the drug concentration within the system evolves in space and time according to the diffusion equation and specified initial and boundary conditions [2][3][4][5][6][7][8][9][10][11][12].Such purely-diffusive models yield a release profile F (t) (cumulative amount of drug released over the time interval [0, t] divided by the initial amount of drug loaded into the device) that increases from zero initially to one in the long time limit (Figure 1).Recent experimental work [13], however, has revealed that it is possible for drug to be trapped within the device and never be released due to drug molecules binding to the carrier material [14,15].This reduces the amount of drug delivered, resulting in a release profile that approaches a value less than one, denoted in this paper by F ∞ , in the long time limit (Figure 1).Several recent mechanistic models have incorporated binding into the governing diffusion model using an irreversible first-order reaction term, where drug molecules become permanently immobilised once bound [14][15][16][17].
In this short paper, we develop a set of analytical expressions for F ∞ for several diffusion-controlled delivery systems with first-order binding reactions.Our study is valid under standard assumptions: the system is radially symmetric (drug concentration is a function of radius and time only), the boundary of the system is stationary (swelling or erosion of the device is negligible) and the drug is initially completely dissolved (drug dissolution is instantaneous).Compact and elegant expressions for F ∞ are provided for both monolithic and core-shell systems of slab, cylinder or sphere geometry (Figures 2-3).The monolithic system consists of a single homogeneous carrier material (with constant diffusivity) while the core-shell system, consists of two homogeneous "layers" of carrier material (core and shell) with distinct diffusivities.For without binding with binding Figure 1: Typical drug release profile, F (t) (cumulative amount of drug released over the time interval [0, t] divided by the initial amount of drug loaded into the device), without binding reactions and with binding reactions.In the long time limit, the fraction of drug released Slab, cylinder and sphere devices of radius R consisting of a single homogeneous carrier material (with constant diffusivity).For the slab, R is assumed to be significantly smaller than both L and W so that drug release is negligible through the minor surfaces (of area 2RL or 2RW ) compared to the two major surfaces of the slab (of area LW ).For the cylinder, R is significantly smaller than H so that the drug release is negligible through the ends (of area πR 2 ) compared to the major surface of the cylinder (of area 2πRH).The initial amount of drug is assumed to be homogeneously distributed throughout the device and drug molecules may bind to the carrier material anywhere within the device.
Slab, cylinder and sphere devices of radius R consisting of two homogeneous "layers" of carrier material (core of radius R c encapsulated by a shell) with distinct diffusivities.For the slab, R is assumed to be significantly smaller than both L and W so that drug release is negligible through the minor surfaces (of area 2RL or 2RW ) compared to the two major surfaces of the slab (of area LW ).For the cylinder, R is significantly smaller than H so that the drug release is negligible through the ends (of area πR 2 ) compared to the major surface of the cylinder (of area 2πRH).The initial amount of drug is assumed to be homogeneously distributed throughout the core only and drug molecules may bind to the carrier material within the shell only.
the monolithic system, the initial amount of drug is assumed to be homogeneously distributed throughout the system while for the core-shell system, the initial amount of drug is assumed to be homogeneously distributed throughout the core only.For the monolithic system, drug molecules may bind to the carrier material anywhere within the device while for the core-shell system, drug molecules may bind to the carrier material within the shell only [15].Both the monolithic and core-shell systems are encapsulated in a thin coating (shell) that is either fully-permeable (no resistance to drug release) or semi-permeable (finite resistance to drug released) [8].In total, 12 distinct expressions for F ∞ are derived, one for each combination of system type (monolithic, core-shell), device geometry (slab, cylinder, sphere) and coating permeability (fully-permeable, semi-permeable).Each expression explains how the value of F ∞ changes for general values of the model parameters (e.g.diffusivity, reaction rate, surface transfer coefficient).
The remaining sections of this paper discuss how the expressions for F ∞ are developed for the monolithic and core-shell systems, respectively.Analytical and numerical evidence is then presented supporting the derived results.

Monolithic System
For the monolithic system, the drug dynamics are assumed to be governed by the reaction-diffusion model: where d specifies the geometry (d = 1, 2, 3 for the slab, cylinder and sphere, respectively), c(r, t) is the drug concentration, R is the device radius (Figure 2), D is the diffusivity, k is the binding rate, c 0 is the initial uniform drug concentration and P is the surface transfer coefficient.Note that the surface boundary condition (at r = R) depends on whether the coating is fully-permeable or semi-permeable with equivalence obtained in the limit as the surface transfer coefficient P tends to infinity.Mathematically, the total fraction of drug released F ∞ is given by the integral of the concentration flux over the release surface(s) divided by the initial amount of drug loaded into the device, which as shown in Appendix A.1 simplifies to under radial symmetry.In Appendix A.2, we show further that F ∞ can be expressed as where C(r) satisfies the boundary value problem: The attraction here is that the above boundary value problem admits exact analytical solutions involving hyperbolic (d = 1, 3) and Bessel functions (d = 2) that can be found by hand or using a computer algebra system.Solving the boundary value problem ( 6)-( 7) separately for d = 1, 2, 3, and evaluating the integral in (5) yields the expressions for F ∞ given below in equations ( 8)- (13).For d = 2, the expressions for F ∞ involve I 0 and I 1 , the modified Bessel functions of the first kind of zero and first order respectively.
Semi-permeable coating Sphere where λ = kR 2 D and P = P R D .
The above results reveal how F ∞ depends on a single dimensionless variable (λ) for the fully-permeable coating and two dimensionless variables (λ, P) for the semi-permeable coating.Studying these expressions, we see that equivalence between the semi-permeable and permeable cases is obtained when P → ∞.This fact together with the observation that the terms independent of P in the denominators of equations ( 11)-( 13) are positive, means that F ∞ is always smaller for the semi-permeable case than the corresponding fullypermeable case.This observation can be explained by the surface resistance resulting in drug molecules remaining in the device longer, increasing the likelihood that they bind to the carrier material and never be released.In summary, the above expressions provide the total fraction of drug released from the monolithic device (total amount of drug released divided by the initial amount of drug loaded into the device) for general values of the diffusivity (D), binding rate (k), device radius (R) and surface transfer coefficient (P ).Similar expressions to those above appear in the classical text [18] for the dual process of absorption into a cylinder/sphere from the surrounding medium.Finally, we note that the actual amount of drug released is obtained by multiplying F ∞ by the initial amount of drug loaded into the device: where L, W and H are defined in Figure 2.

Core-Shell System
For the core-shell system, the drug dynamics are assumed to be governed by the reaction-diffusion model: where d specifies the geometry (d = 1, 2, 3 for the slab, cylinder and sphere, respectively), c c (r, t) is the drug concentration in the core (0 < r < R c ), c s (r, t) is the drug concentration in the shell (R c < r < R), D c is the diffusivity in the core, D s is the diffusivity in the shell, k s is the binding rate (shell only), c 0 is the initial uniform drug concentration (core only) and P is the surface transfer coefficient.For the core-shell system, the total fraction of drug released is given by (see Appendix A.3): In Appendix A.4, we show that F ∞ can be expressed as: where C s (r) satisfies the following boundary value problem: As for the monolithic system, the above boundary value problem admits exact analytical solutions involving hyperbolic (d = 1, 3) and Bessel functions (d = 2) that can be found by hand (quite tediously) or using a computer algebra system.Solving the boundary value problem (21)-(24) separately for d = 1, 2, 3 and evaluating the integral in (20) yields the expressions for F ∞ given below in equations ( 25)-(30).For d = 2, the expressions for F ∞ involve I 0 and I 1 , the modified Bessel functions of the first kind of zero and first order respectively and K 0 and K 1 , the modified Bessel functions of the second kind of zero and first order respectively.

Fully-permeable coating
Sphere where λ s = ksR 2 Ds , and R c = Rc R .
Semi-permeable coating where The above results reveal how F ∞ depends on two dimensionless variables (λ s , R c ) for the fully-permeable coating and three dimensionless variables (λ s , R c and P s ) for the semi-permeable coating.Interestingly (perhaps unexpectedly), F ∞ is independent of the diffusivity in the core D c in all cases.As for the monolithic system, F ∞ is always smaller for the semi-permeable case with equivalence obtained when P → ∞.In summary, the above expressions provide the total fraction of drug released from the core-shell device (total amount of drug released divided by the initial amount of drug loaded into the device) for general values of the core diffusivity (D c ), shell diffusivity (D s ), shell binding rate (k s ), core radius (R c ), device radius (R) and surface transfer coefficient (P ).Finally, we note that the actual amount of drug released is obtained by multiplying F ∞ by the initial amount of drug loaded into the device: , where L, W and H are defined in Figure 3.

Supporting Numerical Experiments
The analytical expressions for F ∞ derived in the previous sections have been verified analytically using MATLAB's symbolic toolbox.Results in Figures 4 and 5 also provide evidence to support the derived expressions for a set of parameter values [16].These figures compare computed values of F ∞ to the long time limiting behaviour of the cumulative fraction of drug released: To compute F (t) we solve the governing reaction-diffusion model, equations ( 1)-( 3) for the monolithic system and ( 14)-( 18) for the core-shell system, numerically by (i) discretising over space [0, R] using a finite volume method with N r uniformly spaced nodes and (ii) discretising over time [0, T ] using the backward Euler method with N t time steps of fixed duration.Here T is chosen to be sufficiently large to capture the long time limiting behaviour of F (t).The numerical method computes discrete approximations to the concentration c(r k , t i ) for k = 1, . . ., N r and i = 1, . . ., N t , where r k = (k − 1)R/(N r − 1) and t i = iT /N t .The numerical approximations to c(r k , t i ) for k = 1, . . ., N r are then used to compute F (t i ) for i = 1, . . ., N t by applying a backward difference approximation to the radial concentration gradient in (31) followed by a trapezoidal rule approximation to the integral over [0, t i ] in (31).Comparisons in Figures 4 and 5 are given for N r = 2001 spatial nodes and N t = 10 4 time steps with excellent agreement reported between the computed values of F ∞ and F (T ) across all combination of system type (monolithic, core-shell), device geometry (slab, cylinder, sphere) and coating permeability (fully-permeable, semi-permeable).Full details on both the analytical and numerical verification are available in our MATLAB code, which is available on GitHub at https://github.com/elliotcarr/Carr2024a.

Conclusions
In summary, we have developed a set of analytical expressions to calculate F ∞ (total amount of drug released divided by the initial amount of drug) for several diffusion-controlled delivery systems in the presence of binding reactions.In total, 12 distinct expressions for F ∞ were derived, one for each combination of system type (monolithic, core-shell), device geometry (slab, cylinder, sphere) and coating permeability (fully-permeable, semi-permeable).Each expression provides analytical insight into the effect that particular physical and geometrical parameters (device geometry, diffusivity, binding rate, surface transfer coefficient) have on the value of F ∞ ; insight which may be helpful for practitioners designing drug delivery systems.While we have considered a wide-range of delivery systems in this work, it is important to note that the results are limited to the specific device configurations presented.Extending the analysis for the core-shell system to explore, for example, binding reactions in the core and/or non-zero initial drug concentration in the shell would yield different expressions for F ∞ .

A Appendices
A.1 Monolithic System: Initial form for F ∞ The total fraction of drug released, defined as the integral of the concentration flux through the release surface(s) (ω) scaled by the initial amount of drug loaded into the device region (Ω), simplifies to (4) as follows due to radial symmetry and the fact that |ω|/|Ω| equals the surface area of the (d − 1)-dimensional sphere of radius R divided by the volume of the d-dimensional ball of radius R (LW and H cancel from both the numerator and denominator for the slab and cylinder, respectively).
A.2 Monolithic System: Alternative form for F ∞ In this appendix, we demonstrate equivalence of the two forms of F ∞ given in equations ( 4) and (5).
Integrating the reaction-diffusion equation (1) over the slab, cylinder and sphere yields due to radial symmetry.Reversing the order of integration and differentiation on the left and evaluating the first integral on the right produces: and hence: Integrating this last equation over the time interval [0, ∞), using the initial condition (2) and the fact that c(r, t) tends to zero at infinite time yields: To complete the demonstration, we need to show that C s (r) satisfies the boundary value problem ( 21)-( 24).The differential equations ( 21)-( 22) are obtained by integrating the reaction-diffusion equations ( 14)-( 15) over the time interval [0, ∞), using the initial condition (16) and the fact that both c c (r, t) and c s (r, t) tend to zero at infinite time: and reversing the order of differentiation and integration as appropriate.