Swelling and Drug Release in Polymers through the Theory of Poisson–Kac Stochastic Processes

Experiments on swelling and solute transport in polymeric systems clearly indicate that the classical parabolic models fail to predict typical non-Fickian features of sorption kinetics. The formulation of moving-boundary transport models for solvent penetration and drug release in swelling polymeric systems is addressed hereby employing the theory of Poisson–Kac stochastic processes possessing finite propagation velocity. The hyperbolic continuous equations deriving from Poisson–Kac processes are extended to include the description of the temporal evolution of both the Glass–Gel and the Gel–Solvent interfaces. The influence of polymer relaxation time on sorption curves and drug release kinetics is addressed in detail.


Introduction
Many fluids and specifically polymeric liquids, suspensions, and gels possess viscoelastic properties, which, in the simplest case, can be characterized through a single relaxation time [1], leading to constitutive equations with memory for the stress tensor as a function of the deformation tensor. Since momentum transport is affected by changes in the constitutive equations, it is natural to assume that the memory effects in viscoelastic materials impact on the transport properties and the associated phenomena, such as glassy-rubbery transition, swelling, and the release of a solute dispersed in the material. Indeed, experiments on swelling and solute transport in polymeric systems indicate that the classical parabolic models (in which the mass flux is proportional to the concentration gradient) fail to predict typical and distinguishing features, such as Case II diffusion and concentration overshoot in sorption experiments [2,3]. A systematic analysis of non-Fickian sorption kinetics in polymer films can be found in Sanopoulou and Petropoulos, 2001 [4], where a list, from 1 to 6, of characteristic deviations from Fickian sorption behaviour is presented.
From the experimental point of view, the effect of polymer relaxation on sorption curves can be investigated in many different ways. A couple of examples among the many present in the literature are discussed here. Focusing on hydrogels, stress relaxation experiments can be carried out to determine the time scale of macromolecular adjustments during the swelling process. Peppas and Brazel [5] measured the characteristic time for stress relaxation in a swollen polymer when subjected to a mechanical stress by the Instron and then related this time to the time required for a dry polymer to adjust to being placed in a solvent, assuming one-dimensional transport. These authors observed that the stress relaxation time constants at 37°C for redox-and thermally initiated free radical solution polymerized samples of P(HEMA-co-MMA) crosslinked by EGDMA and glutaraldehyde-crosslinked PVA samples changes significantly with the amount of crosslinking agent and PVA molecular weight. Therefore, it is possible to design different materials and different experiments to evaluate the contribution of polymer relaxation to the experimental sorption curves. Focusing on organic vapor sorption in polymer film, Sanopoulou and Petropoulos [4] proposed different diagnostic criteria for the physical origin of non-Fickian kinetic features, some strictly related to the presence of viscous relaxation phenomena. Interval sorption as well as integral sorption experiments on films with different thicknesses can be useful to identify the role and contribution of polymer relaxation. By considering that the Deborah number is inversely proportional to the film thickness squared, these authors analyzed the absorption kinetics of CA-Acetone systems and observed significant differences between sorption curves when changing the film thickness from 42 to 112 µm, thus decreasing the Deborah number by a factor of 10.
Form the theoretical point of view, a first coherent approach in this direction was undertaken by Cattaneo [6,7] in order to solve some paradoxes in the parabolic theory of heat transfer. The original articles by Cattaneo deal with the heat transfer (so that the final result of this investigation is usually referred to as the Cattaneo heat equation), but the extension to mass transport is straightforward. Cattaneo's model is indeed a hyperbolic model in which the diffusive flux satisfies a constitutive equation identical to the constitutive equations used in linear viscoelasticity.
As Fickian diffusion corresponds microscopically to a stochastic motion described in terms of Wiener processes, an important conceptual issue is the determination of the kinematic equations underlying the Cattaneo equation. In other words, which stochastic model for particle transport determines macroscopically the concentration field representing the solution of the Cattaneo equation. The ultimate answer to this fundamental question has been given by Mark Kac, at least for one-dimensional spatial problems [8], and his approach provides not only a simple kinematic model associated with the occurrence of memory effects in transport equations, but also an alternative, and physically significant way of describing transport processes via the concept of partial densities (or partial concentrations). Both the Cattaneo hyperbolic transport theory and the Kac stochastic model played a fundamental role in developing more refined theories for non-equilibrium processes, and their stochastic characterization. Cattaneo's model represented the cornerstone in the development of the Extended Thermodynamics proposed by Müller and Ruggeri [9,10] to generalize the classical theory of irreversible processes, and the Kac model is the starting point in the development of a wide class of stochastic processes possessing Markovian transitions, referred to as Generalized Poisson-Kac processes [11][12][13].
The variable surface concentration model proposed by Long and Richman as well as the diffusion-relaxation model by Berens and Hopfenberg (see [14] for a detailed review of these two classical models) represent the first attempts to describe the sorption process in glassy polymers as a linear superposition of phenomenologically independent contributions from Fickian diffusion and polymeric relaxation, the latter influencing the sorptive capacity of the polymer. The extension to more than one distinct structural relaxation process with different relaxation time has been also proposed [15,16] within the framework of a purely Fickian (diffusional) solvent transport with fixed boundaries. A significant step ahead in the field is represented by the work by Cohen and White [17], in which the major effect of a diffusing penetrant on the polymer entanglement network is taken to be the inducement of a differential viscoelastic stress. This couples diffusive and mechanical processes through a viscoelastic response, where the strain depends upon the amount of penetrant. Subsequently, many authors, e.g., Grassi et al. [18,19] first, and Wu and Brasel [20] later, following Camera-Roda and Sarti [21,22] and Cohen [17], considered a solvent flux, different from the simple Fick's law, to account for polymer relaxation in mass transport, thus arriving at a hyperbolic formulation of the transport problem. However, in these works, the effect of the moving boundaries-namely, the Glass-Gel and the Gel-Solvent interfaces-has been neglected as well as the convective contribution to mass transport induced by the presence of a swelling velocity. The "swelling flux" has been accounted for in the model proposed by Lamberti et al. [23,24], but only very low values of the Deborah number have been considered, thus disregarding the effect of polymer relaxation time. In this article, the Poisson-Kac (PK) stochastic approach is adopted and extended to describe and investigate the effect of polymer relaxation time on solvent penetration in glassy polymers. To this end, an extra convective term is included in the classical partial wave transport equations to account for the contribution of the swelling point-wise velocity. Indeed, the model accounts for the three contributions to solvent transport-namely, the diffusive flux, the swelling convective flux, and the relaxation flux. The movements of the two fronts, the Gel-Solvent and the Glass-Gel interfaces, are also accounted for and described in terms of PK partial waves. This enables (I) to study the influence of the Deborah number (from low to high values) on the sorption curves and (ii) to compare the results of the hyperbolic and the classical parabolic transport schemes, accounting for both swelling and moving boundaries. The original PK model has been also generalized to define a broader class of processes possessing Markovian transitions, where the characteristic velocity and the transition rate are continuous functions of the overall density [13].
This extension permits the investigation of the effect of both a polymer relaxation time and a solvent diffusivity that are exponential functions of the local solvent concentration, as established by simplified versions of the free volume theory [25][26][27]. In particular, the introduction of a non-constant polymer relaxation time leads to interesting and unexpected behaviours of the sorption curves, due to the presence of a nonlinear convective term that facilitates solvent penetration and speeds up the outward movement of the Gel-Solvent interface.
In the paper, we also investigate the influence of the polymer relaxation time on the release kinetics of a drug initially loaded in the thin dry film. Release kinetics are strongly influenced by the polymer relaxation time and, unexpectedly, are not necessarily slowed down for intermediate values of the Deborah number.

Poisson-Kac Processes
The counterpart of a linear viscoelastic constitutive equation in mass transport is represented by the Cattaneo equation. Let c(x, t) be the solvent concentration in a polymersolvent solution and J c (x, t) its diffusive flux. We refer here to the concept of "diffusive flux" as the flux associated with pure random molecular motion in the absence of any external macroscopic drift. Below, one-dimensional spatial problems are considered for the sake of simplicity.
In the absence of sources, the balance equation for c(x, t) reads: In the classical Fickian transport theory, the constitutive equation for the flux is: In Cattaneo theory, the constitutive equation for the flux J c (x, t) takes the form: where D is a diffusivity and t r is the characteristic relaxation time. For t r = 0, Equation (3) reduces to the classical Fickian constitutive equation, in which the flux is proportional to the concentration gradient. The introduction of a constitutive equation with memory, such as Equation (3), determines major and significant changes in the nature of the diffusive propagation. To begin with, from Equations (1)-(3), it follows that the evolution equation for the concentration c(x, t) attains the form: In other words, it attains a hyperbolic character due to the second-order derivatives with respect to time. This corresponds to the evolution of solvent concentration in the form of waves with dispersion, possessing the characteristic propagation velocity √ D/t r . For any problem in mass transport involving moving molecules and particles, it is of the highest conceptual and practical interest to determine the structure of the microscopic particle kinematics that originate the emergent macroscopic behavior for the concentration field. Brownian motion, expressed with regard to the long-term properties by Wiener processes, represents the kinematics of Fickian diffusion. For the hyperbolic transport Equation (4) , M. Kac showed that the one-dimensional Cattaneo transport equation is originated from the microscopic particle motion described by the following kinematic equation [8]: where b 0 > 0 is a constant possessing the dimension of a velocity and χ(t, λ) is a Poisson counting process characterized by the transition rate λ > 0. The Poisson process χ(t, λ) attains values 0, 1, 2, . . . spanning the natural numbers. Meanwhile, in the classical applications of probability theory, Prob[χ(0, λ) = 0] = 1, with reference to Equation (5) it is more convenient to consider the following initial conditions Prob[χ(0, λ) = 0] = Prob[χ(0, λ) = 1] = 1/2, Prob[χ(0, λ) > 1] = 0, in order to avoid any initial biasing effects. From Equation (5), it follows that the trajectories associated with this kinematic law correspond to the continuous union of straight segments, with slopes of either +b 0 or −b 0 . The transition time θ between two consecutive velocity switches is controlled by the parity of (−1) χ(t,λ) is a random variable characterized by an exponential probability density function p θ (θ) = λ e −λ θ . The process X(t), the realization of which is x(t) defined by the kinematics Equation (5), is not Markovian, and this feature determines the occurrence of memory effects in the flux of constitutive equations; Equation (3). Nevertheless, it is still possible to construct a Markovian embedding for X(t) by considering the joint stochastic process (X(t), S(t)), where S(t) = (−1) χ(t,λ) is the parity function attaining values ±1, and determining the velocity direction. Since S(t) is a binary process attaining only two values, it is possible to introduce the two partial probability densities p ± (x, t) where: Additionally, the application of the Chapman-Kolmogorov equation to the embedding process (X(t), S(t)) provides the following evolution equations for p ± (x, t): These correspond to a system of two first-order linear partial differential equations. Given the partial densities p ± (x, t), also referred to as partial waves, the overall probability density function for X(t) is expressed by: The associated probability flux J p (x, t) takes the form: A probabilistic terminology is here adopted for describing the statistical properties of an ensemble of particles. That is to say, once applied to a mass transport problem, the mass/molar concentration c(x, t) is related to the density function p(x, t), by a proportionality relation, c(x, t) = K p(x, t), where the constant K admits the dimension of a mass/number-of-moles, depending on the physical meaning of c(x, t).
The use of the partial densities provides a very clear statistical understanding on the evolution of particle concentrations associated with the kinematics (5). The evolution of the overall density p(x, t) can be decomposed into the superposition of two partial waves p ± (x, t) propagating at constant velocity in the two opposite directions along the x-axis, and mutually recombining at the rate λ which defines the statistics of the transition times of the Poisson process. For this reason, the stochastic process (5) can be defined as the elementary Poisson-Kac process.
By summing Equations (7) and (8), the equation for the overall density p(x, t) follows thus: Meanwhile, the constitutive equation for the probability flux J p (x, t), defined by Equation (10), is easily obtained by taking the difference between the two Equations (7) and (8) multiplied by b 0 : Equation (12) coincides with the Cattaneo constitutive Equation (3) upon the identification of the relaxation time t r and the diffusivity D as: The original Poisson-Kac model can be generalized to define a broader class of processes possessing Markovian transitions, and referred to as Generalized Poisson-Kac processes [11,12]. For the scope of the present article, it is interesting to consider a nonlinear expansion of the model where both the characteristic velocity b 0 and the transition rate λ are continuous functions of the overall density p(x, t): In this case, the balance equations for the partial densities p ± (x, t) are altogether similar to Equations (7) and (8): The equation for p(x, t) follows by summing Equations (15) and (16): Observe that, due to the nonlinearities, the density flux is given by: where ψ p (x, t) is an auxiliary field. The evolution equation for the auxiliary field ψ p (x, t) is obtained by taking the difference between the two Equations (15) and (16) multiplied by b 0 : From this, the constitutive equation for the flux is recovered by enforcing definition (18). Equation (19) implies that either the relaxation time t r or the diffusivity D are functions of the space and time coordinates through their dependence on the overall density p(x, t), namely: However, it is important to observe that, in the Kac the equation for p(x, t) is, as expected, a parabolic transport equation: However, a new highly nonlinear convective term appears (the second term in the right-hand side of Equation (21)), that has no counterpart in the corresponding parabolic formulation of a transport model, even accounting for a non-constant concentrationdependent diffusion coefficient.

Parabolic Transport Model for Case II Diffusion
Swelling results from the solvent penetration into the polymer. The solvent enhances the mobility of polymer chains by converting the glassy matrix into a swollen, rubbery material. Sorption in glassy polymers, usually referred to as Case II diffusion, is typically characterized by two moving fronts: (i) a sharp interface between unpenetrated glass and swollen polymer (Glassy-Gel Interface GGI) that propagates inwards into the film and (ii) a Gel-Solvent Interface (GSI) that, in the absence of dissolution, moves outwards and progressively increases the gel layer thickness. A schematic representation of solvent penetration in a glassy film is depicted in Figure 1. The classical 1-d transport equations, not accounting for the relaxation time of polymer, describe the solvent penetration and the temporal evolution of the moving interfaces GGI and and GSI in terms of parabolic partial differential equations for the solvent volume fraction φ(x, t) in the gel layer: where D is the solvent diffusivity and v sw is the point-wise swelling velocity, assumed to be equal (and opposite in sign) to the solvent diffusive (volumetric) flux [28][29][30][31][32][33]: The solvent and polymer are assumed to be incompressible and to mix with no volume change [31].
At the Gel-Solvent interface GSI(t), solvent/polymer thermodynamic equilibrium φ = φ eq is assumed, and the temporal evolution of GSI(t) is described by the Stefan equation [31]: On the Glass-Gel front GGI(t), a threshold concentration to initiate swelling φ = φ G > φ 0 is assumed for the solvent [34], with φ 0 being the initial solvent volume fractions in the dry film. Correspondingly, the temporal evolution of GGI(t) reads as: When GGI(t) reaches the symmetry axis x = 0, the glassy phase disappears and the zero-flux boundary condition applies: J(0, t) = 0.
The initial conditions for φ(x, t), GGI(t), and GSI(t) are GGI(0) = L 0 − , GSI(0) = L 0 , and φ(x, 0) = φ G , with 10 −4 L 0 -i.e., it is assumed that a gel layer of infinitesimal thickness is already formed and, consistently, that the solvent volume fraction φ(x, 0) for L 0 − < x < L 0 has attained the Glass-Gel threshold value φ G . Solvent diffusivity D in the swollen layer can be assumed constant or an increasing exponential function of the point-wise solvent volume fraction: φ [22,35,36] where D G and D eq = D G exp γ are the minimum and maximum solvent diffusivities at the minimum φ G and at the maximum φ eq solvent volume fractions in the gel layer. The amount of absorbed solvent M(t) (per unit area of the film), at each time instant during the course of the swelling process can be evaluated as: where ρ s is the solvent mass density. When equilibrium (asymptotic) conditions are reached, the glassy phase completely disappears-i.e., GGI ∞ = 0-and the half-thickness of the fully swollen film L ∞ = GSI ∞ , in the absence of dissolution, reaches its asymptotic value.
Correspondingly, the asymptotic amount of absorbed solvent M ∞ (per unit area of the film) attains the form:

Poisson-Kac Transport Model for Case II Diffusion
The main goal is to model the Case II diffusion process by including the effects of polymer relaxation time. This can be addressed through the following steps: 1.
identify the solvent volume fraction φ(x, t) with the overall probability density function p(x, t) of the Poisson-Kac stochastic process-i.e., p( include in the dynamics of the partial waves p + (x, t) and p − (x, t) a further convective contribution accounting for the swelling velocity v sw (x, t) due to solvent penetration; 3.
describe the movement of the Glass-Gel interface and Gel-Solvent interface in terms of partial waves p ± (x, t); 4.
account for the effect of polymer relaxation time on the boundary condition at the Gel-Solvent interface.
At first, we consider the simplest case of a constant relaxation time t r and a constant diffusivity Di.e., a constant transition rate λ and a constant characteristic velocity b 0 , Equation (13). By considering that (i) the point-wise swelling velocity v sw (x, t) is assumed to be equal and opposite in sign to the diffusive flux, and that (ii) the "diffusive" flux in the Poisson-Kac formulation is given by Equation (10), it naturally follows that: The evolution equations for the partial densities p + (x, t) and p − (x, t) in the gel layer GGI(t) < x < GSI(t) attain the form: By introducing the dimensionless time τ = t/t r = 2λt and the dimensionless spatial coordinate z = x/L 0 , the transport equation for the partial waves p ± (z, τ) in the gel layer GGI(τ) < z < GSI(τ) reads as: where De is the Deborah number, representing the ratio between the relaxation time t r and the characteristic diffusion time t d = L 2 0 /D. Correspondingly, the Gel-Solvent interface GSI evolves according to the swelling velocity at GGI: The relaxation time t r also controls the relaxation of the solvent volume fraction φ at GSI(t) towards the asymptotic equilibrium value φ eq [22]. This can be accounted for, in the partial wave formulation, as: t r dp eq(t) dt = φ eq − p eq (t) =⇒ dp eq (τ) dτ = φ eq − p eq (τ) This implies a time-dependent boundary condition for the partial wave p − (z, τ)-namely: For the partial wave p + (z, τ), the boundary condition at the Glass-Gel interface GGI(τ) switches from a Dirichlet boundary condition p(z, τ) = φ G to an impermeability condition when the Glassy phase disappears-i.e.,: Meanwhile, the Glass-Gel interface evolves according to the Stefan condition: As for the parabolic transport scheme, the initial conditions for p ± (z, τ), GGI(τ), and GSI(τ) are GGI(0) = 1 − , GSI(0) = 1, and p(z, 0) = φ G that implies for the partial waves p + (z, 0) = p − (z, 0) = φ G /2. This approach can be generalized to include the effect of both a relaxation time t r (p) and a solvent diffusivity D(p) that are exponential functions of the point-wise solvent volume fraction φ(x, t) = p(x, t), as established by simplified versions of the free volume theory [25][26][27]. By setting in Equation (20): the case of an exponentially decreasing relaxation time and an exponentially increasing diffusion coefficient in addressed, i.e.,: By setting γ d = γ r /2 the case of an exponentially decreasing relaxation time and a constant diffusion coefficient D = D 0 = b 2 0 2 λ 0 can be described as well. For both cases, the transport equations for the partial waves and moving boundaries become: where τ 0 = 2λ 0 t and De 0 = b 0 2λ 0 L 0 2 is the Deborah number evaluated at p = φ G , representing the maximum value of the Deborah number during the course of the swelling process.

Analysis of Sorption Curves
The numerical solution of the transport scheme Equations (33)-(41) furnishes the spatio-temporal evolution of the partial waves p + (z, τ) and p − (z, τ), as well as the tempo-ral evolution of the moving interfaces GGI(τ) and GSI(τ) for constant relaxation time t r and solvent diffusivity D.
The transport equations for the partial waves p ± (z, τ) have been numerically solved with a home-made Matlab code implementing a finite-difference approach with (i) N = 5 × 10 3 discretization points for each partial wave, (ii) a second-order upwind representation of convective terms, and (iii) a proper rescaling of the space variablez(τ) = ((z − GGI(τ))/(GSI(τ) − GGI(τ)) to transform the moving boundary problem GGI(τ) ≤ z ≤ GSI(τ) into a fixed boundary problem 0 ≤z(τ) ≤ 1, ∀τ > 0 [37,38]. The resulting set of 2N + 2 ordinary differential equation has been numerically solved with the multistep solver ode15s with specified relative tolerance 10 −5 and absolute tolerance 10 −8 .  Figure 2A-C show a sharp discontinuity of the solvent concentration at the Glass-Gel interface [17] that moves towards the symmetry axis z = 0. When GGI(θ) reaches the symmetry axis, the discontinuity, for larger values of De-i.e., for larger relaxation times (Figure 2A,B)-is reflected back to the external surface while, for smaller values of De ( Figure 2C ) is rapidly smoothed away like in the parabolic case, solution of the transport scheme Equations (22)- (25).
The comparison between the parabolic and the hyperbolic case is presented in Figure 3B, showing the temporal evolution of the amount of absorbed solvent M(θ), normalized with respect to its equilibrium value M ∞ , for De → 0 (parabolic case, dot-dashed black curves) and for increasing values of De. It can be observed that, the larger De-i.e., the larger the relaxation time is, the slower the sorption curve that follows an "anomalous" θ 3/2 scaling before collapsing onto the sorption curve of the parabolic scheme when the boundary value p eq (θ) reaches its asymptotic value φ eq , as shown in Figure 3A. The classical θ 1/2 scaling of the standard parabolic scheme represents the envelope of sorption curves only for smallintermediate values of De, namely De ≤ 10 −2 . For higher values of De, the sorption curves exhibit a smooth transition from the θ 3/2 scaling to the asymptotic saturation behaviour and no θ 1/2 scaling can be detected.
A very similar behaviour is observed for the temporal evolution of the Gel-Solvent interface GSI(θ) shown in Figure 3C. Additionally, for the rescaled Gel-Solvent interface, slower dynamics and a net θ 3/2 scaling are observed for larger values of De. The larger De, the slower the swelling dynamics, and this also influences the dynamics of the Glass-Gel interface GGI(θ), as shown in Figure 3D. Indeed, for this case the disappearance of the Glass phase requires a longer time.
A more complex behaviour of the swelling system can be observed when the case of a non-constant relaxation time is analyzed. The case of an exponentially decreasing relaxation time t r (p) and a constant effective solvent diffusivity D is here addressed, see Equation (43) with γ d = γ r /2. The results of the numerical integration of the transport scheme Equations (44)-(48) for γ r = log 10, log 100 are shown in Figure 4A-D together with the results for a constant relaxation time γ r = 0, data already shown in Figure 3A-D. Sorption curves for γ r > 0 and high values of De clearly show good conformity to case II kinetics, feature V in the classification of sorption kinetics by Sanopoulou and Petropoulos [4] .
It should be observed that, for γ r = log 10, the relaxation time decreases by one order of magnitude from t r (φ G ) to t r (φ eq ) while, for γ r = log 100, it decreases by two orders of magnitude. It is therefore natural to expect that, the larger γ r the faster the sorption dynamics as well as the evolution of the Gel-Solvent interface, as shown in Figure 4B,C. However, it can be readily observed that the sorption curves for γ r > 0 do not settle down onto the asymptotic behaviour of the parabolic scheme Equations (22)- (25), as in the case of a constant relaxation time γ r = 0. The same phenomenon is observed for the temporal evolution of the Gel-Solvent interface, as depicted in Figure 4C.
This feature of the sorption process is intrinsically due to the nonlinear convective term, introduced and discussed in Section 3, Equation (21), arising from a characteristic velocity b 0 and a transition rate λ that are continuous functions of the overall density p, as introduced in the Generalized Poisson-Kac process to model a non-constant relaxation time t r (p). Indeed, in the case under investigation that accounts for the swelling velocity, the Kac limit for the overall density p(z, θ) satisfies the following parabolic equation:  By comparing Equation (49) with the corresponding Equation (22) of the parabolic transport model, we can easily recognize the presence of a new convective term, highly nonlinear, characterized by a point-wise negative velocity, that facilitates solvent penetration and speeds up the outward movement of the Gel-Solvent interface. More specifically, in the case of a constant solvent diffusivity (2γ d = γ r ), such as that presented in Figure 4A-D, Equation (49) attains the form: This convective term (the second term in the right hand side of Equation (50)) can be significant and amplified in interval absorption experiments, where the concentration jump (φ G − φ eq ) can be small. Its presence in the Kac-limit transport equation for the overall density can describe the transition from the S-shaped to the two-stage regime experimentally observed in interval absorption experiments [4].

Drug Release
This section investigates the influence of the polymer relaxation time on the release kinetic of a drug initially loaded in the thin dry film. In the absence of drug-polymer interaction, drug release in a 1-d swelling system can be simply modeled by a 1-d advectiondiffusion equation for the drug concentration c d (x, t), describing drug transport in the gel layer, along the preferential swelling direction x where D d is the drug effective diffusivity in the swelling film. Drug transport is strongly influenced by swelling dynamics not only because of the presence of the swelling convective term but mainly because of the moving boundaries GGI(t) and GSI(t). Specifically, the Glass-Gel interface GGI(t) controls the drug release rate from the Glass to the Gel phase, while the Gel-Solvent interface GSI(t) controls the length of the diffusive path for the drug to be released in the external environment. Indeed, the drug transport Equation (51) must be solved simultaneously with the equations describing the swelling dynamics that furnish, at each time instant, the point-wise swelling velocity v sw (x, t), as well as the position of the GGI(t) and the GGI(t) interfaces.
The boundary condition adopted for drug concentration c d (x, t) at the GGI(t) interface is the Stefan condition: where c 0 d is the initial drug concentration, supposed uniform in the dry film. A perfect sink condition at the Gel-Solvent interface-i.e., c d | GSI(t) = 0-is assumed without loss of generality.
The total amount of drug (per unit surface area) M d (t) released up to time t can be evaluated as: Correspondingly, given the perfect sink condition adopted, the total amount of drug released at equilibrium is M d ∞ = 2 c 0 d L 0 , equal to the total amount of drug initially loaded in the dry film.
For the sake of simplicity, here we analyze the case of a constant drug diffusion coefficient D d . The swelling velocity v sw and the temporal evolution of the GGI and GSI interfaces can be obtained from the solution of the parabolic transport scheme Equations (22)-(25) corresponding to De → 0 or from the solution of the hyperbolic transport scheme Equations (33)-(41) for different values of the Deborah number, corresponding to different polymer relaxation times for a constant solvent diffusivity D. A new parameter α needs to be introduced, representing the drug to solvent diffusivity ratio α = D d /D. Figure 5A-C show drug release curves for α = 1, 0.1, 0.01. When drug and solvent diffusivity are of the same order, i.e., α = 1 ( Figure 5A), drug release curves for De > 0 approach monotonically the drug release curve for De = 0. The larger De, the slower the release and the release curve for De > 0 is always below (slower than) the release curve for This behaviour significantly changes when the drug diffusivity is one order (α = 0.1) or two orders of magnitude (α = 0.01) smaller than solvent diffusivity. Indeed, Figure 5B,C clearly show that, when α = 0.1, the drug release curves for De > 0, for intermediate time scales, become faster and approach "from above" the release curve for De = 0. This effect is amplified in intensity and time duration for the smaller value α = 0.01 analyzed. In order to explain this phenomenon, we need to consider that drug release is strongly influenced by the movement of the two interfaces GGI and GSI. In particular, the Gel-Solvent interface controls the length of the drug diffusive path. The larger De, the slower the GSI velocity at short-intermediate time scale, the shorter the diffusive path and consequently the faster the release kinetics. Moreover, the smaller the drug diffusivity, the more important the effect of a shorter diffusive path on the release.
A shortcut estimation of the drug flux at the outer boundary GSI can be obtained from the macroscopic concentration gradient: shown in Figure 6 for the three values of α analyzed. Dashed lines represent the temporal evolution of the dimensionless macroscopic concentration gradient ∆c d /c 0 d ∆L/L 0 for De = 0, while continuous lines represent the same quantity for De = 10 −2 . It can be readily observed that, for α = 0.01, the macroscopic concentration gradient for De = 10 −2 is significantly larger than that for De = 0 (red curves) exactly in the time interval 10 −4 ≤ θ ≤ 10 −2 during which there is an "inversion" of the release curves (see Figure 5C, azure curve). This effect is still present but less intense for α = 0.1 (blue curves) and actually disappears for α = 1 (orange curves), in agreement with the observed behaviour of the release kinetics shown in Figure 5A,B. As a confirmation of this explanation, Figure 7 shows the behaviour of the drug release kinetics for a polymer relaxation time that is an exponentially decreasing function of the solvent volume fraction. In this case, the faster movement of the Gel-Solvent interface for De > 0 and γ r > 0, already presented and discussed in Section 5, implies a faster increase in the drug diffusion path and a consequent slow down of the release kinetics. Indeed, Figure 7 shows that the larger γ r , the smaller the overshoot of the release curve with a consequent faster collapse onto the parabolic release curve (De = 0).

Conclusions
The article investigates the influence of polymer relaxation time on sorption curves and drug release kinetics. Case II diffusion processes are addresses and analyzed through the theory of Poisson-Kac stochastic processes possessing finite propagation velocity.
The article provides a first, physically significant application of Poisson-Kac processes to moving-boundary problems associated with the swelling dynamics in polymeric matrices. Moreover, it corresponds to a physically relevant example in which the parameters of the process depend in a nonlinear way on the concentration of the diffusing species. This case has been treated theoretically in [13] by considering simple models of exclusion processes. Here, it finds a practically relevant application related to transport in polymeric systems. Moreover, the way of handling boundary conditions represents a generalization of the analysis developed in [39].
A preliminary analysis of sorption curves for a constant polymer relaxation time clearly shows that the larger the Deborah number De, the slower the sorption curve characterized by an anomalous θ 3/2 short-intermediate time-scale behaviour. The classical θ 1/2 scaling of the standard parabolic scheme represents the envelope of sorption curves only for small-intermediate values of De ≤ 10 −2 while, for higher values of De, the sorption curve exhibits a smooth transition from the θ 3/2 scaling to the asymptotic saturation behaviour.
A more complex and unexpected behaviour of the sorption curves has been observed when the case of an exponentially decreasing relaxation time is analyzed. Indeed, sorption curves for De > 0 are, as expected, slower than the corresponding sorption curve for De = 0 at short time-scales. However, in longer time-scales, there is an inversion and the sorption process becomes faster for De > 0 than for De = 0. This phenomenon is intrinsically due a nonlinear convective term, arising from the introduction of a characteristic velocity and a transition rate that are continuous functions of the overall density. This convective term, quantified in Equation (50), facilitates solvent penetration and speeds up the outward movement of the Gel-Solvent interface.
An equally unexpected behaviour has been observed for the drug release curves, when the influence of polymer relaxation time on release kinetics is investigated. Indeed, when drug diffusivity is one order or two orders of magnitude smaller than the solvent diffusivity, the role of polymer relaxation time (De > 0) is to speed up the drug release kinetics at a small-intermediate time scale, with the process being controlled by a smaller length of the diffusive path for the drug to be released in the external environment.
All these observations are extremely important in the analysis of experimental data of sorption and drug release curves and can reasonably explain many of the "anomalous" behaviours that cannot be described by employing the more classical parabolic transport schemes. The model does not introduce new fitting parameters with respect to those introduced in other hyperbolic models and can be applied in a wide range of Deborah numbers De ∈ [10 −4 − 10 3 ] without any stability issues for the numerical solution of the partial waves transport equations. Obviously, being a mechanistic model described by partial differential equations, its application to the analysis of experimental transient sorption data requires the numerical integration of the transport equations and the use of basic optimization tools for the best fit of transport parameters. A similar procedure can be found in [40], where the rate-type viscoelastic model of Camera-Roda and Sarti [21] is applied to reproduce the anomalous sorption of fluoropolymer-solvent systems.
What is relevant in tracing the connection between hyperbolic transport theory and stochastic processes is that the model equations used in this article are suitable for a direct stochastic description, meaning that all the results obtained in this article can in principle be obtained from the stochastic Lagrangian kinematics of solute particles. Within the scientific community interested in transport properties in polymeric systems, the stochastic interpretation of the constitutive equations with memory associated with the glassy-rubbery transition and with the influence of viscoelasticity to mass transport is missing. This article provides a rational way to fill this gap.
The PK approach can be naturally extended to two/three dimensional problems as well as to include more than one distinct structural relaxation process with different relaxation times.