Predicting Drug Release From Degradable Hydrogels Using Fluorescence Correlation Spectroscopy and Mathematical Modeling

Introduction

Predicting release from degradable hydrogel networks is highly significant for a multitude of applications: designing multilayered or multicomponent drug delivery devices; controlling drug or other molecule delivery by pore size, pore size distribution, affinity, or other interactions (Lin and Metters, 2006), facilitating nutrient and gas exchange in three-dimensional cell scaffolds for tissue engineering (Leddy et al., 2004), controlling solute diffusivity via crowding and confinement to enhance extracellular matrix production by cells to create tissues ex vivo (Chen et al., 2011), and many others.

For example, hydrogels are at the forefront of an ever-growing drug delivery industry, estimated to become a $90 billion market by 2025 (Finn, 2016). Drug delivery methods range from oral, topical, intravenous, to—more recently, localized and targeted delivery. The increasing number and diversity of medications, such as small molecules, biologics, or cells, require more efficient and customizable delivery. Localized and controlled delivery approaches involving degradable hydrogels have proven particularly beneficial (Kurisawa et al., 2005; Bhattarai et al., 2010; Jain et al., 2017b).

Hydrogels are excellent candidates for drug delivery because of their biocompatibility, bioinertness, and ability to preserve the activity of biomacromolecules (Peppas et al., 1999; Lee and Mooney, 2001; Lutolf and Hubbell, 2003; Zustiak and Leach, 2011; Zustiak et al., 2013; Jain et al., 2017c). Degradable hydrogels are particularly advantageous because drug release can be controlled via degradation and there is no need to remove a device once the payload is depleted (Zhao and Milton Harris, 1998; Kim and Park, 2002). A common hydrogel used in drug delivery is polyethylene glycol (PEG), which is inherently non-degradable. However, various degradable moieties in the form of reactive end groups may be added to PEG to make it degradable (Lutolf and Hubbell, 2003; Zustiak et al., 2013; Jain et al., 2017a). Additionally, PEG can be synthesized with various chain lengths and multifunctionalities to further control degradation (Jain et al., 2017a).

While degradable hydrogels are increasingly sought, most studies developed to predict solute release from hydrogels focus on non-degradable polymeric matrices, or for simplicity ignore degradation and the associated changes in release. For example, some studies have focused on a “moving boundary problem,” where hydrogels are swellable (Brazel and Peppas, 2000; Fu and Soboyejo, 2011). However, these hydrogels, although swellable, do not degrade and thus the model does not fully predict release. There are also studies on predictive release from oral delivery systems such as hydroxypropyl methylcellulose compacted to form a tablet (Siepmann et al., 1999; Siepmann and Peppas, 2012). However, these systems are based on surface dissolution rather than bulk degradation, as would be the case of a hydrogel matrix. In another example, Mason et al. (2001) developed an exponential analytical model for release of model protein bovine serum albumin (BSA) from a degradable PEG-polylactic acid co-polymer hydrogel. However, this model is limited to the hydrogel described in the study, as one of the input model parameters depended on experimental measurements of hydrogel degradation products.

In this study, we developed a mathematical model to predict solute release profiles from degradable hydrogels using a simple, analytical approach based on a Fickian diffusion model. We experimentally measured diffusivity of two model proteins, namely BSA and immunoglobulin G (IgG) using fluorescence correlation spectroscopy (FCS) and the thickness of a PEG hydrogel over 5 h. These components were then input to a MATLAB algorithm to extract a predicted release profile. The predicted profile was compared to a release profile obtained using traditional bulk release methods and showed good agreement. We further discuss simplifications to the model, which could be implemented without sacrificing prediction accuracy.

Theoretical and Computational Model

We assumed that the mechanism by which a protein solute is released from the hydrogel is Fickian diffusion. Due to the flat geometry of the hydrogel slab, a 1D partial differential equation (PDE) model for the concentration c of the protein was deemed adequate. This led to the PDE:

where D(t) is the time-varying diffusivity of the protein solute. The diffusivity is time-varying since, as the hydrogel degrades with time (see Figure 1 for a schematic of bulk degradation), the mesh size increases, likewise increasing diffusivity.

Figure 1. (A) Chemical structures of 4-arm PEG-acrylate macromer and the PEG-dithiolglycolate crosslinker. Schematic depicting (B) the formation of a PEG hydrogel with an embedded protein solute and (C) its degradation over time. The embedded protein is entrapped within the mesh network of the gel and as bulk degradation occurs, the mesh size of the hydrogel increases, allowing the proteins to diffuse more freely.

In addition to diffusivity, the geometry of the gel slab is also time-varying as the gel degrades. This is reflected in the time-varying spatial domain of the PDE, which is −h(t) ≤ x ≤ h(t), where h(t) is the time-dependent half-thickness of the gel.

The PDE was subject to the boundary conditions, where solute concentration was assumed to be zero at gel boundaries, c(−h(t), t) = 0 = c(h(t), t). This assumption was based on the fact that the releasate solution surrounding the hydrogel was stirred during bulk release experiments (preventing formation of concentration gradients), and gel volume was ~100-times smaller than the volume of the releasate solution (creating an infinite sink environment). In addition, the PDE was subject to an initial condition, where the protein solute was uniformly distributed throughout the gel at t = 0, c(x, 0) = c₀, −h(0) ≤ x ≤ h(0). To simplify the time-varying geometry of the gel, a change of variables of the form $y=\frac{x}{h(t)}$ was used to transform the domain from −h(t) ≤ x ≤ h(t) to −1 ≤ y ≤ 1, changing the PDE in Equation (1) to:

Subsequently, release rate $\dot{Q}(t)$ was formulated in terms of the concentration of the protein solute as:

The solution of the PDE model requires knowing D(t) and h(t) as functions of time. This can be accomplished in two ways. One approach is to develop a theoretical model to predict D(t)and h(t). This involves modeling how the hydrogel swells and degrades in time and how this in turn affects evolution of mesh size as well as hydrogel geometry as a function of time. Knowledge of mesh size at a given time may then be used to predict diffusivity at that time. See (Donovan et al., 2016) for instance on how to relate mesh size to diffusivity, described by us previously. However, reliably modeling the gel degradation is challenging. An alternative approach is to measure D(t) and h(t) at certain snapshots in time and then use curve fitting to obtain these as continuous functions of time for input to the PDE model. In this work, we tested this second approach.

To obtain continuous functions D(t) and h(t) from time snapshots of the solute diffusivity (measured by FCS) and the hydrogel slab thickness, time snapshot data can be fit with splines or some other family of curves using a weighted least-squares method. We decided to try using both quadratic splines as well as the exponential family:

determined by three positive parameters, namely a, b, and c.

The numerical solution of the PDE in Equation (2) was accomplished as follows. First a spatial discretization using centered finite difference approximations was employed to obtain a system of ordinary differential equations (ODEs). Then this system was solved via the MATLAB ODE solver ode23s. The spatial integration in Equation (3) was approximated by Simpson's rule, and $\dot{Q}(t)$ was integrated in time with the MATLAB ODE solver ode23s along with the system of ODEs. Software was developed in MATLAB for the weighted least-squares fit with splines and the lsqnonlin program in MATLAB was used for the weighted least-squares fit with the exponential family.

To test whether the time-varying nature of diffusivity and hydrogel thickness affected release profiles, we also computed release profile predicted by the PDE model corresponding to constant D- and h-values, where the constant values were calculated as an average of all time points. Moreover, we decided to assess feasibility of predicting release profiles from measurements based on far fewer time snapshots. We only used the FCS measured diffusivities and hydrogel thickness measurements obtained at three carefully chosen time points, fitted the exponential family of curves for D(t) and h(t), and used these as inputs to the PDE model to compute predicted release.

Materials and Methods

Materials

4-arm PEG-acrylate (4-arm PEG-Ac; 10 kDa) was obtained from Laysan Bio Inc. (Arab, AL). The dithiol crosslinker PEG-dithiolglycolate (PEG-DD1; 3.4 kDa) was synthesized via Fischer esterification reaction as reported by us previously (Zustiak and Leach, 2010; Jain et al., 2017a). Fluorescent Dye Removal Columns were purchased from Thermo Scientific (Waltham, MA). Perfusion chambers were purchased from Sigma Aldrich (St. Louis, MO). Silicone isolator sheets (0.5 mm thick) from Grace Bio Labs (Bend, OR) were used as spacers. Bradford protein assay was obtained from Bio-Rad (Hercules, CA). BSA (66 kDa), Immunoglobulin G (IgG; 150 kDa), Atto 655-NHS ester, and all other reagents were procured from Millipore Sigma (St. Louis, MO) unless stated otherwise.

Hydrogel Preparation

Hydrolytically degradable hydrogels were prepared by combining 4-arm PEG-Ac and the PEG-DD1 crosslinker in an equimolar ratio of reactive groups Ac:SH. A hydrogel matrix was then formed via Michael-type addition—a highly specific and mild gelation chemistry—between the reactive Ac and SH groups. Briefly, 4-arm PEG-Ac and PEG-DD1 were dissolved in a 0.3 M, pH 7.4 triethanolamine (TEA) buffer to make a 10% w/v PEG hydrogel precursor solution (Figures 1A,B). The precursor solution was then pipetted onto a parafilm-covered glass slide in 50 μL droplets. One millimeter thick silicon spacers were placed at the ends of the glass slide and a second parafilm-covered glass slide was used to sandwich the droplets to form thin slabs. The slabs were incubated at room temperature for 30 min to ensure complete gelation. Note that gelation time for gels made with 4-arm PEG-Ac and PEG-DD1 crosslinker at pH 7.4 is ~2 min (Jain et al., 2017a).

Measurement of Hydrogel Swelling Ratio and Mesh Size

Hydrogels were characterized for swelling ratio (Q_M) and mesh size (ξ) as a function of time. Hydrogel mass was measured immediately after fabrication to obtain the relaxed mass (M_R). Afterwards, gels were stored in phosphate buffered saline (PBS) and re-weighed at 2, 4, and 6 h to obtain the swollen mass (M_S). Gels were then placed in an oven at 60°C overnight and re-weighed to obtain the dry mass (M_D). Q_M was calculated at the time points mentioned above as M_S/M_R. ξ was calculated according to the Flory-Rehner theory as follows (Canal and Peppas, 1989; Zustiak and Leach, 2011):

where υ_2,s is the polymer volume fraction in the swollen state, C_n is the characteristic ratio for PEG, M_C is the average molecular weight between two adjacent crosslinks, M_r is the molecular weight of a PEG repeat unit, and l is the average bond length between the C-C and C-O bonds in a PEG repeat unit.

Note that the Flory-Rehner theory was developed for equilibrium swollen gels. Here, equilibrium swelling could not be established as the gels swelled continuously with degradation. Nevertheless, we expect the theory gives good indication of how ξ changed with time.

Measurement of Hydrogel Geometry

Hydrogels were placed in microcentrifuge tubes with 1x PBS for swelling. At specified time points gels were removed from PBS and the diameters were measured using calipers, while gel thickness was measured using a micrometer. Change in hydrogel geometry upon swelling is depicted in Figure 1C.

Measurement of Protein Bulk Release and Diffusivity

Hydrogels were prepared as mentioned previously with the following modification: 2% w/v of protein (BSA or IgG) was encapsulated by adding it directly to the hydrogel precursor solution prior to gelation (Figure 1B). For release experiments, the fabricated protein-loaded gels were transferred to a centrifuge tube with 5 mL of pre-warmed 1x PBS and placed on a shaker platform (Labquake Shaker, ThermoFisher, Waltham, MA) in an incubated environment at 37°C. A 1 mL releasate sample was taken at specific time points until degradation and replaced with 1 mL of fresh PBS to maintain a sink volume of 5 mL. The releasates were stored at −20°C until all time points were collected. Releasates' protein content was determined using Bradford protein assay following the manufacturer's protocol. An effective bulk diffusion coefficient (D_e) was calculated using the following relation for short release times (M_i−/M_∞ < 0.6) (Ritger and Peppas, 1987):

where M_i is protein content at each time point, M_∞ is protein content at degradation, t is time, and δ is the half-thickness of the gel slab.

Measurement of Protein Diffusivity Using Fluorescence Correlation Spectroscopy

For FCS measurements proteins were labeled with Atto 655-NHS ester fluorophore as per the manufacturer's protocol. Unbound fluorophores were removed using Fluorescent Dye Removal Columns as per the manufacturer's protocol with 95% efficiency. For FCS measurements, all gels were loaded with fluorescently labeled protein as described previously. Hydrogels were placed in a solution of Atto 655-labeled protein with the same concentration as in the gel (2% w/v) for the entirety of the experiment to avoid concentration gradients. At specific time points, hydrogels were removed from the soaking solution, gently blotted to remove excess solution, and placed on a #1 coverslip for FCS measurements. To avoid evaporation during measurement, a few drops of the soaking solution were pipetted on top of the gel, which was then covered with a custom lid.

The FCS instrument was calibrated using 0.2 nM Atto 655 dye in PBS. A 648 nm ps pulsed laser was used at an optical power of ~11.4 μW. A 640/LP filter was used along with a 655 dichroic mirror to obtain measurements for 120 s per sample. A minimum of 3 measurements at different locations were performed per hydrogel.

The obtained autocorrelation function G(τ) for each measurement was fitted by the following equation (Magde et al., 1974):

where N is the number of fluorescent particles, τ_D is the diffusion time, p = r_o/z_o is an instrumental constant, r_o is the radius of the focused laser beam spot, and z_o is the axial length of the focused laser beam spot.

For two non-interacting, diffusing solutes Equation (7) can be rewritten as (Michelman-Ribeiro et al., 2007; Zustiak et al., 2010):

where m₁ and m₂ are related to the quantum yield and average number of each diffusing species and τ₁ and τ₂ are their respective diffusion times.

Here, the two-component autocorrelation function provided a more accurate fit as the gel precursor solutions were prepared with protein solutes that were exogenously labeled. Hence, they contained some un-reacted fluorophores in addition to fluorophore-labeled proteins. Note that τ_D for the free fluorophore was measured separately (and determined from the single component fit in Equation 7) and used as a fitted parameter (τ₂) in Equation (8) to determine the diffusion time of the protein (τ₁). Additionally, the autocorrelation function was fit using a Triplet Extended (3D) model to account for the possible excitation of molecular triplet states at higher laser intensities. Lastly, the autocorrelation function was normalized as follows:

where G(τ_D) is the value of the Equation (8) at each time point and G(τ₀) is the value of Equation (8) at the initial time point. The effective tracer diffusion coefficient for each protein was calculated from τ_D as (Zustiak et al., 2010):

Statistical Analysis

Results are reported as mean averages with error bars of ± one standard deviation of triplicate samples from three independent experiments. A two-tailed Student's t-test was used to compare two samples. Differences between data sets were considered significant when p < 0.05.

Results

Hydrogel Swelling and Mesh Size as a Function of Time

As expected, the Q_M and mesh size of the hydrogel continuously increased with degradation (Table 1). Note that degradation was complete at ~8 h; however, at ~7 h the gels lose their shape and start fracturing. Hence, all measurements were conducted until 6 h, when the gels still had sufficient integrity to be handled. The Q_M increased by 7 times over the 6-h measurement. Note that for similar but non-degradable PEG gels, equilibrium swelling is achieved around 2 h (Zustiak et al., 2010). The lack of swelling saturation and the continuous increase in Q_M indicate continuous bulk degradation. Similarly, mesh size increased continuously by 1.3-times over the course of the measurement.

Table 1. Hydrogel characterization properties over the lifetime of the hydrogel. Errors represent one standard deviation.

Measurements of Solute Diffusivity Using FCS

FCS was used to enable real-time in situ measurements of solute diffusivity in the hydrogel as a function of time. Figure 2 shows representative autocorrelation curves for both labeled model protein solutes used in this study, namely BSA and IgG. As expected, for both proteins the autocorrelation curves at 5 h were shifted to the left compared to the curves for 0 h, indicating substantially smaller diffusion times at the later time points (Figures 2A,C). The residuals indicate an excellent agreement between the raw autocorrelation data and the fit (Equation 8). Overall, diffusion times were higher at all time points for the larger IgG protein, compared to the smaller BSA. For BSA, the τ_D-values were 3.3 ± 1.3 ms at 0 h and 1.4 ± 0.5 ms at 5 h. For IgG, the τ_D-values were 4.2 ± 0.7 ms at 0 h and 2.2 ± 0.5 ms at 5 h. Note that τ_D is inversely proportional to diffusivity.

Figure 2. Representative normalized autocorrelation curves 0 and 5 h for fluorescently labeled (A) BSA and (C) IgG. The normalized residual curves indicate agreement between the autocorrelation curve and the fit for (B) BSA and (D) IgG (n = 9).

Measurements of Bulk Release and Diffusivity

The fractional bulk release profiles for both proteins embedded in the hydrogel are shown in Figure 3A. BSA released faster from the gel than IgG, exemplified by the higher slope of the BSA fractional release profile. This was expected, because BSA is smaller (hydrodynamic radius = 3.5 nm) than IgG (hydrodynamic radius = 4.7 nm) (Zustiak et al., 2010). Both proteins' hydrodynamic radii are smaller than the mesh size of the hydrogel (13–17 nm; Table 1), facilitating release. Calculation of the bulk diffusion coefficients using Equation (6) showed the D_e for BSA to be 1.8-times higher than that for IgG (Figure 3B).

Figure 3. (A) Fractional release profiles of BSA and IgG from PEG hydrogels for small diffusion times (M_i/M_∞ < 0.6). The dotted lines represent the trendlines for the fractional release profiles of BSA (R² = 0.999) and IgG (R² = 0.983) (n = 6). (B) Calculated bulk effective diffusivities for BSA and IgG. A student's two-tailed t-test was performed. ^*Designates significant differences (p < 0.05, n = 6).

Validation of the Mathematical Model With Experimental Data

Tracer diffusivities (measured by FCS and computed from Equation 8) of labeled BSA and IgG continuously increased with gel degradation (Figures 4A,B). Diffusivity of labeled BSA exhibited a higher increase of 2.63-times during the 5 h of measurements, while the diffusivity of labeled IgG increased 1.93-times. Gel thickness increased 1.28 times (Figures 4C,D). The increase in thickness was abrupt in the first hour but mostly unchanged for the rest of the measurement time. Note that the gel also increased in diameter upon swelling. However, change in diameter was not used in the model, which was based on a 1D diffusion as explained earlier.

Figure 4. The spline fits calculated using MATLAB for the diffusion coefficient values obtained via FCS for fluorescently labeled (A) BSA and (B) IgG, and (C) the measured hydrogel thickness as a function of time (n = 9). (D) Images showing the increase in diameter and thickness of the hydrogel with time (blue coloring was used for visualization).

The time snapshots of the diffusivity values acquired via FCS measurements and the measured hydrogel thickness values were both fitted with splines as well as the exponential family in Equation (4). Quadratic splines with break points at 0, 80 and 300 min were used to fit thickness data, while quadratic splines with break points at 0, 150, and 300 min were used to fit diffusivity data. Note that the break points chosen were far fewer than the data points. This was done to obtain smooth curves. Since available MATLAB programs did not allow for choosing the break points to differ from the data points, a MATLAB program was written for this purpose (Barnard, 2018). Exponential fit was accomplished using the lsqnonlin program in MATLAB. Results were input to the mathematical model to obtain release profiles. Figure 5 shows only the results obtained via the spline fit, as both the spline fit and the exponential fit gave similar results (Figure S1). The model-computed release profile of BSA followed the experimental values very well until 90 min, beyond which point the model slightly under-predicted compared to the experimental release values. The model-computed release profile of IgG over-predicted at all time points compared to the experimental release values. However, the trend of the curve of the release profile was similar to the experimental values even with the over-prediction.

Figure 5. Comparison of release profiles of (A) BSA and (B) IgG obtained from the mathematical model and bulk release experiments (n = 9).

Comparison of Time-Varying vs. Averaged Constant Diffusivity and Gel Height

As explained earlier, the constant D- and h-values were calculated as an average of all time points. The constant average diffusivity for BSA was 2.43 × 10⁻⁷ and for IgG was 1.72 × 10⁻⁷ (both in cm²/s) and constant average gel thickness was 1.80 mm. The comparison of release profiles predicted using time-varying D and h as well as constant D and h against experimentally measured release profiles is shown in Figure 6. It can be seen that prediction from the constant model is quantitatively similar to that of the time-varying D and h model for BSA and somewhat different for IgG. In the case of IgG, the constant model better agreed with the experimental release profile.

Figure 6. Release profiles calculated using a variable and a constant D and h model compared to experimental values of (A) BSA and (B) IgG (n = 9).

Prediction of Release Profile Based on Measurements at Three Time Points

To assess the feasibility of accurately predicting release profiles based on a small number of time snapshots, we decided to use the diffusivities and hydrogel thickness obtained at only three points in time. After some trial and error, we found that the release profiles predicted based on fitting the exponential family of curves for D and h measurements taken at 0, 60, and 240 min were very similar to the release profiles predicted by fitting splines for the data from all the time points. This is shown in Figure 7. We also tried fitting splines to the data from these time points. The results were not good and hence are not shown here.

Figure 7. Release profiles calculated using a variable and a three-point D and h model compared to experimental values of (A) BSA and (B) IgG (n = 9).

Discussion

Predicting release from degradable hydrogels has value in many applications such as drug delivery (Lin and Metters, 2006) and tissue engineering (Leddy et al., 2004). Degradable hydrogels are especially suitable for drug delivery because they don't require removal once the drug is delivered. In tissue engineering, they can degrade as native matrix is produced by cells. However, most studies developed to predict solute release from hydrogels focus on non-degradable polymeric matrices, or for simplicity ignore degradation and the associated changes in release (Brazel and Peppas, 2000; Mason et al., 2001; Fu and Soboyejo, 2011).

Here we chose to predict release from degradable hydrogels by focusing on the change in solute diffusivity and gel geometry as a function of time (both are caused by degradation). To monitor diffusivity changes, we could choose from two approaches—a mathematical model to predict the change in diffusivity with time or direct measurement. We measured directly, for better accuracy and because of the difficulties and errors associated with predicting variable solute diffusivity in hydrogels. While many models currently exist for predicting diffusivities in hydrogels, each has unique drawbacks and assumptions that could lead to inaccuracies (Davidson and Deen, 1988; Amsden, 1998a,b, 2002; Masaro and Zhu, 1999; Brazel and Peppas, 2000; Donovan et al., 2016). Also, models typically focus on predicting a constant diffusivity (as opposed to time-varying one as required in our approach) and require information on how degradation changes gel parameters such as mesh size. Our approach of using experimental data on diffusivity allowed us to treat the gel as a “black box.”

To measure diffusivity in situ as a function of time we used FCS. FCS measures the fluctuations emitted from a picoliter illuminated volume containing nanomolar concentrations of fluorescent solutes, typically induced by the motion of the solutes moving in and out of the volume. This powerful technique has been increasingly employed to assess solute diffusion in complex environments such as crowded media (Dauty and Verkman, 2004; Weiss et al., 2004; Zustiak et al., 2011), polymer gels (Michelman-Ribeiro et al., 2004; Stylianopoulos et al., 2010; Zustiak et al., 2010), and biological tissues (Lee et al., 2011). Importantly, being a single molecule technique, FCS can measure diffusivity of solutes at low concentrations with unique sensitivity and precision (Zustiak et al., 2010). The nanomolar working concentrations enable studies of rare and expensive solutes (such as growth factors), which could be cost-prohibitive for other approaches. For comparison, the bulk release experiments used here to validate the model, require mM solute concentrations and are labor-intensive and time-consuming.

However, while we chose FCS for measuring diffusivities, our approach is broader and only requires that a method is used to measure diffusivity as a function of time. New methods to measure diffusivities non-invasively in real time are continuously being developed, but one needs to be aware that each method has advantages and drawbacks. For example, dynamic light scattering (DLS) is well-established for measuring diffusivities of dilute solutes in small volumes of liquid or sol-gel samples (Joosten et al., 1990). However, it requires a solute and polymer of disparate sizes, since both contribute to the scattering signal. Hence, DLS is most suitable for measuring proteins, viruses, micellesk and other microparticles, but not for nanosolutes such as salts and small molecule drugs. DLS has also been shown to overestimate diffusion coefficients (Annunziata et al., 2005). Rayleigh interferometry has been suggested as a technique of superior accuracy, but requires macroscopic concentration gradients, which are not practical for many applications (Annunziata et al., 2005). Nuclear magnetic resonance microscopy (NMR) can measure diffusion of solutes in liquids or gels, but has limits related to spin-spin relaxation times, magnetic gradient strength and eddy currents, which could skew results (Massaro and Zhu, 1999). Further, NMR that allows studying solid samples also requires expensive and dedicated facilities not readily available to most investigators. Fluorescence recovery after photobleaching (FRAP) (Berk et al., 1993) has a great potential, but it is not suitable for low solute concentrations and has the additional drawback of occasional reversible photobleaching, which can skew diffusion results.

For this study we chose to use a PEG hydrogel that would degrade in a matter of hours, to more easily follow changes in solute diffusivity and gel thickness with gel degradation. Hence, we used a hydrogel formed from 4-arm PEG-Ac crosslinked with the custom-synthesized crosslinker PEG-DD1 (Figure 1), which degrades by hydrolysis in ~8 h (Jain et al., 2017a). Experiments were only conducted for 5 h, as the gel lost integrity and was hard to handle after that time. Q_M and mesh size continuously increased, as expected for a bulk degrading hydrogel (Table 1).

To underscore the utility of our approach, we tested two model proteins of varying sizes. Because BSA (66 kDa; R_h = 3.5 nm) is smaller than IgG (150 kDa; R_h = 4.7 nm), we expected BSA to diffuse and release from the gel more rapidly (Zustiak and Leach, 2011; Jain et al., 2017b). Because both proteins are smaller in hydrodynamic radius than the gel mesh (13–17 nm), a diffusion and degradation-controlled release was expected. The changes in these proteins' diffusivities as a function of hydrogel degradation was fitted via a spline or an exponential fit to obtain input to the mathematical model. The data showed a general trend of increasing diffusivity throughout degradation (Figures 4A,B), which was expected based on the increasing mesh size of the gel (Table 1). Further, in Figure 2, we noted that the 5 h timepoint curve was noticeably shifted to the left compared to the 0 h timepoint, indicating a smaller diffusion time (Zustiak et al., 2010; Jain et al., 2017b). There was a more pronounced shift for the smaller BSA compared to IgG; because IgG is larger it will not be able to diffuse as freely as BSA inside the hydrogel.

With the diffusivity and thickness data points, a curve of best fit was generated using MATLAB to be used as the input variables for the mathematical model (Figure 4). The gel thickness did increase ~29%, however the diameter also increased ~33%. We chose not to take the diameter of the gel into account because we are assuming 1D diffusion of the solute in the z-direction. Since the gel can be modeled as a thin cylinder, the gel's two faces are the only diffusing surfaces; as gel thickness increases, the path for diffusion also increases.

To validate the mathematical model, we compared its predictions to release profiles acquired through bulk release experiments. Although the mathematical model's predicted release profile was very similar to the experimental one, there were differences, especially at the later time points. For BSA, model predictions seemed to align well until after ~90 min, at which point the model began to slightly under-predict release (Figure 5A). However, model prediction appeared to be within the error range of the experimental values. With IgG, the model continuously over-predicted release, but the release trends aligned well (Figure 5B). We suggest over-prediction is due to the way diffusivity was measured via bulk release or through FCS. Overall, we noted that diffusivities measured by FCS (tracer diffusivity, Equation 10) were slightly higher than the ones obtained from bulk release experiments (bulk diffusivity, Equation 6). For example, the bulk diffusivity for BSA was 1.71 × 10⁻⁷ cm²/s, corresponding to the initial tracer diffusivity at <100 min (Figures 3B, 4A). However, for the larger IgG, the bulk diffusivity was 0.89 × 10⁻⁷ cm²/s compared to the tracer diffusivity, which started at >1 × 10⁻⁷ cm²/s (Figures 3B, 4B). We expected to see similar, but not exactly the same, diffusivities obtained through FCS and bulk release experiments for several reasons. First, we expected similarities because of the way we performed the FCS measurements. With FCS we measured at multiple gel locations to ascertain that we were not measuring only local diffusivity, but that our measurements represented the entire gel. Further, the illuminated FCS volume had a diameter in the micrometer range, while the hydrogel was nanoporous (mesh size of ~13–17 nm). Hence, even the local diffusivities measured should represent the bulk hydrogel. However, it has also been shown that tracer (as measured by FCS) and bulk diffusivities could differ, with the tracer diffusivity being slightly higher than the bulk (Anderson and Reed, 1976). This could partially explain why, especially for IgG (Figure 4B), the model which relied on FCS data slightly over-predicted release compared to data from bulk release experiments. However, it is important to note that the model overall followed closely the trend predicted by bulk experiments.

We next explored ways to simplify the model or the experimental measurements required for input to the model, to make our approach more broadly accessible to others. Prediction accuracy was compared to the time-varying approach shown in Figure 5 for all cases. In one simplification, the diffusivity and thickness values were averaged across all time points and one constant value was used for each input (Figure 6). Using an averaged constant diffusivity and gel thickness (averaged over the entire process of gel degradation) simplifies the mathematical model used (the PDEs) and does not require a change in variable. For both BSA and IgG, the variable and constant model predictions were similar, with the constant model fitting the experimental data better. We believe the similarity between release predictions by the variable and the constant models may be partially explained by how D and h change as gel degrades. It must be noted that an increase in D increases the rate of release, while an increase in h decreases the rate of release. Since both D and h are increasing with time, the effects of these increases are somewhat but not entirely canceling each other. The constant model may sometimes better predict the release profile because the constant D and h are computed as averages over all time points. Thus, there is greater cancelation of experimental measurement errors, which improves accuracy of the values of D and h used.

In another simplification, we considered decreasing the number of points required to profile the change in diffusivity and gel thickness to only three. Reducing the number of time points simplifies experimental data acquisition. We found that fitting the exponential family of curves to the diffusivity and thickness data at 0, 60, and 240 min was sufficient to predict release profile, as shown in Figure 7. While a trial and error approach was used to select the points, we made sure to include the 0 min time point to represent the gel prior to swelling, a second early time point and one late time point.

Lastly, in fitting continuous smooth curves to the time snapshots to obtain D(t) and h(t), two approaches were explored: quadratic splines and an exponential family of curves given by Equation (4). The advantage of the splines is that the resulting weighted least squares minimization problem is linear, while the corresponding minimization problem for the exponential family is non-linear. The disadvantage of using splines is that they may yield an unrealistic oscillatory curve; this can be avoided with the exponential family. This is particularly useful in fitting the gel thickness h(t), as the derivative h′(t) is not expected to become negative. Moreover, using splines involves several choices which include the degree of the splines, degree of smoothness as well as choice of the break points. Figure S1 demonstrates that results of the spline fit and exponential fit are very close. When the three time points fit is used, the better method of fit is the exponential.

Conclusions

We have developed a mathematical and computational model which uses time snapshots of diffusivity and geometry data as inputs to predict the release profile of proteins embedded in degrading hydrogels. The mathematical model based on Fickian diffusion was described by a 1D PDE with time-varying diffusivity and hydrogel thickness. The time snapshots of diffusivity were measured experimentally by FCS. The overall model was validated for the two model proteins BSA and IgG embedded in PEG hydrogels by comparing our model-predicted release profiles to conventional bulk release profiles obtained through experimentation. Our approach of predicting the release profile bypasses the difficult task of modeling gel degradation to predict time-varying diffusivity and geometry. Further, we showed that the model could be simplified further without loss of accuracy by using either an averaged constant diffusivity and gel thickness (for simpler mathematical modeling) or a three-point measurement (for simpler experimental measurements). The developed approach could be valuable in various applications of degradable hydrogels, such as drug delivery to predict release profiles; or tissue engineering to predict nutrient and metabolite exchange within scaffolds.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

SZ and MR conceived, designed, and supervised the study. SS performed all experimental procedures. EB, BH, and MR developed and implemented the mathematical model. SZ, MR, and SS wrote the manuscript with input from all authors.

Funding

This project was supported by Start-up funds from Saint Louis University awarded to SZ. SS was supported by a Barta Graduate Scholarship awarded from Parks College of Engineering, Aviation and Technology, Saint Louis University.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to thank Dr. Nicola Pozzi for access to and assistance with the FCS instrument.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2019.00410/full#supplementary-material

References