Modelling functionalized drug release for a spherical capsule

Advances in material design has led to the rapid development of novel materials with increasing complexity and functions in bioengineering. In particular, functionally graded materials (FGMs) offer important advantages in various fields of application. In this work, we consider a heterogeneous reaction-diffusion model for an FGM spherical drug releasing system that generalizes the multi-layer configuration to arbitrary spatially-variable coefficients. Our model proposes a possible form for the drug diffusivity and reaction rate functions exhibiting fixed average material properties and a drug release profile that can be continuously varied between the limiting cases of a homogeneous system (constant coefficients) and two-layer system (stepwise coefficients). A hybrid analytical-numerical solution is then used to solve the model, which provides closed-form expressions for the drug concentration and drug release profiles in terms of generalized Fourier series. The resulting concentration and mass profiles show how the release rate can be controlled and continuously varied between a fast (homogeneous) and slow (two-layer) release.


Introduction
Spherical drug carriers are among the most common formulations for a controlled release system.In particular, microcapsules are small spherical particles produced by coating templates constituted of different polymers and using various fabrication strategies [1].Although they can be made of a variety of sizes and materials, capsules of interest for most bio-applications have diameters ranging from some nanometers to a few micrometers.Specific examples include liposomes, pellets, nanocontainers, and others [2].The effectiveness of polymeric delivery systems can be improved by designing structures with modified material properties that are capable of responding to specific pre-set conditions that prescribe the release of the loaded drug.
One of the approaches generally recognized as effective in the assembly of polymer particles is the layer-by-layer technique.These drug carriers are considered as challenging releasing devices because of their unique multi-layer structural properties [3].Another family of drug delivery systems is constituted Homogeneous Two-layer g FGM Figure 1: Schematic representation of an FGM with continuous variation of porosity/density compared to an homogeneous material (no variation) and a two-layer material (stepwise variation).by stimuli-responsive capsules that control the release of the therapeutic active agents in response to external triggers such as temperature, pH and many others [4,5].Among other concurrent effects, such as dissolution, polymer swelling and possible degradation, diffusion remains the most important mechanism used to control the release rate from drug delivery systems [6,7].Some experiments, however, show that a fraction of the initial drug loaded is retained within the shell and is never released, due to the specific capacity of the polymer to permanently bind the drug molecules [8].It is common to model this observed phenomenon through first order reaction kinetics [4,8].
Mathematical models of drug release from spherical carriers provide insights of mass transport and drug kinetics involved in drug delivery as well as the effect of design parameters, such as the device geometry and drug loading distribution, on the release mechanism and can significantly reduce the number of experimental studies [9].However, these models depend on so many variables and parameters that, if not appropriately simplified, can raise more questions than useful answers.Analysis of diffusion-controlled system are confined to homogeneous spheres where an exact solution is available [10], or to layered capsules [11,12], where various mathematical models have been proposed to describe the drug release from this system over the years [1,6,9,13,14].
While such configurations are well understood, there is still room for improvement in mechanistic models to control the release mechanism from a drug-loaded sphere.For example, the effect of nonhomogeneity represents an important feature that can influence greatly the release properties.Functionally graded materials (FGMs) are a variety of composite materials in which the material properties vary smoothly and continuously (Fig 1).This is in contrast to previous approaches for achieving varying material properties, such as layer-by-layer assembly, where there is an abrupt change in properties from one layer to the next.FGMs, i.e composite materials that have a progressive compositional gradient, are already currently used in a wide range of applications [15,16].Using today's micro-engineering potential, it is possible to manufacture and control the material properties of the substrate to have the desired smart release properties [2].For example, new possibilities are derived from 3D printing technology to manufacture material micro-porosity and density in non-homogeneous substrates [17,18].
A mathematical model of drug release from a thin film FGM has been recently proposed and solved numerically [19].In the current work, we propose a reaction-diffusion continuum model to describe drug transport within, and release from, a drug-loaded FGM spherical capsule and develop a hybrid analytical-numerical eigenfunction expansion solution [20,21,22,23,24] to handle the spatiallyvariable coefficients.Our model proposes a possible form for the drug diffusivity and reaction rate functions, which exhibits the same average material properties and can be continuously varied between the limiting cases of a homogeneous system (constant coefficients) and two-layer system (stepwise coefficients)

Max
The rest of the article is organized as follows.In the next section, we present the model equations and boundary conditions that govern the drug mass release from a non-homogeneous FGM spherical system.In section 3, we present the hybrid analytical-numerical solution methodology leading to a closedform solution.In section 4, drug concentration and drug release profiles are presented for two distinct cases: pure diffusion and reaction-diffusion.Through extensive simulations, we explore the effect of FGM systems on the drug release mechansim by comparing the release profiles to those obtained from standard homogeneous and two-layer systems.Finally, section 5 provides general conclusions and some perspectives for future studies.

Using FGM in releasing spherical particles
Drug nanocontainers and releasing microcapsules are the subject of considerable research effort because of their structural and morphological properties, allowing the synthesis of materials capable of responding to biochemical alterations of the environment [4].Particularly, layer-by-layer polymeric releasing particles have gained increasing interest for their ability to control and tune the release of one or more therapeutic drugs [25].Here, the layers are constituted of different materials having specific physico-chemical characteristics and are customized to allow a selective diffusion and better control the transfer rate [1].In the layer-by-layer configuration, a semi-permeable external shell (coating) is often designed to shield and preserve the encapsulated drug from degradation and chemical aggression, and guarantee a more controlled and sustained release [26].With the aim of overcoming and generalizing the layered structure, we explore the potential of a material with continuously changing properties.
Recently, more attention has been paid to the class of functionally graded materials (FGMs), in several fields of application [16].FGMs are a special kind of composite materials in which the microstructural properties vary smoothly and continuously in space [15].In a purely diffusive model, the continuously varying nature of FGMs naturally lends itself to different functional forms of the diffusion coefficient D(x) in the domain Ω: where c(x, t) is the mass volume-averaged concentration of drug.We assume the diffusivity is higher at lower polymer density (inner region) and lower at higher polymer density (outer region), to account for a material that gradually thickens outwards.
From experiments, however, it is observed that a fraction of the initial loaded drug is retained and never released.A possible explanation of this phenomenon is a chemical reaction due to polymerdrug interaction.In other words, due to long polymeric chains and possible electrostatic interactions, a small percentage of the initial loaded drug is entrapped without being released [4].We model this phenomenon by using reaction kinetics, where the drug molecules travelling through the polymer can potentially be permanently bound with a rate k(x) [8].This generalizes equation (2.1) to include a first-order reaction term: where k(x) [s −1 ] is a space-dependent reaction rate.

Drug release from a FGM sphere
We consider a reservoir-type drug carrier with an active agent loaded in a spherical polymeric matrix, which is one of the most common formulations for a controlled release system.The spherical carrier is assumed to have radius R giving domain Ω = {x ∈ R 3 | x < R} and outer surface ω = {x ∈ R 3 | x = R}, where • is the Euclidean norm.The interior of the spherical capsule is made of a nonhomogeneous FGM, reflecting a customized composition that allows for selective spatially-dependent diffusion and reaction to better control the drug transfer rate (Fig 2 ).The case of a core-shell capsule is also included in this model, through stepwise diffusivity and reaction functions.As for in-vitro experiments, the sphere is immersed in an external ambient medium of a large extent (relative to size of the sphere), taken as semi-infinite.
In the case of an isotropic sphere centered on the origin with a boundary condition on its outer surface, we can assume that net drug diffusion occurs along the radial (r) direction only, and thus we restrict our study to a one-dimensional model, as follows: The model model permits a space dependent initial concentration (2.4) and accounts for a flux resistance (with mass transfer coefficient P ) at the external surface (2.6) due to the semi-permeable coating [12].
The choice of D(r) and k(r) A specific form for D(r) is needed to characterize the material heterogeneity of the FGM.In particular, we assume that the physical medium properties (density or porosity) may change continuously along the radius, being softer in the core and harder towards the surface, leading to a decreasing function D(r) (cf.Fig 2) that varies from a maximum possible value of D max at the center (r = 0) to a minimum possible value of D min at the surface (r = R).Among a variety of feasible continuous diffusivity functions, we consider This choice is a standard smooth approximation to a two-layer stepwise diffusivity where H(•) is the Heaviside function at r = σ with α > 0 [-] inversely related to the width of the transition layer and σ [cm] denoting the location of the transition center.
Our aim is to understand the effect of varying α and σ on the drug release profile.To maintain an FGM with the same average diffusive properties, for a specified choice of α, we calculate a corresponding value of σ so that the average value of D(r) over the spherical capsule is constant: in spherical coordinates, when using radial symmetry and V (Ω) = 4πR 3 /3.In our results, we set which is the unique value that yields σ → R/2 in the limiting case of a two-layer stepwise medium (α → ∞).In summary, for a specified choice of α, we calculate σ by solving the nonlinear equation: As a consequence of the decreasing diffusivity towards the external surface (due to the thickening material) the drug reaction rate k(r), which is typically proportional to the polymer density/porosity, undergoes a similar radial variation as D(r), but in the opposite direction (cf.Fig 2), resulting in an increasing function from a minimum possible value of k min at the centre (r = 0) to a maximum possible value of k max at the surface (r = R) : where the values of α and σ are the same as those used for D(r).As for D(r), the average value of k(r) over the full spherical capsule is constant:

Solution methodology
We solve the heterogeneous diffusion model (2.3)-(2.6)by introducing the following non-dimensional variables: where c 0 (r) and considering the resulting non-dimensional analogue of equations (2.3)-(2.6): Equations (3.4)-(3.7)constitute a linear problem with spatially-variable coefficients and homogeneous boundary conditions.To solve this problem, we use a hybrid analytical-numerical approach where ĉ(r, t) is expanded in terms of eigenfunctions: T n ( t)X n (r). (3.8) The space function X n (r) In the solution expansion (3.8), we let X n (r) be the eigenfunctions associated with the following Sturm-Liouville problem: The eigenvalues, denoted by λ n for n ∈ N + , are defined as the positive values of λ satisfying (3.14) with the corresponding eigenfunctions given by: which are orthonormal on the interval [0, 1]: The time function T n ( t) In the solution expansion (3.8), the time functions T n ( t) are computed by imposing that (3.8) satisfy the actual governing equation (3.4) with space dependent D(r) and k(r).We first substitute (3.8) into (3.4) and differentiate to give T n ( t)X n (r).
Next, multiplying both sides of this equation by r2 X m (r) and integrating from r = 0 to r = 1, we see that T m ( t) satisfies the following differential equation after making use of the differential equation (3.9) and orthogonality (3.16).Equation (3.17) identifies where X n (r), appearing in the definition of A mn (3.18), is given by: The appropriate initial condition for the differential equation (3.17) is identified by combining the initial condition (3.5) with the expansion (3.8) and then applying orthogonality (3.16): Assembling the differential equations (3.17) and initial condition (3.19) for m ∈ N + together yields a matrix system: where the entries of A are defined in equation (3.18) and T 0 = [T 1 (0), T 2 (0), . ..]T .The exact solution of (3.20) is expressed in terms of a matrix exponential with the nth entry of T( t) defining the time function T n ( t) in the solution expansion (3.8).

Fraction of drug released
To characterise the release process, we calculate the cumulative fraction of drug released as a function of time, M (t).This quantity is obtained by integrating the concentration flux over the outer surface of the spherical capsule and normalizing by the initial mass of drug loaded in the capsule: where n is the unit vector normal to ω directed outward from Ω.In spherical coordinates under radial symmetry, M (t) simplifies to Finally, inserting the solution expansion (3.8) gives the final form for the fraction of drug released where the nth entry of the vector U( t) = t 0 T(s) ds = (e tA − I)A −1 T 0 defines U n ( t) = t 0 T n (s) ds.

Numerical results
We now present numerical results for two distinct cases: pure diffusion (k(r) = 0) and reaction-diffusion (k(r) > 0).All results are calculated using a uniform initial concentration c 0 (r) = C 0 and the parameter values given in Table 1.We have implemented the analytical solution in MATLAB truncating all infinite series at N = 150 terms and using a simple bisection method to solve the eigenvalue equation ( 3 1 and four different choices of α and σ.Observe that the smallest value of α accurately captures the case of a homogeneous medium (e.g.D(r) = D avg ), while the largest value of α accurately captures the case of a two-layer medium (e.g. For the case of pure diffusion, drug is released differently depending on the parameter α (Figs 4, 8 and 9) and these results demonstrate the wide variety of concentration profiles using FGMs.Actually, in some circumstances, maintaining local drug concentrations within some defined therapeutic range is desirable, while in other cases a burst of drug is required.If the release is not controlled appropriately, this can lead to periods where toxic and/or sub-therapeutic concentrations are achieved.The inclusion of the spatially-varying diffusivity (2.7) provides greater control over the drug release profile (through the single tuning parameter α) while maintaining the same average material properties (2.8).In  we see that the fraction of drug released increases more rapidly when decreasing the value of α, giving rise to a quicker overall drug delivery for small values of α and a more-sustained release for large values of α.The initial burst of dose at small α may be beneficial when a rapid delivery, rather than a delayed sustained release, is desired.
Let us now consider the case of reaction diffusion.When reaction is included in the model, the capsule retains a fraction of mass that remains bound to the polymer and is never released.This yields concentration profiles (Figs 6 and 10) that decrease more rapidly than the corresponding profiles obtained for pure diffusion (Figs 4 and 9).Since the total cumulative released mass is less than the initial mass of drug loaded in the capsule (cf.equation (3.22)), the cumulative fraction of drug released asymptotes to a value less than one ( Fig 7).Comparing the drug release profiles for reaction-diffusion (Fig 7 ) with the corresponding profiles for pure diffusion (Fig 5) it is clear that (i) the presence of reaction results in decreased release times and (ii) the total released mass decreases for increasing values of α (Fig 7).Hence, the inclusion of reaction in the model provides further control of the drug release profile through the tuning parameter α.In summary, our model suggests that different functional forms for the diffusivity and reaction rate are able to provide a variety of drug release characteristics different from those provided by the limiting cases of a homogeneous system (constant coefficients) and two-layer system (stepwise coefficients).

Conclusions
Polymeric engineered materials have been exploited in a range of different application in biomedicine, aerospace and material science and can be useful in the pharmaceutical industry in the field of drug delivery.With recent advances in bioengineering, novel functionally graded materials (FGMs) have been introduced for the development of drug releasing devices and systems.They contribute to the tailoring of material for optimal drug administration including targeted release and customizability.The goal of the present study was to elucidate the potential transport mechanism and the drug kinetics behavior due to the diffusion and reaction shape-material functions, providing guidance for design micro-structure of polymer platforms and capsules development.A variety of space-dependent reactiondiffusion FGM shapes has been proposed and implemented.The hybrid analytical-numerical solution improves the understanding of the mass transfer from an FGM capsule, including in the presence of a binding reaction.
However, it is important to recognize some limitations of the present one-dimensional model.Drug dynamics in the release medium outside the capsule are ignored, so that the interactions between the capsule and the medium are represented entirely by the boundary condition at the external surface.Diffusion coefficients are also assumed to be independent of concentration.For most practical applications such assumptions are reasonable, and therefore, the model results may be helpful in the evaluation of drug kinetics and may provide new pathways for smarter delivery systems.The strategy proposed here, once calibrated, can be utilized in a predictive way to limit the number of experiments.Thus, by showing the correlation between properties of the drug kinetics and material function variables, our model can be used to determine and optimize the processing parameters to ensure a controlled drug delivery within a certain time.

Figure 3
Figure 3 displays the diffusivity function D(r) and reaction rate function k(r) for the parameter values in Table1and four different choices of α and σ.Observe that the smallest value of α accurately captures the case of a homogeneous medium (e.g.D(r) = D avg ), while the largest value of α accurately captures the case of a two-layer medium (e.g.D(r) = D max if 0 < r < R/2 and D(r) = D min if R/2 < r < R).For the case of pure diffusion, drug is released differently depending on the parameter α (Figs 4, 8 and 9) and these results demonstrate the wide variety of concentration profiles using FGMs.Actually, in some circumstances, maintaining local drug concentrations within some defined therapeutic range is desirable, while in other cases a burst of drug is required.If the release is not controlled appropriately, this can lead to periods where toxic and/or sub-therapeutic concentrations are achieved.The inclusion of the spatially-varying diffusivity (2.7) provides greater control over the drug release profile (through the single tuning parameter α) while maintaining the same average material properties (2.8).InFig 5, Figure 3 displays the diffusivity function D(r) and reaction rate function k(r) for the parameter values in Table1and four different choices of α and σ.Observe that the smallest value of α accurately captures the case of a homogeneous medium (e.g.D(r) = D avg ), while the largest value of α accurately captures the case of a two-layer medium (e.g.D(r) = D max if 0 < r < R/2 and D(r) = D min if R/2 < r < R).For the case of pure diffusion, drug is released differently depending on the parameter α (Figs 4, 8 and 9) and these results demonstrate the wide variety of concentration profiles using FGMs.Actually, in some circumstances, maintaining local drug concentrations within some defined therapeutic range is desirable, while in other cases a burst of drug is required.If the release is not controlled appropriately, this can lead to periods where toxic and/or sub-therapeutic concentrations are achieved.The inclusion of the spatially-varying diffusivity (2.7) provides greater control over the drug release profile (through the single tuning parameter α) while maintaining the same average material properties (2.8).InFig 5,

Figure 4 :Figure 5 :Figure 6 :
Figure 4: Drug concentration as a function of radius (3.8) for the pure diffusion case.Profiles are shown at four distinct times, t = 10 −2 , 10 −1 , 10 0 , 10 1 , with the black arrow indicating the direction of increasing time.

Figure 7 :
Figure 7: Cumulative fraction of drug released as a function of time (3.23) for the reaction-diffusion case.Vertical dashed lines indicate approximate release times corresponding to the time when 99% of the total released mass has been released, i.e., the time when M ( t) = 0.99 lim t→∞ M ( t).

Table 1 :
133).The in-built MATLAB functions integral, expm and fzero are used, respectively, to (i) evaluate the integrals in equations (3.18) and (3.19) (ii) compute the matrix exponential e tA and (iii) solve the nonlinear equation (2.10) for σ.Further implementation details are available in our code, which can be downloaded from https://github.com/elliotcarr/Carr2023a. Parameters of the problem.
Figure 3: (left) Diffusivity functions D(r) (2.7) used for both the pure diffusion and reaction-diffusion cases.(right) Reaction rate functions k(r) (2.11) used for the reaction-diffusion case.For each value of α, the table gives the corresponding value of σ satisfying equation (2.10).