{'fraction_non_alphanumeric': 0.04734005155847199, 'fraction_numerical': 0.01886571361612374, 'mean_word_length': 4.342863105175292, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modelling mass diffusion for a multi-layer sphere immersed in a semi-infinite medium: application to drug delivery

We present a general mechanistic model of mass diffusion for a composite sphere placed in a large ambient medium. The multi-layer problem is described by a system of diffusion equations coupled via interlayer boundary conditions such as those imposing a finite mass resistance at the external surface of the sphere. While the work is applicable to the generic problem of heat or mass transfer in a multi-layer sphere, the analysis and results are presented in the context of drug kinetics for desorbing and absorbing spherical microcapsules. We derive an analytical solution for the concentration in the sphere and in the surrounding medium that avoids any artificial truncation at a finite distance. The closed-form solution in each concentric layer is expressed in terms of a suitably-defined inverse Laplace transform that can be evaluated numerically. Concentration profiles and drug mass curves in the spherical layers and in the external environment are presented and the dependency of the solution on the mass transfer coefficient at the surface of the sphere analyzed.


Introduction
Models of mass transfer from spheres are commonly used from both a theoretical and applicative point of view. For example, studies on drug delivery from microsphere-shaped capsules or from lipid vesicles as liposomes are currently experiencing a growing interest in regards to the role played by the carrier's geometry, in its loading, stability, toxicity and, ultimately, release performance [1]. The purpose of these systems is to maintain a desired drug concentration in the blood or in the tissue for as long as possible. Among other concurrent effects, such as dissolution and possible degradation, diffusion remains the most important mechanism used to control the release rate from drug delivery systems [2,3].
On the other hand, encapsulation with multiple durable concentric layers enhances the mechanical stability and biocompatibility, protecting the sphere from the external environment and premature degradation. For some specific applications, a thin coating film is required to envelop the whole spherical structure to protect it from chemical aggression and mechanical erosion [4]. The applicative goal is to accurately predict the drug release profile from a spherical capsule and improve the overall therapeutic efficacy and safety of these drug carrier systems.
Diffusion-driven mass transfer is normally described by Fick's first and second laws [5] and transport in porous media are governed by mass diffusion and convective flow models such as Darcy and the Brinkman models [6]. So far, several exact and approximate solutions have been developed to analyze the kinetics of a dispersed solute from a polymeric matrix having a spherical shape. The mass diffusion problem is analogous to the problem of heat transfer from a sphere that has been solved by many authors in the past with a large amount of published models and approaches. A number of configurations of heat diffusion have been treated in the pioneering work of Carslaw and Jaeger [7]: the spherical matrix can have different surface conditions, with a prescribed inward or outward flux or be in contact with a well-stirred medium. The case of a composite sphere with no contact resistance has been solved with Laplace transform in analogy with that of contiguous slabs. Since the 60s, Higuchi [8] derived analytic solutions for a single sphere in a perfect sink using pseudo-state approximations, without a boundary layer effect. In the classical book of Crank [5], the diffusion in a sphere from a well-stirred medium is solved by Fourier expansion in the case of constant or time varying surface concentration and the case of a constant flux at the surface. More recent work includes empirical, semi-empirical and mechanistic diffusion models [9]. An exact solution for diffusional release of a dispersed solute from a spherical polymer matrix into both semi-infinite and finite external mediums has been developed [10,11]. Simulations with Monte-Carlo techniques have also been used, where an exponential expression for drug release is prescribed [12]. For a comprehensive review of existing mathematical models for mass transfer from polymeric microspheres and transport in tissues, the reader is referred to [13].
In a recent work, the problem of a releasing spherical composite capsule has been solved [14]: therein the external medium has been confined by a finite length (release distance), say a cut-off beyond which the concentration remains sufficiently small and constant. However, in in-vitro experiments or in in-vivo cases, micro-or nano-spheres are immersed in a bulk ambient medium of size several orders of magnitude larger than that of the sphere itself. Due to the difference of scales, this surrounding medium is considered semi-infinite.
In the present paper, a semi-infinite medium is considered around a bi-layered sphere made of an inner core and an outer shell of different drug diffusion coefficients and a rigorous analytic solution based on the Laplace transform is proposed. This approach, used for other biologicaloriented models [15], in this paper combines ideas presented previously [17,18] with some novel features. With the assumption of continuity of diffusive flux between layers, the basic idea of our solution approach is to set the diffusive flux at each of the interlayer surfaces to be equal to an unknown function of time [16,18]. This allows the multi-layer problem to be reformulated into a series of coupled single layer problems, which are then solved using the Laplace transform subject to a constraint that enforces that the solutions in each layer satisfy the second specified interlayer condition (after continuity of diffusive flux). The novelty of our approach is that we consider a model consisting of a number of finite layers (spherical shells) together with the semi-infinite outermost layer (ambient medium). The method also avoids the computation of eigenvalues or orthogonal eigenfunction expansions (as used in previous works [14,17,18]), which means that our solution expressions do not require the truncation of infinite series.
The remaining sections of this paper are organized in the following way. In the next section, we present the general mathematical model for mass transfer from (into) a multi-layer sphere into (from) a semi-infinite medium. In Section 3, we consider the special case of a two-layer sphere under the assumption of radial symmetry. This one-dimensional model is then solved using the Laplace transform in Section 4 for both the desorbing and absorbing cases. Section 5 includes extensions of the solution procedure to an arbitrary number of layers and to non-uniform initial data. Numerical results and discussion are given in Section 6.

Mass transfer from/into a composite sphere
Consider a multi-layer sphere made of an internal core or depot (Ω 0 ) enclosed by a number of durable shells (Ω i , i = 1, 2, ..., n, see Fig. 1) constituted of different materials and having specific physico-chemical characteristics. These layers are customized to allow a selective diffusion and better control the transfer rate [3]. The last shell is immersed in the external ambient medium Ω e of a large extent (relative to size of the sphere), taken as semi-infinite. In most cases, diffusion is the dominant mechanism of drug transport and, due to the composite nature of the medium, drug kinetics is hard to model and predict. Here we want to study: (i) the mass diffusion from the composite sphere into the ambient medium (outward flux, releasing/desorbing sphere) and (ii) the physically dual process of the absorption from the environment within the sphere (inward flux, absorbing sphere).
In Ω 0 , we assume that the drug dissolution occurs instantaneously compared with that of diffusion [2]. For diffusion-controlled spheres, the drug release profile is obtained by solving Fick's second law of diffusion: where c 0 is the concentration field and D 0 is the diffusion coefficient. Analogously, in the surrounding shells (Ω i , i = 1, 2, ..., n), and in the external medium (Ω e ), the following diffusion equations govern the drug transport: with c i and D i denoting the concentration and diffusivity in Ω i (possible convection or reaction terms in Ω e are considered negligible here). The above mechanistic model has been recently introduced in [14]: in such a study the semi-infinite medium has been truncated at a release distance, an artificial and, in a sense, arbitrary finite cut-off length beyond which all drug is assumed extremely small at a given time and where a perfect sink condition is applied. Similarly to the model presented in [14], in this work, we assume: (a) the composite sphere is made of homogeneous enveloping concentric layers; (b) the diffusivity is constant in each layer; (c) the process is diffusion dominated; However, in contrast to [14] we consider: (d) the release medium as semi-infinite and unstirred.
On the other hand, with respect to the classical single-layer approaches [5,7], we remove some unphysical hypotheses: (i) c 0 is kept constant at ∂Ω 0 ; (ii) c e is spatially constant in Ω e as in a well-stirred medium or in a perfect sink condition.
Modelling interfaces and external coating At the interfaces between adjacent concentric layers in the capsule, flux continuity holds: with n denoting the surface external unit normal. Additionally, due to partitioning, non-perfect contact exists at the interfaces with the constant ratio of concentration determining the partition coefficient [5]: where σ i is the drug partition coefficient between layer i and i + 1.
To prevent fast release, the sphere's outmost shell Ω n is protected with a thin coating Ω m having diffusivity D m and a small finite thickness h (Fig. 1). This coating shields and preserves the encapsulated drug from degradation and fluid convection, and guarantees a more controlled release [4]. To model the drug dynamics in Ω m , we use a simple interface condition between the outermost shell (Ω n ) and the external medium (Ω e ) that incorporates the physical properties of the coating as in [14]: where P ∝ D m /h (ms −1 ) is the coating mass transfer coefficient and σ n is related to the coating partition coefficients. Note that Eq. (2.5) includes two limit cases for P : when P = 0 the case of impermeable coating (∇c n = 0) is obtained, and if P → ∞ (coating in perfect contact) the form of the other interlayer conditions (2.4) is recovered, namely, c n = σ n c e .

The core-shell model
The two-layer sphere is by far the most representative configuration of absorbing/desorbing system and in this section we consider this special case of n = 1 (Fig. 2): an internal core (Ω 0 ) encapsulated by a single polymeric shell (Ω 1 ) surrounded by a "large" external medium (Ω e ). In other words, the core-shell system is comprised of two concentric spheres of increasing radius immersed in Ω e . Although we consider n = 1, it is fairly straightforward to extend the solution methodology presented here to any number of concentric layers as we discuss later in Section 4.2.
We assume the net drug transport occurs in the radial direction only, and therefore we consider a radially-symmetric one-dimensional model (Fig. 2) as in [14]. In this case, the general formulation Figure 2: Schematic representation of the cross-section of the radially symmetric two-layer sphere, comprising an internal core Ω 0 , the concentric layer Ω 1 and the thin coating layer Ω m (in red, zoomed on the right). This sphere together with the semi-infinite external medium Ω e , constitutes a three concentric layer system.
of Section 2 is reduced to a three-layer problem, which in 1D radial symmetry reads: where r is the radial coordinate.
The initial conditions are: where C 0 , C 1 and C e are specified constants. For a desorbing/releasing sphere we take C 0 > 0 and C 1 = C e = 0 while for an absorbing sphere we take C 0 = C 1 = 0 and C e > 0.

Solution procedure
We now present our solution approach for solving the core-shell model (3.1)-(3.10). As a first step, we normalize the variables, the parameters and the equations via the change of variables: and by redefining the nondimensional constants: where C max := max(C 0 , C 1 , C e ) and D max := max(D 0 , D 1 , D e ).

Computing the concentration
Let us consider first the solution in the internal core (Ω 0 ). The Laplace transform of (4.5) yields the following boundary value problem for c 0 (r, s) := L{c 0 (r, t)}: The general solution of (4.8) is where µ 0 (s) := s/D 0 . Now, consider: The boundary condition (4.9) requires B = 0. Using (4.12) with B = 0 and applying (4.10) allows A to be identified: In summary, the solution of (4.8)-(4.10) is given by: where Carrying out a similar process for the polymeric shell (Ω 1 ) and external medium (Ω e ), we obtain the following expressions for the Laplace transforms of c 1 (r, t) and c e (r, t): where: , (4.18) ,  Applying the inverse Laplace transform to Eqs (4.14), (4.16) and (4.17) yields the concentration in each layer:  To evaluate the solutions (4.21)-(4.23) at a given time t, the following quadrature formula [19] is used to calculate the inverse Laplace transforms: for i = 0, 1, e and j = 0, 1, where s k := z 2k−1 /t and Re{·} denotes the real part. The constants f 2k−1 and z 2k−1 , k = 1, . . . , N/2, are defined as the residues and poles, respectively, of the best (N, N ) rational approximation to the exponential function on the negative real line, computed via the Carathéodroy-Fejér method [19]. All results generated in this paper are computed by setting N = 14 [17]. The computation (4.24) requires evaluation of the unknown functions g 0 (s) and g 1 (s) at s k = z 2k−1 /t. To this aim, we solve the following linear system derived by taking Laplace transforms of the (as yet unused) interface conditions in Eqs (3.5) and (3.6): Substitution of (4.14), (4.16) and (4.17) into (4.25) and (4.26) produces a linear system of equations: where: Solving the linear system (4.27) for x yields the evaluations g 0 (s k ) and g 1 (s k ), appearing in the numerical inverse Laplace transforms (4.24), which completes the solution.

Computing the mass
The drug mass in all concentric layers as a function of time is defined as the volume integral of the concentration, which simplifies to the following formulas under the assumption of radial symmetry [14]: To evaluate the integrals, the solution expressions (4.21)-(4.23) are substituted into (4.31)-(4.32) and the resulting integrals computed numerically. Using the initial conditions (3.8)-(3.10), we have: (4.33) For the desorbing case with initial data C 0 = 1 and C 1 = C e = 0: while for the absorbing case with initial conditions C 0 = C 1 = 0 and C e = 1: Note that in the desorbing case, as described in Eqn (4.34), all the initial (finite) mass is transferred outside to the external medium (Ω e ), while in the absorbing sphere an initial (infinite) mass is initially given in Ω e and residual mass remains there for all time.

Extension of the solution procedure
Arbitrary number of concentric layers We now discuss extension of the analytical solution derived in Section 4.2 to the case of a spherical core enveloped by an arbitrary number n of concentric spheres with increasing radius, such that R 0 < R 1 < . . . < R n . In this case, the analogue of the model (3.1)-(3.10) is given by: subject to the boundary, interface and initial conditions: . . , n, c e (r, 0) = C e in (R n , ∞).
The general problem above can be solved using an identical procedure to that described for the case study in Section 3. Firstly, the problem is reformulated into a series of n + 2 single layer problems in a similar manner to that described in Section 4.1 by first defining the following n + 1 unknown mass fluxes:  The inverse Laplace transforms are computed using (4.24) with the evaluations g j (s k ) (j = 0, 1, . . . , n) satisfying an extended version of the system (4.27) with dimension (n + 1) × (n + 1), where A is a tridiagonal matrix. The entries of A, x and b are defined as follows: for j = 1, ..., n + 1, A j,j−1 = a j−1,0 (R j , s k ) for j = 2, ..., n + 1, for j = 1, ..., n, for j = 1, ..., n + 1, 0, 1, . . . , n − 1 and a n (R n , s k ) := a n,1 (R n , s k ) − σ n a e,0 (R n , s k ) − P −1 .
Non-uniform initial data So far, in both the mathematical modelling and solution procedure, we have assumed uniform initial conditions in each concentric layer, i.e., C 0 , C 1 and C 2 in Eqs (3.8)-(3.10) are constants. In this section, we revisit the two-layer sphere and expand the solution procedure outlined in Section 4 to spatially-dependent initial conditions. To explain the process, consider the model described in Section 3 with Eq. (3.8) replaced with the non-uniform initial condition: In this case, Eq. (4.8) becomes: The general solution of (5.7) is now: where the particular solution, c (p) 0 (r, s), can be derived using the method of variation of parameters [20]: Applying the boundary conditions (4.9)-(4.10) produces the modified form of (4.14) for the case of non-uniform initial data in the first layer: c 0 (r, s) = c 0 (r, s) + a 0,1 (r, s)g 0 (s), (5.10) where a 0,0 (r, s) is as defined in Eq. (4.15) and The inverse Laplace transform yields the modified solution in the first layer: where L −1 { c 0 (r, s)} can be evaluated using the quadrature formula (4.24) as follows: The evaluations g 0 (s k ) and g 1 (s k ) are computed as in Section 4.2 with the exception that (4.30) is replaced with: (5.14)

Solution verification for a homogeneous sphere
We first assess our Laplace transform solution (4.21)-(4.23) using the homogeneous analogue of the core shell-model (3.1)-(3.10): where we have a single-sphere of radius R and take f (r) = C 0 for r < R C e for r > R, with C 0 = 1 and C e = 0 for the desorbing case and C 0 = 0 and C e = 1 for the absorbing case. The exact solution of this simplified problem is well-known for an arbitrary initial function f (r) [21]. Substituting the specific forms of f (r) into this exact solution yields for the desorbing case: and for the absorbing case: where the solutions are valid for r > 0 in both cases. In Fig. 3, we compare the above classical solutions to our Laplace transform solution (4.21)-(4.23) for dimensionless values of D = 1 and R = 0.4. To solve the homogeneous model, we use the following choices of parameters in the core shell-model (3.1)-(3.10): D 0 = D 1 = D e = D, R 1 = R and P → ∞ (i.e. P −1 = 0 in Eq. (4.28), with C 0 = C 1 = 1 and C e = 0 for the desorbing case and C 1 = C 2 = 0 and C e = 1 for the absorbing case. Note that the concentration profiles are depicted on a truncated finite domain, 0 ≤ r ≤ 1, as beyond r = 1 the solution is effectively constant. Clearly, at all times shown, both the Laplace transform solution (4.21)-(4.23) (continuous lines in Fig. 3) and classical solutions (6.1)-(6.2) (markers in Fig. 3) are in excellent agreement. In summary, these results confirm numerically that our analytical solution correctly reduces to the exact solution of the homogeneous model when D 0 = D 1 = D e = D, σ 0 = σ 1 = 1 and P → ∞.
Desorbing case Among many applicative fields of releasing spheres, we focus here on layer-by-layer coated capsules, as controlled drug carriers. They have attracted significant attention for therapeutic applications and deserve special interest because of their potential for sustained release. For the desorbing case, the drug is transported from inner core, via the intermediate shell, to the release medium: each layer receives mass from the layer beneath it and transfers it to the layer above it, in a cascading sequence until the drug is completely released from the capsule. The coating mass transfer coefficient P constitutes the distinctive parameter that controls the flux exiting the capsule. Fig. 4 shows the concentration profiles in the case of two different values of P : P → ∞ (uncoated sphere, Fig. 4a) and P = 5 · 10 −8 (coated sphere, Fig. 4b). Concentration is decreasing inside each layer and is discontinuous at the interlayer interfaces for finite P , with the mass flux continuity preserved (Fig. 4ab). Excellent agreement is achieved when comparing our results to those presented in [14]. Due to the sink boundary condition (3.7), all mass eventually accumulates in the external release medium. In other words, due to the condition (3.7), all drug mass is released into the environment after a sufficiently long time and the total mass is preserved. The drug mass monotonically decreases in the core (layer 0), while at the same time increasing up to some peak before decaying asymptotically in the hydrogel layer (layer 1) (Fig. 4cd, with the mass normalized by its initial value 4 3 πR 3 0 C max ). In the release medium, the mass progressively increases at a rate depending on the diffusive properties of the two-layer materials. The simulation indicates that the time and the size of the mass peak in the hydrogel layer (layer 1) is related to the releasing properties of the core, on the one hand, and to the diffusivity of the release medium, on the other hand, together with the mass resistance of the coating. The thin hydrogel layer retains a negligible mass due to its thickness, and the core is completely emptied after roughly 10 hours, in the case of P → ∞. After that time, all the mass is transmitted to the external medium. A much more sustained release occurs in the case of a coating having a finite and small mass transfer coefficient (P = 5 · 10 −8 ). After 10 hours, a substantial amount of drug remains in the core and hydrogel layers, with the core not completely empty until approximately 22.5 hours.
We now investigate the sensitivity of the solution to the value of coating mass transfer coefficient P as an effective rate-controlling parameter. It turns out that for the above parameters (6.3), the sensitive values of P are in the range: 10 −8 ≤ P ≤ 10 −3 . For P = 10 −8 , the coating almost acts as an impermeable barrier with a very small transfer rate from the capsule to the external medium evident in Fig. 5b. Interestingly, decreasing P by a factor of 5 from the case P = 5 · 10 −8 considered earlier (Figs. 4bd) has a huge effect on the release rate: after 30 hours roughly 20% of the mass still remains in the capsule (core and hydrogel layers, Fig. 5bd). Setting P = 10 −3 produces results that are indistinguishable from P → ∞ as there is no observable difference between Figs 4ac and 5ac. In this case, the capsule surface is in perfect contact with the external ambient medium as evident by the continuity in concentration (Fig. 5a).
Absorbing case The process of drug absorption from a saturated solution determines, in part, its bioavailability and is the basis of in-vitro experiments. Drug kinetics from environment into an absorbing sphere is very similar to the desorbing case, except that the initial mass is present in the external semiinfinite medium and the drug transport direction is reversed. Fig. 6 (with mass normalized by 4 3 πR 3 0 C max ) is the counterpart of Fig. 4 and shows how, and to what extent, the drug diffuses into the two-layer sphere. In the case of a coated sphere, the diffusion rate is lower and the drug reaches saturation after a longer period of time. In contrast to the desorbing case, the polymeric shell (layer 1) fills up to a maximum concentration and receives diffused mass from the external source, taken as a large reservoir, and transfers it to the inner core. The normalized saturation mass in layer 1 (Fig. 6cd) depends only on the geometrical configuration (see Eq. (4.35)):

Conclusions
A better understanding of the mass transfer from a drug carrier or a vehicle in a living tissue for therapeutic purposes constitutes an important challenge in medicine nowadays. Mathematical modelling helps in predicting the drug release rates and diffusion behavior from these delivery systems, thereby reducing the number of experiments needed.
In the current work, inspired by the above biomedical application, we have presented a mathematical model and developed an analytical technique to study the diffusion-controlled mass desorption-absorption systems in microspheres. We have focussed on drug release from twolayer spherical capsules, consisting of an inner core and an outer shell protected by a thin coating, immersed in an external semi-infinite medium. This type of geometry, introducing additional de-sign parameters to the formulation (i.e. relative sizes and relative drug diffusion and partition coefficients of inner and outer structures), enriches the possibilities in terms of pharmacokinetics, controlled release rate and, ultimately, delivery performance. As in many biological systems, the model contains a number of parameters, subject sometimes to high variability and uncertainty, that need to be identified before it can be used in a predictive way to provide the drug kinetics. Once the parameters are identified, the proposed methodology provides a simple tool that can be used to quantitatively characterize the drug diffusion, improve the technological performance and optimize the release rate for therapeutic purposes. By virtue of the one-to-one analogy of mass diffusion and heat conduction problems, the presented approach can be successfully applied to the similar model of heat transfer from/in a sphere immersed in a large ambient medium.