A Multiscale Model for Solute Diffusion in Hydrogels

The number of biomedical applications of hydrogels is increasing rapidly on account of their unique physical, structural, and mechanical properties. The utility of hydrogels as drug delivery systems or tissue engineering scaffolds critically depends on the control of diffusion of solutes through the hydrogel matrix. Predicting or even modeling this diffusion is challenging due to the complex structure of hydrogels. Currently, the diffusivity of solutes in hydrogels is typically modeled by one of three main theories proceeding from distinct diffusion mechanisms: (i) hydrodynamic, (ii) free volume, and (iii) obstruction theory. Yet, a comprehensive predictive model is lacking. Thus, time and capital-intensive trial-and-error procedures are used to test the viability of hydrogel applications. In this work, we have developed a model for the diffusivity of solutes in hydrogels combining the three main theoretical frameworks, which we call the multiscale diffusion model (MSDM). We verified the MSDM by analyzing the diffusivity of dextran of different sizes in a series of poly(ethylene glycol) (PEG) hydrogels with distinct mesh sizes. We measured the subnanoscopic free volume by positron annihilation lifetime spectroscopy (PALS) to characterize the physical hierarchy of these materials. In addition, we performed a meta-analysis of literature data from previous studies on the diffusion of solutes in hydrogels. The model presented outperforms traditional models in predicting solute diffusivity in hydrogels and provides a practical approach to predicting the transport properties of solutes such as drugs through hydrogels used in many biomedical applications.


■ INTRODUCTION
Hydrogels are three-dimensional cross-linked polymer networks mainly composed of water (typically 70%−99%). 1 The polymeric network within the liquid provides hydrogels unique properties that make them particularly attractive for biomedical engineering applications. 2−5 Clinical use of hydrogels in regenerative medicine and drug delivery critically relies on the diffusion of solutes across the hydrogel matrix, where entrapped cells in 3D tissue engineering scaffolds must be supplied with oxygen and nutrients, 6,7 and drug delivery requires controlled release mechanisms. 1,8,9 While the translation of these three-dimensional networks to biological systems depends on the diffusion of solutes, quantitatively understanding this dynamic movement of these solutes within a hydrogel network has been challenging. 10,11,20−24,12−19 The diffusion of solutes in hydrogels is commonly modeled by one of the following three theoretical frameworks: (i) hydrodynamic theory, 25 which considers friction between the solute and the surrounding hydrogel matrix; (ii) f ree volume theory, 26 which assumes that the solute is transported via dynamic empty spaces between molecules; and (iii) obstruction theory, 27 which models the polymer net as a barrier for the diffusion of the solute with the liquid. While each of these models is successful at capturing features of experimental diffusion data under some circumstances, they are often not predictive because they are applied outside of experimental regimes for which they are applicable and/or contain unspecified empirical parameters. 28,29 Despite the substantial progress made during the past 40 years, a standard model that accurately predicts mass transport in hydrogels is still lacking. As a consequence, time and capital-intensive trial-and-error procedures are often employed to assess the viability tests of specific hydrogel applications. In this paper, a comprehensive model combining three primary theories for the diffusion of solutes in hydrogels at multiple length scales is presented. We call this framework the multiscale diffusion model (MSDM).
To test this framework experimentally, we employed a series of hydrogels made of poly(ethylene glycol) (PEG) as a model system on account of their broad use in both tissue engineering 30 and drug delivery applications. 31,32 To validate our model, we needed to quantify the size of the (sub)nanoscopic free volume holes present within these hydrogels. For that purpose, we used positron annihilation lifetime spectroscopy (PALS), a unique technique capable of measuring these molecular pores in biomaterials under wet conditions. To further test our model, which describes the solute as a hard sphere similar to most traditional models, we perform a meta-analysis of 66 distinct solute/hydrogel combinations reported previously by us and other research groups (Table S1). The model proposed here suggests that diffusion occurs through several complementary mechanisms, described independently in traditional models, and that all mechanisms must be considered together for accurate prediction of solute diffusion.
■ RESULTS Model Formulation. The diffusivity (D 0 ) of a solute with hydrodynamic radius (r s ) in a pure liquid is given by the hydrodynamic theory, as defined by the Stokes−Einstein equation: 33 where k b is the Boltzmann constant, T is the absolute temperature, and η is the dynamic viscosity of the liquid.
In a hydrogel, solute diffusivity is altered by the presence of the polymer chains, which form a network with open spaces between the chains. The size of these spaces, typically termed the characteristic "mesh size" of the network, can be nano-or microscopic, and the space between polymer chains is filled by aqueous solution. Here, ξ is defined as the correlation length, also known as the average mesh size of the network, 34 which can be experimentally measured (e.g., neutron scattering). Another property of relevance is the free volume, formed by dynamic empty voids between all molecules forming the hydrogel structure (e.g., between the water molecules, between the polymer molecules, and at the water/polymer interface). The average size of these pockets is quantified by the free volume void radius (r FV ), is atomic in scale, and hence is typically several orders of magnitude smaller than the mesh size ( Figure 1).
In this paper, it is postulated that the solute diffusivity normalized by the diffusivity in pure solvent, D/D 0 , will be dictated by the probability of diffusing (i) via free volume (FV) voids and/or (ii) alongside aqueous solution through a mesh of size ξ. These events are essentially mutually exclusive as they occur at distinct length scales and thus the probability that a solute will diffuse by both mechanisms simultaneously is essentially zero. Thus Based on free volume theory, solute diffusivity in a hydrogel with a given polymer volume fraction, φ p , normalized by the solute diffusivity in liquid, was developed by Lustig and Peppas 35 as follows: Here C is a sieving factor, and Y is the ratio between the critical volume required for the solute diffusion (approximated as the volume of the solute, V S ) and the available free volume per molecule in the aqueous solution inside the hydrogel (approximated as the average volume of the free volume voids in the water, V FVW ). Based on obstruction theory, which is derived from Fick's law, the probability of finding an open space of size ξ between polymer chains of radius (r f ) forming the hydrogel network was described by Amsden et al. 36 as follows: where the polymer mesh within the hydrogel is modeled as a motionless physical obstacle for diffusion of the solute. Bringing these probabilities together, the following equation is obtained: In this equation, a weighting factor A is used to combine solute diffusivities via free volume holes and/or through the mesh such that the total probability for the normalized diffusion to occur by any mechanism is quantified as a number between 0 and 1. Here, the prefactors account for the assumption that the diffusion via free volume dominates when the average radius of the free volume voids is comparable to the hydrodynamic radius of the solute (r FV ∼ r s ), whereas the mesh size becomes the limiting factor when the size of the solute is much larger than the free volume voids (r FV ≪ r s ). In addition, intermolecular forces are considered negligible in these systems, which is a reasonable assumption based on previous treatments of numerous hydrogel systems and distinct solutes. 35,36 Figure 1. Scale effects in solute diffusion in hydrogels. The diffusion of a solute within a hydrogel occurs via aqueous solution and through liquid-filled, nano-to microscopic open spaces between the polymer fibers or free volume (dynamic, subnanoscopic, empty voids between the molecules). Which of these mechanisms dominates diffusion depends on the ratio between the hydrodynamic radius of the solute and the radius of the free volume voids in the hydrogel.

Article
In this work we made use of the Gaussian error function, which is classically encountered in many problems in diffusion, as the prefactor. It is well-known that the radius of the free volume holes in many hydrogels and diverse biological tissues follow a Gaussian distribution with an average r FV and a distribution breadth denoted by σ VF . 37−41 For a particle with radius r s and free volume holes of r FV that are normally distributed, the corresponding error function is = ( ) This function describes the probability that a diffusing particle of r s will find a single intermolecular space of r FV , which falls in the range of 0−1 since r r FV s is always a positive value.
Consequently, the probability of the particle diffusing through other mechanisms (i.e., not via free volume holes) is captured by the complementary error function − = ( ) The weighting factor we made use of here is therefore not based on empirical observations and is instead based on a physical theory, encountered in integrating the normal distribution of the radii of the free volume voids normalized by the solute radius. With the inclusion of these functions, the previous equation can be rewritten as follows: In this equation, the solute is modeled as a hard sphere and the free volume voids in the hydrogel as empty spheres with radius r FV (r FVW in the case of water). While many techniques have been developed to characterize the mesh size and radius of the solute, allowing these parameters to be easily obtained, the free volume of a system is more difficult to characterize and therefore frequently overlooked. In a previous study, the mesh size in our hydrogels was estimated by spherical indentation. 42 Though, to validate our MSDM model, shown in eq 6, it was necessary to measure r FV within the PEG hydrogels. For this end, we used positron annihilation lifetime spectroscopy. Positron Annihilation Lifetime Spectroscopy (PALS). PALS is a nondestructive technique capable of probing atomicsized open spaces. The lifetime of positrons implanted in materials is measured, which yields information about the electron density. Particularly, orthopositronium (o-Ps) is a system with a radius of 1.59 Å, consisting of an implanted positron and an electron from the material with parallel spins, formed in the free volume voids (Figure 2A). The o-Ps lifetime is correlated to the size of the free volume voids inside a material, including hydrogels. This difference in size is reflected in the positron lifetime spectra obtained (see representative spectra for PEG5 and PEG25 in Figure 2B). By deconvoluting these raw spectra, we obtain the characteristic average o-Ps lifetime and distribution ( Figure 2C). Finally, the o-Ps lifetime is correlated to the average free volume radius (r FV ) and distribution (σ VF ) by a well-established model (Table 1). 40,43,44 The average radius of the free volume voids, r FV , was similar for the three most dilute PEG hydrogels tested and was found to be very similar to the free volume in water r FVW (2.69 Å). 42,45 In contrast, the r FV determined for the hydrogel system with the highest polymer volume fraction, PEG25, differed by more than 20% from the average radius of the free volume in water. Presumably, the higher polymer content in the latter system introduces significantly more free volume as it contains a greater amount of surface area both between two polymer chains and between water and the polymer chains.
Validation of the Model. To validate our approach, we used the MSDM model to predict the diffusivity (D) of solutes of various sizes in the series of PEG-based hydrogels with distinct mesh sizes described in the Materials Section ( Figure  3). The values for the cross-sectional radius of the hydrated PEG fiber (r f ) and the average free volume void size in water (r FVW ) were taken from the literature as 5.1 and 2.69 Å, respectively. 45,46 It is clear that the predicted diffusivity values decrease monotonically with increasing solute radius ( Figure  3), in accordance with physical expectations. We compared these predictions to experimental diffusivity values previously obtained 42 by fluorescence recovery after photobleaching (FRAP) for dextran solutes of different sizes (r s = 1.9, 3.5, and 19.4 nm, corresponding to molecular weights of 4, 20, and 2000 kDa, respectively) in these hydrogels ( Figure 3).
Additionally, we compared the MSDM model to (i) free volume theory as defined by Lustig and Peppas 35  ). 36 Figure 4 shows the predictions of the three different models compared to experimental data discussed previously, normalized by their diffusivity in pure

Macromolecules
Article liquid (D/D 0 ). 47 Again, the values for the cross-sectional radius of the hydrated PEG fiber (r f ) and the average free volume void size in water (r FVW ) were taken from the literature as 5.1 and 2.69 Å, respectively. 45 To further validate the MSDM model and quantify the differences in the accuracy between models, we performed an analysis of data obtained by other research groups in previous studies. 36,46,48,49 In total, 66 distinct solute/hydrogel combinations using two different hydrogel systems (PEG and alginate) were analyzed (Table S1). Figure 5 shows experimentally obtained solute diffusion coefficients in the hydrogels, plotted versus the prediction from eqs 3, 4, and 6. For this analysis, as free volume hole size data was not available for many of the data sets, two approximations were made: (i) Y = 1 for the free volume theory predictions as previously suggested by the Lustig and Peppas 35 and (ii) r FV ≈ r FVW for the MSDM model. When access to PALS equipment is not available, the approximation that r FV ≈ r FVW is appropriate because the free volume holes in hydrogels deviate from those in water only at high polymer concentrations (Table 1). With this, the MSDM model takes into account four independent parameters, while obstruction theory has three and free volume theory has two independent parameters. The MSDM model (R 2 = 0.885, calculated with respect to the perfect prediction line) outperforms obstruction theory (R 2 = 0.743) and free volume theory (R 2 = 0.787) in predicting the experimentally obtained diffusivity more accurately. Moreover, a study of the residuals ( Figure 3B) shows that the MSDM mean square error value (MSE) is ∼4 times smaller (185.7) than that for the obstruction theory (710.0) and 3 times smaller than that for free volume theory (564.8).
Typical drugs range in diameter from 1 nm (small molecules) to 5 nm (antibodies). For potential drug delivery applications of our model, we analyzed the residuals for the Values of the intensity of the orthopositronium signal (I o-Ps ), orthopositronium Lifetime (τ o-Ps ) and distribution (σ o-Ps ), average radius of the free volume voids (r FV ), average volume of the free volume voids (V F ) and distribution (σ VF ), and for the different PEG hydrogels tested. χ 2 are the χ 2 test values.

Macromolecules
Article predicted diffusivity values with different models versus solutes with radii within this typical range (r s from 0.5 to 5 nm; Figure  S1A). The outcome of this analysis highlights the accuracy of the MSDM model for solutes with sizes that lie within the typical range of size range for drug molecules: 1−5 nm.

■ DISCUSSION
Here, a predictive model for the diffusion of solutes in hydrogels was formulated combining three traditional theoretical frameworks: hydrodynamic theory, free volume theory, and obstruction theory. It was demonstrated that the MSDM model successfully predicts the diffusion of dextran solutes of multiple sizes in chemically cross-linked PEG and physically cross-linked alginate hydrogels with multiple mesh sizes (66 distinct combinations in total), with an increased qualitative and quantitative accuracy compared to traditional models. Previous studies 1,50 have employed hydrodynamic theory in the form of the Stokes−Einstein equation to predict the diffusivity in hydrogels with a mesh size larger than the solute size. However, this approach appears to be inaccurate, as it overpredicts the experimental diffusivity in hydrogels even if the mesh size exceeds the solute radius by a factor of 18 (in Figure 4, this case would be represented by a horizontal line of D/D 0 = 1). On the other hand, free volume theory tends to underestimate the diffusion coefficient unless the solute size is large. Moreover, when using the approximation Y = 1 proposed by Lustig and Peppas, 35 the model predicts negative (and, as a consequence, unphysical) values for the diffusion coefficient of solutes with radius 19.4 nm (Table S2). In addition, the MSDM model is notably more accurate (R 2 = 0.885, calculated with respect to the perfect prediction line) than the equation developed by Amsden et al. 36 based on the obstruction theory (R 2 = 0.743).
Most drug molecules have a hydrodynamic radius that falls between 0.5 and 5 nm. In this solute size regime of paramount importance for drug delivery applications, the diffusion behavior is fundamentally not captured by traditional models ( Figure S1). By combining these traditional models to create a comprehensive description of solute diffusion, the MSDM model is remarkably more accurate than the any of these models individually. The mean square error is 3−4-fold smaller than traditional theories.
As in traditional frameworks, the diffusivity (D) predicted in the MSDM model decreases monotonically with increasing solute radius (Figure 3). The MSDM theory does not contain any term nor fitting parameter that is empirical. Our approach is entirely probabilistic in nature and begins from the conservation of probability (eq 2), which captures two independent modes of diffusion: (i) via free volume voids and/or (ii) alongside water through the polymer mesh. Both of these diffusion modes contribute to the total probability. The model presented here describes the competition between different mechanisms of diffusion, and three regimes emerge. First, when r s ≈ r FV (and therefore r s ≪ ξ), the MSDM model reduces to the free volume theory: In this regime, the solute diffuses via dynamic free volume voids formed in the hydrogel. Here, the probability of finding a free volume void decreases with an increase in the solute size. Second, for larger solutes, the diffusion via this mechanism becomes increasingly unlikely, and a local minimum in the D/ D 0 prediction is reached (r s ∼ 1 nm for the hydrogels tested; see Figure 4). In this intermediate zone, the probability of the solute to diffuse within the hydrogels via the liquid (instead by free volume voids) starts to gain importance. As the mesh size is still much larger than the solute, it allows its passage. As the solute size increases, both the values for D and D 0 decrease, but the D/D 0 ratio increases, as there exists a superposition of two diffusion mechanisms in this regime. Finally, a local maximum is reached when the solute size starts to be comparable to the mesh size (r s ≈ ξ and therefore r s ≫ r FV ), and the equation presented in this study reduces to the obstruction theory: Here, the probability of finding an aperture in the mesh decreases with a further increase in solute size. In other words, this local minimum and maximum indicate a transition in the hinder mechanism of the diffusing particle. Importantly, this local minimum and maximum has been previously shown to be a physical phenomenon, as observed by Hoh and Zia. 51,52 Thus, depending on the size scale of solute, the MSDM accounts for different mechanisms of diffusion, as summarized in Figure 6.
Overall, the MSDM model describes different dominant types of molecular transport depending on the solute size. This model provides valuable information not only for a theoretical understanding of the diffusion in hydrogels but also for practical applications in tissue engineering or drug delivery applications, by potentially predicting more accurately the diffusion of solutes such as growth factors, nanoparticles, or drugs in hydrogels. It is also important to note that solute−gel  42,46,48 and alginatebased hydrogels 36,49 more accurately than eqs 3 and 4free volume theory (triangles) and obstruction theory (circles). Each color of data points represents a different study of the meta-analysis. In the residuals versus experimental diffusivity plot, the perfect prediction is shown as a black dotted line.

Macromolecules
Article interactions, including hydrophobic effects, H-bonding, electrostatic, and van der Waals interactions, are present and of considerable importance in several systems. While few studies have explored drug−hydrogel interactions, this could include up to a 25% decrease in diffusivity. 53 However, inclusion of such interactions would be incredibly specific to each material, which is not the intention of our model that we hope can be broadly applicable to many systems. 46,54 We also want to highlight that for drug delivery systems the process of drug release is dominated by diffusion when the mesh is larger than the drug due to domination of steric interactions; in these cases, drug−gel interactions are not even necessary. Thus, many principles introduced here could be applied to drug delivery. 1 In addition, the MSDM model could be used in future studies as a baseline to quantify gel−solute interactions.
The favorable experimental validation of the model presented here is limited to dextran solutes−PEG hydrogels (considered as inert). To fully understand the strengths and weaknesses of the MSDM model, future work across a wide range of hydrogels and solutes is necessary. To apply this model when PALS data are not available, it could be assumed that the available free volume voids in the aqueous solution are much more numerous than those surrounding the polymer network in hydrogels. In other words, the average free volume void size in the hydrogels could be approximated as the one in water. The MSDM model was tested including this approximation (see Figure S2) and obtained good results (R 2 = 0.93, calculated with respect to the perfect prediction line for the 12 experimental points, whereas without the approximation the obtained value was R 2 = 0.94). The main limitations of the MSDM model, similar to each of the traditional models described above, are that the solute is described as a hard sphere (thus precluding possible shape effects) and that it does not consider possible intermolecular forces between the solute and the polymer. Future studies will investigate the utility of the MSDM model in more complex aqueous systems such as living cells as well as in organogels.

■ CONCLUSION
The MSDM model presented here is capable of integrating three different mechanisms for the solute diffusion in hydrogels occurring on various size scales. Compared to previous frameworks, this combined model leads to more accurate predictions of the experimentally measured diffusivity, determined by dextran solutes of different sizes in PEG hydrogels and alginate-based hydrogels with different mesh sizes and free volume properties. Development of tunable and well-defined materials in clinical use, regenerative medicine, and drug delivery relies critically on the ability of solutes (e.g., nutrients or drugs) to move through the material scaffold. The model presented could offer a way to design hydrogels for tissue engineering or drug delivery applications with tailored solute transport properties. The dramatic reduction in costs of fabrication and testing may pave the way to a larger use of hydrogels in biomedical commercial applications. This study can also be a fundamental advance in understanding the physics behind diffusion of solutes in hydrogels.
■ METHODS Hydrogel Preparation. Poly(ethylene glycol) dimethacrylate (PEGDMA) monomers with a molecular mass of M n = 1000 g mol −1 as specified by the manufacturer (Polysciences, Inc., USA), were used to fabricate PEG hydrogels with various mesh and free volume void sizes. PEGDMA was dissolved in phosphate-buffered saline (PBS) in weight concentrations of 5, 7, 10, and 25 wt % called PEG5, PEG7, PEG10, and PEG25, respectively. Ammonium persulfate (APS) was used as a free-radical initiator, and N,N,N′,N′tetramethylethylenediamine (TEMED) was used as an accelerator. APS volumes of 150 μL (10 wt % in H 2 O) and 75 μL of TEMED were added simultaneously to 10 mL of the precursor solution, vortexed for 10 s to ensure thorough mixing, poured into cylindrical molds 35 mm in diameter, and left overnight to cross-link. Prior to testing, the gels were placed in PBS for 24 h at room temperature to reach swelling equilibrium and to allow the removal of nonreacted monomers.
Positron Annihilation Lifetime Spectroscopy. Two samples from each hydrogel sandwiched the 22 NaCl positron source of about 15 μCi from PerkinElmer (USA), protected by two 7.5 μm Kapton foils and sealed by a CAPLINQ (Canada) double-sided sticky Kapton tape. The PAL spectrometer is composed by two H1949-50 Hamamatsu (Japan) photomultipliers added to two BC-422 plastic scintillators from Saint Gobain (USA) performed vertically inside a Radiber S.A. (Spain) FFD-1402 refrigerator. All the electronic modules were purchased from ORTEC (USA). The temperature of When the solute size is comparable to the free volume pockets, the solute will diffuse though free volume. (B) When the solute is substantially larger than the free volume, it will diffuse with the liquid within the hydrogel and cross the mesh. (C) When the solute is larger than the mesh size, the diffusion will be hindered by the polymer fibers of the mesh.

Macromolecules
Article the samples was kept constant at 25°C by a SALICRU (Spain) variable power supply, a 3508 temperature controller from Eurotherm (United Kingdom), a Watlow Europe (Germany) 100 W FIREROD cartridge heater, and a TC S.A. (Spain) PT-100 CS5 (1|5) temperature sensor introduced in an aluminum sample holder. The resolution function was 258 ps, and the detection rate was ∼70 counts s −1 . Each positron lifetime spectrum was obtained by ∼3 million counts and analyzed with the LT_polymers software. The Kapton contribution was extracted (19.8%, 0.382 ns) prior to the decomposition of each spectrum into three components.
where r 0 = r VF + Δr VF and Δr VF is an empirical parameter fitted as