Derivation of an Analytical Solution to a Reaction-Diffusion Model for Autocatalytic Degradation and Erosion in Polymer Microspheres

A mathematical reaction-diffusion model is defined to describe the gradual decomposition of polymer microspheres composed of poly(D,L-lactic-co-glycolic acid) (PLGA) that are used for pharmaceutical drug delivery over extended periods of time. The partial differential equation (PDE) model treats simultaneous first-order generation due to chemical reaction and diffusion of reaction products in spherical geometry to capture the microsphere-size-dependent effects of autocatalysis on PLGA erosion that occurs when the microspheres are exposed to aqueous media such as biological fluids. The model is solved analytically for the concentration of the autocatalytic carboxylic acid end groups of the polymer chains that comprise the microspheres as a function of radial position and time. The analytical solution for the reaction and transport of the autocatalytic chemical species is useful for predicting the conditions under which drug release from PLGA microspheres transitions from diffusion-controlled to erosion-controlled release, for understanding the dynamic coupling between the PLGA degradation and erosion mechanisms, and for designing drug release particles. The model is the first to provide an analytical prediction for the dynamics and spatial heterogeneities of PLGA degradation and erosion within a spherical particle. The analytical solution is applicable to other spherical systems with simultaneous diffusive transport and first-order generation by reaction.

Introduction applied directly to PLGA because the reaction kinetics in those models were independent of the drug diffusion dynamics. The reaction and diffusion processes for PLGA microspheres are coupled and should not be treated independently.
A reaction-diffusion model has been proposed previously [24] to treat autocatalytic degradation and erosion in the polymer plates and cylinders composed of chemically similar poly (lactic acid) (PLA), and the numerical predictions were compared to experimental data from a study on the degradation of PLA plates of different thicknesses [25]. While the model in [24] can be extended to other geometries using finite element analysis to solve the equations numerically, the influence of the relationship between diffusion and reaction parameters on the zones where diffusion and degradation have strong or weak influences was only presented for onedimensional planar and cylindrical geometries. We aim to derive an analytical solution to a reaction-diffusion model specifically for the case of spherical geometry that can assess the relative dominance of the reaction and diffusion phenomena through a single dimensionless parameter.
Reaction-diffusion equations have commonly been used to model spherical catalyst pellets that experience simultaneous reaction and diffusion, and analytical solutions to the equations are available [26,27]. Despite the similarities in the reaction-diffusion equations for the two systems, the catalyst pellets and PLGA microspheres differ in ways that are critical for their mathematical treatment: (1) the reaction term in PLGA is a generation term instead of a consumption term and (2) the directions of diffusion are reversed for PLGA microspheres and catalyst pellets. In catalyst pellets, the reactant diffuses into the sphere where it is consumed by a reaction. In contrast, the autocatalytic carboxylic acid end groups in PLGA microspheres are distributed throughout the sphere where more are generated by the degradation reaction, and some fraction of the autocatalyst diffuses out of the sphere. These differences between the physics of reaction and diffusion in spherical catalyst pellets and PLGA microspheres must be carefully accounted for during the analogous mathematical analyses of the reaction-diffusion equations used to model the systems.
Here, we derive an analytical expression for the transient, radial concentration of the autocatalytic carboxylic acid end groups of PLGA by simultaneously treating degradation and erosion of the polymer. Quantifying the transient, spatial distribution of the autocatalyst within the polymer is important for understanding how different conditions may contribute to accelerating hydrolysis in the interior of large microspheres and for preventing adverse effects of the acidic conditions within a microsphere on the drug stability or bioactivity [28]. Additionally, such a mathematical treatment provides insights into how the coupling between autocatalysis and drug diffusion may trigger transitions between diffusion-controlled and erosion-controlled drug release regimes.

Mathematical model for the autocatalyst concentration
The conservation equation for a chemical species subject to reaction and diffusion within a radially symmetric sphere is where c(r, t) is the concentration, 0 r 1 is the normalized radial position defined as r ¼r=R,r is the the radial distance from the center of the sphere, R is the radius of the sphere, t ! 0 is time, D is the diffusion coefficient, and R V is the net rate of generation of species per volume. We assume that the microsphere volume is constant for bulk-eroding PLGA, neglecting any effects of possible microsphere swelling. The chemical species of interest is the autocatalytic carboxylic acid end groups of the polymer chains or autocatalyst. PLGA degradation is often treated with pseudo-first-order kinetics [13,18], where the autocatalyst undergoes first-order growth while the concentrations of water, c H 2 O and ester bonds in the polymer, c E , are assumed to remain constant, giving where the rate constant k for random ester bond hydrolysis incorporates the constant concentrations c H 2 O and c E (the applicability of the assumption of constant c E is checked in the Results and discussion section). PLGA erosion is treated by assuming that all carboxylic acid end groups have a uniform constant diffusion coefficient D independent of the length of the polymer chain to which they are attached; with this assumption, diffusion of degradation products can be included in the analytical expression for autocatalyst concentration. Substituting α = D/R 2 and (Eq 2) into (Eq 1) yields for k > 0, and α > 0. The boundary conditions are @cð0; tÞ @r ¼ 0; t ! 0 ð4Þ and cð1; tÞ ¼ c r 1 ; t ! 0; ð5Þ and the initial concentration distribution is where c r 1 ! 0 and c t 0 > 0. Note that c r 1 c t 0 for flux toward the exterior of the sphere.

Analytical solution for the autocatalyst concentration
By substituting v(r, t) = rc(r, t), the partial differential equation (PDE) (Eq 3) and its initial and boundary conditions are transformed to a linear, nonhomogeneous second-order PDE with a source term and nonhomogeneous Dirichlet boundary conditions: with boundary conditions The PDE for v(r, t) can be solved by superposition of the solutions to two related problems [29]: (i) the complementary steady-state boundary value problem, with boundary conditions u(0) = 0 and u(1) = c r 1 , where u(r) is the equilibrium steady-state distribution, and (ii) the PDE for the function w(r, t) v(r, t) − u(r), with homogeneous boundary conditions and initial condition The steady-state solution, u(r), and the solution to (Eq 12) for w(r, t) by the method of eigenfunction expansion are derived below and then combined to obtain the solutions for v(r, t) and c(r, t).
Steady-state solution for u(r). The steady-state boundary value problem (Eq 11) can be rewritten as where the Thiele modulus, F, for this first-order reaction-diffusion system is The Thiele modulus characterizes the relative importance of diffusion and reaction and is defined as the ratio of the characteristic times for diffusion (1/α) and reaction in the absence of mass transfer limitations (1/k). The solution is of the form [30] u ¼ A cos ðFrÞ þ B sin ðFrÞ: ð17Þ The boundary condition u(0) = 0 for finite values of c requires that A = 0. At the surface, so B = c r 1 / sinF. The solution to the steady-state ordinary differential equation (ODE) with nonhomogeneous boundary conditions is Method of eigenfunction expansion solution for w(r, t). The linear, nonhomogeneous PDE with homogeneous boundary conditions given by (Eq 12) can be solved using the method of eigenfunction expansion [29], which consists of expanding the unknown solution w(r, t) in a series of the eigenfunctions for the related homogeneous problem: where a n (t) are the time-dependent, generalized Fourier coefficients and ϕ n (r) are the eigenfunctions of the homogeneous PDE for diffusion without a source term with homogeneous Dirichlet boundary conditions, The eigenfunction expansion of the source term is where b n (t) = ka n (t) since the source term is first-order in the linearized concentration.
The Fourier sine series can be differentiated term by term since w(r, t) and sin(nπr) satisfy the same boundary conditions [29]. Inserting the eigenfunction expansions for w(r, t) from (Eq 20) and the source term from (Eq 22) into the PDE in (Eq 12) yields For each n = 1, 2, . . ., The solution to (Eq 24) is a n ðtÞ ¼ a n ð0Þ exp ðÀðn 2 p 2 À F 2 ÞatÞ; ð25Þ where the a n (0) can be derived by multiplying (Eq 20) by ϕ m for t = 0, considering the orthogonality of the eigenfunctions, and integrating over the spatial domain to give a n ð0Þ ¼ In the specific case of uniform initial distribution, c t 0 (r) = c t 0 and a n ð0Þ ¼ Substituting the expressions for ϕ n from (Eq 21), a n from (Eq 25), and a n (0) from (Eq 26) into (Eq 20) yields the solution to the nonhomogeneous PDE with homogeneous boundary conditions (Eq 12): Solution for v(r, t). Substituting the expressions for u from (Eq 19) and w from (Eq 28) into v = u+w yields the solution to the nonhomogeneous PDE with nonhomogeneous boundary conditions (Eq 7): Solution for c(r, t). The concentration of the autocatalyst in radial coordinates is The indeterminate form at r = 0 is resolved by where x denotes a constant. The autocatalyst concentration at r = 0 is and at 0 < r < 1 is 2rc t 0 ðrÞ sin ðnprÞdr þ 2c r 1 npðÀ1Þ n n 2 p 2 À F 2 Â exp ðÀðn 2 p 2 À F 2 ÞatÞ sin ðnprÞ r : With uniform initial distribution, (Eq 32) and (Eq 33) become Numerical solution for the autocatalyst concentration The partial differential equation for the autocatalyst concentration Eqs (3)-(6) can also be solved by numerical methods. The method of lines [31,32] reduces the PDE to a system of ODEs by discretizing the radial dimension onto a finite grid with equal spacing Δr and coordinates r i = iΔr for i = 0, 1, . . ., N. Using the boundary conditions and the second-order centered finite difference approximations for spherical coordinates [33] to approximate the spatial derivatives, the semi-discrete ODE for the concentration at each grid point, c i (t), is > > > > > > < > > > > > > : The system of ODEs was solved using the function ode45 in MATLAB. The numerical solution gives suitable accuracy when compared with the analytical solution.

Transient autocatalyst profiles
Our analytical expressions for the transient radial concentration of the autocatalyst Eqs (32)-(35) quantify the importance of the Thiele modulus, F, for determining whether PLGA degradation and erosion are enhanced by accumulation or diminished by diffusion. The results presented here are all for the case of uniform initial distribution c t 0 (r) = c t 0 given by Eqs (34)-(35).
For small values of F (Fig 1a-1b), diffusion dominates the conservation equation, so any amount of the autocatalyst generated by the reaction diffuses away quickly and a steady state can be reached. Such a microsphere experiences homogeneous erosion without substantiallycatalyzed degradation (Fig 2a). Small microspheres, fast diffusion, or slow reaction can give small values of F. For intermediate values of F on the order of π (Fig 1c-1d), the autocatalyst accumulates in the center of the microsphere forming a hollow interior (Fig 2b). For large values of F (Fig 1e), the first-order reaction generation dominates the conservation equation, so the autocatalyst accumulates in the interior faster than it can diffuse from the sphere. Autocatalytic degradation is severe throughout the microsphere (Fig 2c), and the entire particle may erode rapidly. The predicted autocatalyst concentration profiles (Fig 1) and the illustrations of the spatial distribution of the autocatalyst in microspheres of different sizes where small sizes correspond to small values of F (Fig 2) are consistent with experimental visual evidence of the microspheres-size-dependent spatial distribution of the acidic microclimate within degrading PLGA microspheres [34].

Transition between regimes
To assess the bounds on the concentration profiles in different ranges of F, the maximum values of the concentration profiles, which occur at the microsphere center r = 0, are compared (Figs 3 and 4). The upper and lower bounds are the reaction-dominant limit and the diffusiondominant limit, respectively. Considering only the degradation reaction without diffusion, the autocatalyst concentration is described by which has the solution where c rxn is referred to as the reaction-dominant limit. This limit is the commonly treated case for autocatalytic hydrolysis kinetics [13], where the concentration at each position is independent of its neighbors. As F ! 1, c(r, t) ! c rxn (t) (Fig 3), so values in the range F ! π are in the erosion-controlled regime.   In a PLGA microsphere, the degradation reaction cannot continue infinitely because the hydrolysis reaction only proceeds as long as ester bonds can be cleaved in the polymer. The maximum time for the reaction occurs when all of the ester bonds have been converted to monomers, so c rxn (r, t max ) = c t 0 + c E and where M n 0 is the initial number-average molecular weight and M 1 is the average weight of a monomer. As t ! t max the assumption of constant ester bond concentration is violated. At t = 1/k, t/t max % 1-2% for PLGA 50:50 with M 1 = 65.05 Da, M n 0 = 10-1000 kDa, and k = 0.012-0.08 day −1 [7], thus the constant ester concentration assumption is valid within the characteristic time for the reaction (kt = 1 gives t % 12-80 days). Considering only diffusion without the degradation reaction, the autocatalyst concentration is described by with boundary and initial conditions given by Eqs (5)- (6), which has the solution [33] for r = 0 and for 0 < r < 1 where c diffn is referred to as the diffusion-dominant limit. As F ! 0, c(r, t) ! c diffn (r, t) (Fig 4), so values in the range F < π are in the diffusion-controlled regime. For intermediate values of F, the autocatalyst concentration is bounded by the diffusion-and reaction-dominant limits; the concentration peaks due to early accumulation and enhanced degradation in the center and later decays due to diffusion. F = π is the tipping point between the diffusion-and erosion-controlled regimes characterized by exponential decay and exponential growth of the autocatalyst, respectively.

Steady-state autocatalyst profiles
The autocatalyst concentration, c(r, t), exhibits exponential growth due to the first-order hydrolysis reaction and exponential decay due to diffusion. The steady-state limit is approached as the rates of diffusion and generation by reaction offset each other or when the reacting species has diffused out of the system completely. For the concentration to reach a steady state, the time derivative of c(r, t) must be zero; either the exponential term must (i) be constant or (ii) approach 0 as t ! 1. Therefore, In the most restrictive case of n = 1, F π.
The constant surface boundary condition c r 1 determines the magnitude of the diffusion driving force and indicates the concentration of species in the medium assuming no mass transfer limitations at the surface; this is consistent with a buffered medium where the autocatalyst released from the microsphere is perfectly absorbed by the medium. If c r 1 = 0, values of F < π have trivial steady-state solutions c(r, t)/c t 0 ! 0 as t ! 1 where the species diffuses completely out of the sphere preventing runaway of the first-order, autocatalytic reaction, while F = π has a stable nontrivial solution c(r, t)/c t 0 ! 2 as t ! 1, balancing the contributions from the reaction and diffusion. If c r 1 > 0, the steady-state solutions for F = mπ are nontrivial (Fig 5), and neither the steady-state linearized solution, u, nor the analytical solution, c, are defined for F = mπ, m = 1, 2, . . .. The transient profiles for F = mπ can be approximated by the numerical solution even though steady-state solutions do not exist in these cases. The case of c r 1 > 0 is interesting from a mathematical perspective but may only be physically appropriate for media that have high concentrations of strong acid that augment the carboxylic acid end groups in catalyzing the degradation. The boundary condition c r 1 = 0 is valid assuming no autocatalyst in a weakly acidic buffered medium and is used for the results shown in Figs 1, 3, and 4.

Conclusions
An analytical expression was derived for the transient, radial concentration of a species undergoing simultaneous diffusion and first-order reaction generation with constant, but not necessarily zero, surface boundary concentration. The expression differs from the common reaction-diffusion case treated in the spherical catalyst literature as we treat a generation rather than consumption reaction term. To our knowledge, this is the first application of such an analytical expression for the reaction-diffusion equation to PLGA microspheres to model degradation and erosion of the polymer. Treating the diffusive transport of the autocatalytic species to capture spatial heterogeneities during degradation is a unique contribution of this work.
A limitation of the analytical expression for the autocatalyst concentration is the assumption that the diffusion coefficient is independent of the lengths of the polymer chains to which the autocatalytic carboxylic acid end groups are attached. In reality only a fraction of the polymer chains are small enough to be water-soluble and able to diffuse through the aqueous pores. To account for the overestimation of autocatalyst mobility in the assumption needed to simplify Analytical Solution to Polymer Microsphere Reaction-Diffusion Model the mathematical formulation, a smaller value for the diffusion coefficient than the value for just the soluble chains should be used to represent the average diffusivity of soluble and insoluble polymer chains. A more detailed model distinguishing between the diffusivities of soluble and insoluble autocatalyst populations would require a numerical solution to the PDE for each population.
The analytical expression for autocatalyst concentration indicates that the Thiele modulus is a key parameter for predicting the transition between the diffusion-and erosion-controlled release regimes. With the coupling between reaction and diffusion of the autocatalyst treated by this model, size-dependent effects on autocatalysis can be explored and incorporated into detailed predictive models for drug delivery from autocatalytic PLGA microspheres, which could ultimately contribute to the in silico optimal design of controlled drug release particles.