A Diffusion-Reaction Model for Predicting Enzyme-Mediated Dynamic Hydrogel Stiffening

Hydrogels with spatiotemporally tunable mechanical properties have been increasingly employed for studying the impact of tissue mechanics on cell fate processes. These dynamic hydrogels are particularly suitable for recapitulating the temporal stiffening of a tumor microenvironment. To this end, we have reported an enzyme-mediated stiffening hydrogel system where tyrosinase (Tyrase) was used to stiffen orthogonally crosslinked cell-laden hydrogels. Herein, a mathematical model was proposed to describe enzyme diffusion and reaction within a highly swollen gel network, and to elucidate the critical factors affecting the degree of gel stiffening. Briefly, Fick’s second law of diffusion was used to predict enzyme diffusion in a swollen poly(ethylene glycol) (PEG)-peptide hydrogel, whereas the Michaelis–Menten model was employed for estimating the extent of enzyme-mediated secondary crosslinking. To experimentally validate model predictions, we designed a hydrogel system composed of 8-arm PEG-norbornene (PEG8NB) and bis-cysteine containing peptide crosslinker. Hydrogel was crosslinked in a channel slide that permitted one-dimensional diffusion of Tyrase. Model predictions and experimental results suggested that an increasing network crosslinking during stiffening process did not significantly affect enzyme diffusion. Rather, diffusion path length and the time of enzyme incubation were more critical in determining the distribution of Tyrase and the formation of additional crosslinks in the hydrogel network. Finally, we demonstrated that the enzyme-stiffened hydrogels exhibited elastic properties similar to other chemically crosslinked hydrogels. This study provides a better mechanistic understanding regarding the process of enzyme-mediated dynamic stiffening of hydrogels.


Introduction
Hydrogels are hydrophilic and crosslinked water-swollen polymers [1][2][3] particularly suitable for mimicking extracellular matrix (ECM) in human tissues [4,5]. The effects of ECM compositions and degradability on cell fate processes have been extensively studied [6,7]. In recent years, mechanical properties of the ECM are also deemed as a crucial factor regulating tissue regeneration and disease progression [8,9]. As such, hydrogels with spatiotemporally regulated mechanics are increasingly utilized for studying mechanotransduction in healthy and diseased cells. Mechanical properties of a water-swollen hydrogel have been described by Anseth and colleges using Flory-Rehner and rubber elasticity theories [10,11]. In general, hydrogel elasticity, viscosity, and plasticity are characterized via dynamic mechanical analysis or shear rheometry [12,13]. Mechanical properties of a swollen hydrogel are directly related to gel crosslinking density, which are determined by macromer functionality, precursor compositions, polymerization conditions, and degree of gel swelling [14]. Understanding

Correlation of Gel Crosslinking Density, Mesh Size, and Enzyme Diffusivity
Hydrogel crosslinking density has a significant impact on the diffusion of soluble molecules in the network [32]. More specifically, mesh size of a gel network is the primary factor determining diffusivity of any solute in a highly swollen hydrogel. During the process of enzyme-mediated gel stiffening, hydrogel crosslinking density increases not only with time, but may also vary spatially. Consequently, the diffusivity of Tyrase may be impacted by the stiffening network. The diffusivity of any solute in a highly swollen hydrogel can be estimated by the classical Lustig-Peppas relationship, which describes solute diffusivity (i.e., DE gel ) using a correlation of hydrodynamic radius of the solute (RE), mesh size of the network (ξ), and the diffusivity of the solute in a solution (see section 3.3.). Many of these parameters can be obtained from the literature or be determined experimentally. To gain insight into the impact of gel crosslinking on enzyme diffusion in hydrogels, we prepared hydrogels with different macromer compositions that led to varying shear moduli (G' ~0.5 to ~10 kPa). At any given formulation (and modulus), the mass (q) and volumetric (Q) swelling ratio, as well as the corresponding mesh size (ξ), of the resulting hydrogels could be readily determined. These results were then used to establish a correlation between gel stiffness and DE gel (see section 3.3.).
In an earlier work, we characterized the shear moduli of Tyrase-mediated stiffening hydrogels, which varied from ~0.5 to ~5 kPa [28]. This range of gel stiffness is relevant to changes of tissue stiffness during tumor progression [33]. Therefore, we experimentally determined Q and ξ of hydrogels within shear moduli from ~0.5 to ~5 kPa. As polymer content increases, Q was decreased from ~30 to ~18, while ξ was correspondingly decreased from ~19 to ~15 nm ( Figure 2A). It is important to note that this range of mesh size is much larger than the hydrodynamic radius of Tyrase (RE = 4.5 nm) [34]. Next, we estimated DE gel using the Lustig-Peppas relationship [35]. Clearly, diffusivity of Tyrase in solution (DE solution ) was not significantly affected by changes of gel mesh size (i.e., both ξ and Q approach infinity). In a soft gel (e.g., G' ~ 0.5 kPa), DE gel was 3.80 × 10 −11 m 2 /s. In a stiff gel (e.g., G' ~5 kPa), it was decreased slightly to 3.58 × 10 −11 m 2 /s ( Figure 2B). Since DE gel in a stiffer gel is only ~5.8% smaller than that in a softer gel, the gradually increasing gel crosslinking during the stiffening process should not impose a significant diffusion hindrance to Tyrase.

Correlation of Gel Crosslinking Density, Mesh Size, and Enzyme Diffusivity
Hydrogel crosslinking density has a significant impact on the diffusion of soluble molecules in the network [32]. More specifically, mesh size of a gel network is the primary factor determining diffusivity of any solute in a highly swollen hydrogel. During the process of enzyme-mediated gel stiffening, hydrogel crosslinking density increases not only with time, but may also vary spatially. Consequently, the diffusivity of Tyr ase may be impacted by the stiffening network. The diffusivity of any solute in a highly swollen hydrogel can be estimated by the classical Lustig-Peppas relationship, which describes solute diffusivity (i.e., D E gel ) using a correlation of hydrodynamic radius of the solute (R E ), mesh size of the network (ξ), and the diffusivity of the solute in a solution (see Section 3.3). Many of these parameters can be obtained from the literature or be determined experimentally. To gain insight into the impact of gel crosslinking on enzyme diffusion in hydrogels, we prepared hydrogels with different macromer compositions that led to varying shear moduli (G'~0.5 to~10 kPa). At any given formulation (and modulus), the mass (q) and volumetric (Q) swelling ratio, as well as the corresponding mesh size (ξ), of the resulting hydrogels could be readily determined. These results were then used to establish a correlation between gel stiffness and D E gel (see Section 3.3).
In an earlier work, we characterized the shear moduli of Tyr ase -mediated stiffening hydrogels, which varied from~0.5 to~5 kPa [28]. This range of gel stiffness is relevant to changes of tissue stiffness during tumor progression [33]. Therefore, we experimentally determined Q and ξ of hydrogels within shear moduli from~0.5 to~5 kPa. As polymer content increases, Q was decreased from~30 to~18, while ξ was correspondingly decreased from~19 to~15 nm (Figure 2A). It is important to note that this range of mesh size is much larger than the hydrodynamic radius of Tyr ase (R E = 4.5 nm) [34]. Next, we estimated D E gel using the Lustig-Peppas relationship [35]. Clearly, diffusivity of Tyr ase in solution (D E solution ) was not significantly affected by changes of gel mesh size (i.e., both ξ and Q approach infinity). In a soft gel (e.g., G'~0.5 kPa), D E gel was 3.80 × 10 −11 m 2 /s. In a stiff gel (e.g., G'~5 kPa), it was decreased slightly to 3.58 × 10 −11 m 2 /s ( Figure 2B). Since D E gel in a stiffer gel is only~5.8% smaller than that in a softer gel, the gradually increasing gel crosslinking during the stiffening process should not impose a significant diffusion hindrance to Tyr ase .

Prediction of Enzyme Diffusion in Hydrogels with Different Crosslinking Density
Correlations of gel modulus, mesh size, and enzyme diffusivity as shown in Figure 2 have provided critical information regarding the extent to which DE gel was affected by an increasing gel crosslinking density. To establish the premise that the gradually stiffened hydrogel would not impose significant diffusion hindrance to Tyrase, we predicted Tyrase distribution (Equation (1)) within the hydrogels using a constant DE gel , which yields the following equation: here, 3.80 × 10 −11 m 2 /s and 3.58 × 10 −11 m 2 /s were used to represent DE gel in a soft and stiff gel, respectively. If distributions of Tyrase in hydrogels with these two diffusivities show negligible differences within a relevant time scale, it can be safely assumed that the stiffening network only exhibits a minimal hindrance on enzyme transport. Equation (1) was solved numerically using the initial and boundary conditions listed in section 3.3 [36]: The computational results of Equation (2) represent time-and space-dependent Tyrase diffusion into a soft ( Figure 3A) or a stiff gel ( Figure 3B). We plotted the results from 0 to 6 h, a timeline previously used for Tyrase-mediated hydrogel stiffening [27,28]. From the prediction results, it is clear that, regardless of enzyme diffusivity, the entire hydrogel can be equilibrated (i.e., CE/CE0 ≈ 1) with the infiltrating enzyme after only 2 h of diffusion. Furthermore, a symmetrical Tyrase distribution can be clearly seen along the thickness of the gel owing to the bi-directional diffusion condition. While significant variations of Tyrase distribution as a function of time and space are observed in the first 2 h, there is no discernable differences between enzyme diffusion in softer and stiffer hydrogels, suggesting that the stiffening process will not significantly hinder enzyme diffusion in these hydrogels. Finally, a gradient of enzyme concentration can be expected near the surface of the gel within the first 2 h. These predictions have justified that, regardless of gel-network crosslinking density, a period of 6 h is sufficient for Tyrase-mediated hydrogel stiffening.

Prediction of Enzyme Diffusion in Hydrogels with Different Crosslinking Density
Correlations of gel modulus, mesh size, and enzyme diffusivity as shown in Figure 2 have provided critical information regarding the extent to which D E gel was affected by an increasing gel crosslinking density. To establish the premise that the gradually stiffened hydrogel would not impose significant diffusion hindrance to Tyr ase , we predicted Tyr ase distribution (Equation (1)) within the hydrogels using a constant D E gel , which yields the following equation: here, 3.80 × 10 −11 m 2 /s and 3.58 × 10 −11 m 2 /s were used to represent D E gel in a soft and stiff gel, respectively. If distributions of Tyr ase in hydrogels with these two diffusivities show negligible differences within a relevant time scale, it can be safely assumed that the stiffening network only exhibits a minimal hindrance on enzyme transport. Equation (1) was solved numerically using the initial and boundary conditions listed in Section 3.3 [36]: The computational results of Equation (2) represent time-and space-dependent Tyr ase diffusion into a soft ( Figure 3A) or a stiff gel ( Figure 3B). We plotted the results from 0 to 6 h, a timeline previously used for Tyr ase -mediated hydrogel stiffening [27,28]. From the prediction results, it is clear that, regardless of enzyme diffusivity, the entire hydrogel can be equilibrated (i.e., C E /C E0 ≈ 1) with the infiltrating enzyme after only 2 h of diffusion. Furthermore, a symmetrical Tyr ase distribution can be clearly seen along the thickness of the gel owing to the bi-directional diffusion condition. While significant variations of Tyr ase distribution as a function of time and space are observed in the first 2 h, there is no discernable differences between enzyme diffusion in softer and stiffer hydrogels, suggesting that the stiffening process will not significantly hinder enzyme diffusion in these hydrogels. Finally, a gradient of enzyme concentration can be expected near the surface of the gel within the first 2 h. These predictions have justified that, regardless of gel-network crosslinking density, a period of 6 h is sufficient for Tyr ase -mediated hydrogel stiffening.  Note that gel thickness was set at 1 mm with the assumption that Tyrase diffuses symmetrically from the surfaces (x = 0 and x = 1) to the center of the hydrogel (x = 0.5).

Verification of Enzyme Diffusion in Non-Stiffening Hydrogels
In addition to model predictions, we obtained experimental Tyrase diffusion results through imaging Tyrase distribution in a hydrogel strip cast in a channel slide connected by two reservoirs filled with solutions containing the enzyme (CE0). We prepared the hydrogel with high shear moduli (G' ~5 kPa), which represents a stiffened hydrogel network with the most hindrance to solute transport. Figure 4A illustrates the progression of bi-directional Tyrase transport into the thin hydrogel strip (thickness = 1 mm), where Tyrase concentration in the hydrogel (CE) increases as more enzyme molecules infiltrate the hydrogel. Figure 4B shows the experimental Tyrase diffusion profiles at 1, 3, and 6 h (i.e., solid symbols), as well as the diffusion model predictions (i.e., dashed lines). After 1 and 3 h of bi-directional diffusion, both experimental data and the model prediction exhibited symmetrical Tyrase distribution along the diffusion path. After 6 h of diffusion, CE at the middle of the gel reached equilibrium with CE0 in both experiment and model prediction. Since the gels were formed at a higher crosslinking density (i.e., G' ~5 kPa) in this example, a 6-h period should be sufficient for hydrogels with lower stiffness (higher mesh size) to reach equilibrium with CE0. While these assessments have not yet taken enzyme reactions into account, they complement the diffusion model prediction shown in Figure 3. These results are also supported by experimental stiffening results reported previously, in which gel stiffening was completed within 6 h of Tyrase incubation [28].  Note that gel thickness was set at 1 mm with the assumption that Tyr ase diffuses symmetrically from the surfaces (x = 0 and x = 1) to the center of the hydrogel (x = 0.5).

Verification of Enzyme Diffusion in Non-Stiffening Hydrogels
In addition to model predictions, we obtained experimental Tyr ase diffusion results through imaging Tyr ase distribution in a hydrogel strip cast in a channel slide connected by two reservoirs filled with solutions containing the enzyme (C E0 ). We prepared the hydrogel with high shear moduli (G'~5 kPa), which represents a stiffened hydrogel network with the most hindrance to solute transport. Figure 4A illustrates the progression of bi-directional Tyr ase transport into the thin hydrogel strip (thickness = 1 mm), where Tyr ase concentration in the hydrogel (C E ) increases as more enzyme molecules infiltrate the hydrogel. Figure 4B shows the experimental Tyr ase diffusion profiles at 1, 3, and 6 h (i.e., solid symbols), as well as the diffusion model predictions (i.e., dashed lines). After 1 and 3 h of bi-directional diffusion, both experimental data and the model prediction exhibited symmetrical Tyr ase distribution along the diffusion path. After 6 h of diffusion, C E at the middle of the gel reached equilibrium with C E0 in both experiment and model prediction. Since the gels were formed at a higher crosslinking density (i.e., G'~5 kPa) in this example, a 6-h period should be sufficient for hydrogels with lower stiffness (higher mesh size) to reach equilibrium with C E0 . While these assessments have not yet taken enzyme reactions into account, they complement the diffusion model prediction shown in Figure 3. These results are also supported by experimental stiffening results reported previously, in which gel stiffening was completed within 6 h of Tyr ase incubation [28].  Note that gel thickness was set at 1 mm with the assumption that Tyrase diffuses symmetrically from the surfaces (x = 0 and x = 1) to the center of the hydrogel (x = 0.5).

Verification of Enzyme Diffusion in Non-Stiffening Hydrogels
In addition to model predictions, we obtained experimental Tyrase diffusion results through imaging Tyrase distribution in a hydrogel strip cast in a channel slide connected by two reservoirs filled with solutions containing the enzyme (CE0). We prepared the hydrogel with high shear moduli (G' ~5 kPa), which represents a stiffened hydrogel network with the most hindrance to solute transport. Figure 4A illustrates the progression of bi-directional Tyrase transport into the thin hydrogel strip (thickness = 1 mm), where Tyrase concentration in the hydrogel (CE) increases as more enzyme molecules infiltrate the hydrogel. Figure 4B shows the experimental Tyrase diffusion profiles at 1, 3, and 6 h (i.e., solid symbols), as well as the diffusion model predictions (i.e., dashed lines). After 1 and 3 h of bi-directional diffusion, both experimental data and the model prediction exhibited symmetrical Tyrase distribution along the diffusion path. After 6 h of diffusion, CE at the middle of the gel reached equilibrium with CE0 in both experiment and model prediction. Since the gels were formed at a higher crosslinking density (i.e., G' ~5 kPa) in this example, a 6-h period should be sufficient for hydrogels with lower stiffness (higher mesh size) to reach equilibrium with CE0. While these assessments have not yet taken enzyme reactions into account, they complement the diffusion model prediction shown in Figure 3. These results are also supported by experimental stiffening results reported previously, in which gel stiffening was completed within 6 h of Tyrase incubation [28].

Effect of Enzyme Concentration on Reaction Velocity
In addition to predicting enzyme diffusion in gels with different crosslinking densities, we investigated catalytic reactions of Tyr ase using different phenolic substrates. Tyr ase is known to catalyze oxidization of Tyr, tyramine, and other phenolic derivatives (e.g., hydroxyphenylacetic acid) into DOPA, DOPAquinone, and finally to DOPA dimers [34]. The catalytic reaction of Tyr ase involves several steps, including a monophenol cycle, a diphenol cycle, and substrate inhibition ( Figure 5A) [37]. First, the non-activated deoxy-tyrosinase (Tyr ase Deoxy ) binds with the oxygen molecule (O 2 ) to generate an activated form of tyrosinase (Tyr ase Oxy ). In the presence of L-Tyr, Tyr ase Oxy initiates the monophenol cycle that produces DOPA. Since DOPA and L-Tyr are both substrate of the Tyr ase Oxy , the excess L-Tyr can react with both Tyr ase Oxy and the DOPA-Tyr ase Oxy complexes, and thus reaction may move toward diphenol cycle or substrate inhibition [38]. To gain more insight into these reaction kinetics, we monitored the concentration of dissolved O 2 , as its disappearance is the first step in the activation of Tyr ase . As expected, concentration of dissolved O 2 decreased only after the addition of Tyr ase (at the 2-min mark) to the L-Tyr containing solution ( Figure 5B). O 2 content dropped rapidly, except for the lowest enzyme concentration used (i.e., 0.3 µM). Furthermore, dissolved O 2 was completely depleted within 5, 3.5, and 2 min when the solution was added to 1.5, 2.25, and 3 µM Tyr ase , respectively. At 0.3 µM Tyr ase , only~3% decrease in dissolved O 2 was detected after 6 min of enzyme addition, potentially because the rate of oxygen consumption by the small amount of enzyme was much slower than its replenishment from the air.

Effect of Enzyme Concentration on Reaction Velocity
In addition to predicting enzyme diffusion in gels with different crosslinking densities, we investigated catalytic reactions of Tyrase using different phenolic substrates. Tyrase is known to catalyze oxidization of Tyr, tyramine, and other phenolic derivatives (e.g., hydroxyphenylacetic acid) into DOPA, DOPAquinone, and finally to DOPA dimers [34]. The catalytic reaction of Tyrase involves several steps, including a monophenol cycle, a diphenol cycle, and substrate inhibition ( Figure 5A) [37]. First, the non-activated deoxy-tyrosinase (Tyrase Deoxy ) binds with the oxygen molecule (O2) to generate an activated form of tyrosinase (Tyrase Oxy ). In the presence of L-Tyr, Tyrase Oxy initiates the monophenol cycle that produces DOPA. Since DOPA and L-Tyr are both substrate of the Tyrase Oxy , the excess L-Tyr can react with both Tyrase Oxy and the DOPA-Tyrase Oxy complexes, and thus reaction may move toward diphenol cycle or substrate inhibition [38]. To gain more insight into these reaction kinetics, we monitored the concentration of dissolved O2, as its disappearance is the first step in the activation of Tyrase. As expected, concentration of dissolved O2 decreased only after the addition of Tyrase (at the 2-min mark) to the L-Tyr containing solution ( Figure 5B). O2 content dropped rapidly, except for the lowest enzyme concentration used (i.e., 0.3 µM). Furthermore, dissolved O2 was completely depleted within 5, 3.5, and 2 min when the solution was added to 1.5, 2.25, and 3 µM Tyrase, respectively. At 0.3 µM Tyrase, only ~3% decrease in dissolved O2 was detected after 6 min of enzyme addition, potentially because the rate of oxygen consumption by the small amount of enzyme was much slower than its replenishment from the air.  Figure 5B suggests that a faster Tyrase reaction was accompanied by a rapid consumption of dissolved oxygen. However, the consumption of dissolved oxygen represents only the first step in the Tyrase reaction cycle (i.e., from Tyrase Deoxy to Tyrase Oxy , Figure 5A). In order to understand the kinetics of subsequent reaction steps, it is necessary to determine the amount of actual product formation. To  Figure 5B suggests that a faster Tyr ase reaction was accompanied by a rapid consumption of dissolved oxygen. However, the consumption of dissolved oxygen represents only the first step in the Tyr ase reaction cycle (i.e., from Tyr ase Deoxy to Tyr ase Oxy , Figure 5A). In order to understand the kinetics of subsequent reaction steps, it is necessary to determine the amount of actual product formation. To this end, we utilized MBTH assay to monitor the production of DOPA [39,40]. In principle, MBTH reacts with oxidized substrates through both monophenol and diphenol cycles and produces a visible complex with a pink color [41]. Although higher DOPA contents were detected at higher doses of Tyr ase (Figure 5C), all reactions (except for 0 µM Tyr ase ) were not completed within 6 min of testing, suggesting that the catalytic step of Tyr ase /L-Tyr reaction was slower than the rate of oxygen binding to Tyr ase Deoxy . Additionally, the kinetics of DOPA production appeared to be in a linear relationship with respect to time. After plotting reaction velocity as a function of enzyme concentration, a linear correlation was established ( Figure 5D). This linear relationship might be due to a relatively high substrate concentration (i.e., 10 mM) when compared with the high binding affinity between L-Tyr and Tyr ase . In the case of a much smaller K M compared with C S , K M can be omitted in either the regular or the modified Michaelis-Menten equation (see Section 3.5), which yields the following equations: as shown in Equations (3) and (4), omitting K M results in a linear correlation between V p and C E regardless of the status of substrate inhibition, which can be characterized by a decreasing reaction rate at high substrate concentrations [42]. Using the general or modified Michaelis-Menten equation, we obtained k cat and K M for Tyr ase -mediated reactions (Table 1). Clearly, all three substrates exhibited at least 10-fold lower K M than the substrate concentration used in the experiments (i.e., 10 mM, Figure 5), thus justifying the omission of K M in Equations (3) and (4).

Effect of Substrate Concentration of Enzymatic Reaction
As mentioned earlier, utilization of Tyr for Tyr ase may exhibit substrate inhibition that reduces catalytic activity of Tyr ase . To evaluate whether such an effect exists when a peptide substrate is used, we treated a model peptide CYGGGYC with Tyr ase and used other substrates as controls (e.g., L-Tyr or L-DOPA). As shown in Figure 6A, Tyr ase exhibited the highest reactivity for L-DOPA (V P = 8.6 µM/s) among all substrates. The reaction rates for L-Tyr and CYGGGYC were 0.88 and 0.70 µM/s, respectively. Tyr ase /L-DOPA reaction appeared to proceed through the diphenol cycle without discernable substrate inhibition, as the reaction velocity reached a plateau value at high substrate concentration. Furthermore, we noticed a slightly lower maximum reaction rate at a higher concentration (2-10 mM) of L-Tyr, which was indicative of substrate inhibition ( Figure 6A). Therefore, the modified Michaelis-Menten equation (see Section 3.5) was utilized to obtain kinetic constants for Tyr ase /L-Tyr reaction. Interestingly, there was no significant substrate inhibition when tyrosine-containing peptide (i.e., CYGGGYC) was used as substrate for Tyr ase . It is likely that the peptidyl tyrosine residues exhibited different affinity (K M ) for Tyr ase . Indeed, Marumo et al. have investigated the effect of peptide sequence on Tyr ase reaction efficiency [43]. Compared to Tyr, some tyrosine-containing peptide sequences (e.g., Gly-Tyr-Gly and Lys-Glu-Thr-Tyr-Ser-Lys) have a higher DOPA conversion ratio when reacting with Tyr ase due to a higher binding affinity between Tyr ase and oxidized Tyr. Upon fitting the reaction velocity data with the general or modified Michaelis-Menten model (see Section 3.5), we found that, compared with soluble L-Tyr, peptidyl Tyr exhibited a higher binding affinity (i.e., lower K M ) with Tyr ase (Table 1).
DOPA dimer crosslinking). At any given enzyme and substrate concentration, the time needed to convert all substrates into products can be approximated through dividing the initial substrate content (CS0) by the reaction velocity (VP) and catalytic efficiency (kcat). Regardless of substrate concentration, there is a hyperbolic relationship between reaction time and enzyme concentration. Naturally, at a higher enzyme concentration, the time needed to convert all substrates (10 mM) into products would be much faster than using a lower enzyme concentration. For example, at 3 µM of Tyrase, it would take about 2 h to convert all 10 mM of substrates into products. Without considering the enzyme diffusion, Figure 6B provides a first pass estimation of the time needed to achieve stiffening at any given enzyme and substrate concentration.

Numerical Simulation of Diffusion-Reaction in Hydrogel
We have separately characterized/analyzed Tyrase diffusion (Figures 2-4) in hydrogels and studied enzymatic reactions of Tyrase with tyrosine-containing peptide crosslinker (Figures 5 and 6). Next, we considered both enzyme diffusion and reaction to obtain profiles of peptide substrate consumption and product formation within the gel network over time and space. First, we generated computational data using Tyrase diffusion profile CE(x,t), which served as input for the Michaelis-Menten equation. With experimentally obtained kcat and KM for the peptide linker CYGGGYC (Table  1), we employed the Lambert W function to numerically solve the Michaelis-Menten equation (see section 3.5.) [44,45]. First, we derived the substrate consumption rate (Equation (5)) from the general Michaelis-Menten equation: where the space-and time-dependent CS(x,t) can be expressed as: The space-and time-dependent product (i.e., oxidized tyrosine or DOPA) concentration can be expressed using the following equation, since it is equivalent to the amount of substrate (i.e., tyrosine) consumed: Another key aspect in designing an enzymatic stiffening hydrogel is to understand how fast the enzyme fully converts the substrates (e.g., peptidyl tyrosine residues) into products (i.e., additional DOPA dimer crosslinking). At any given enzyme and substrate concentration, the time needed to convert all substrates into products can be approximated through dividing the initial substrate content (C S0 ) by the reaction velocity (V P ) and catalytic efficiency (k cat ). Regardless of substrate concentration, there is a hyperbolic relationship between reaction time and enzyme concentration. Naturally, at a higher enzyme concentration, the time needed to convert all substrates (10 mM) into products would be much faster than using a lower enzyme concentration. For example, at 3 µM of Tyr ase , it would take about 2 h to convert all 10 mM of substrates into products. Without considering the enzyme diffusion, Figure 6B provides a first pass estimation of the time needed to achieve stiffening at any given enzyme and substrate concentration.

Numerical Simulation of Diffusion-Reaction in Hydrogel
We have separately characterized/analyzed Tyr ase diffusion (Figures 2-4) in hydrogels and studied enzymatic reactions of Tyr ase with tyrosine-containing peptide crosslinker (Figures 5 and 6). Next, we considered both enzyme diffusion and reaction to obtain profiles of peptide substrate consumption and product formation within the gel network over time and space. First, we generated computational data using Tyr ase diffusion profile C E (x,t), which served as input for the Michaelis-Menten equation. With experimentally obtained k cat and K M for the peptide linker CYGGGYC (Table 1), we employed the Lambert W function to numerically solve the Michaelis-Menten equation (see Section 3.5) [44,45]. First, we derived the substrate consumption rate (Equation (5)) from the general Michaelis-Menten equation: where the space-and time-dependent C S (x,t) can be expressed as: The space-and time-dependent product (i.e., oxidized tyrosine or DOPA) concentration can be expressed using the following equation, since it is equivalent to the amount of substrate (i.e., tyrosine) consumed: We plotted solution of Equation (7) in Figure 7, which describes the formation of DOPA dimers owing to the Tyr ase -mediated reaction within the gel network. The prediction results demonstrate that the majority of network-immobilized tyrosine residues will be converted to product within the first 6 to 8 h of Tyr ase (3 µM) diffusion/reaction. Furthermore, comparing simulation results shown in Figure 3 (i.e., enzyme distribution in hydrogel) with Figure 7 (product formation in hydrogel), it is clear that, while enzyme diffusion may not be critically affected during the stiffening process (Figure 3), the rate of product formation lags behind enzyme diffusion. Hence, it is critically important to take both diffusion and reaction into account when predicting the degree of additional crosslink formation.
We plotted solution of Equation (7) in Figure 7, which describes the formation of DOPA dimers owing to the Tyrase-mediated reaction within the gel network. The prediction results demonstrate that the majority of network-immobilized tyrosine residues will be converted to product within the first 6 to 8 h of Tyrase (3 µM) diffusion/reaction. Furthermore, comparing simulation results shown in Figure 3 (i.e., enzyme distribution in hydrogel) with Figure 7 (product formation in hydrogel), it is clear that, while enzyme diffusion may not be critically affected during the stiffening process ( Figure  3), the rate of product formation lags behind enzyme diffusion. Hence, it is critically important to take both diffusion and reaction into account when predicting the degree of additional crosslink formation.

Correlation of Hydrogel Mechanical Property and Its Microstructure
In our previous experimental results, we have shown that a period of 6 h is sufficient to induce enzymatic hydrogel stiffening [27,28]. Further increasing the enzyme incubation time to 12-48 h only marginally increased gel stiffness. Prior experimental observations were largely in agreement with the diffusion-reaction modeling results (Figure 7). To gain a deeper understanding of the crosslinking of enzyme-stiffened hydrogels, we prepared additional PEG8NB-CYGGGYC hydrogels (G0 ' ~1 kPa) and performed enzyme-induced stiffening using different concentrations of Tyrase (0-3 µM). After 6 h of stiffening (and overnight washing to remove residue enzyme), gel moduli and swelling ratios were measured. As shown in Figure 8A, Tyrase treatment led to increased gel shear moduli (from ~1 to ~4 kPa) and decreased volumetric swelling ratio (from ~36 to ~12), a correlation commonly observed in chemically crosslinked hydrogels [14,46]. Next, we examined whether the correlation of G' and the polymer volume fraction (ν2,s, a reciprocal of Q) of the enzyme-stiffened hydrogels could be described by the classical rubber elasticity theory, which predicts a linear dependency between G' and ν2,s [10,47]. As shown in Figure 8B, the correlation of experimentally obtained G' and ν2,s follows a power law with an exponent of 1.95. This value is significantly higher than the linear dependence (i.e., exponent = 1) predicted by the theory of rubber elasticity, which suggests the presence of

Correlation of Hydrogel Mechanical Property and Its Microstructure
In our previous experimental results, we have shown that a period of 6 h is sufficient to induce enzymatic hydrogel stiffening [27,28]. Further increasing the enzyme incubation time to 12-48 h only marginally increased gel stiffness. Prior experimental observations were largely in agreement with the diffusion-reaction modeling results (Figure 7). To gain a deeper understanding of the crosslinking of enzyme-stiffened hydrogels, we prepared additional PEG8NB-CYGGGYC hydrogels (G 0 '~1 kPa) and performed enzyme-induced stiffening using different concentrations of Tyr ase (0-3 µM). After 6 h of stiffening (and overnight washing to remove residue enzyme), gel moduli and swelling ratios were measured. As shown in Figure 8A, Tyr ase treatment led to increased gel shear moduli (from~1 tõ 4 kPa) and decreased volumetric swelling ratio (from~36 to~12), a correlation commonly observed in chemically crosslinked hydrogels [14,46]. Next, we examined whether the correlation of G' and the polymer volume fraction (ν 2,s , a reciprocal of Q) of the enzyme-stiffened hydrogels could be described by the classical rubber elasticity theory, which predicts a linear dependency between G' and ν 2,s [10,47]. As shown in Figure 8B, the correlation of experimentally obtained G' and ν 2,s follows a power law with an exponent of 1.95. This value is significantly higher than the linear dependence (i.e., exponent = 1) predicted by the theory of rubber elasticity, which suggests the presence of network non-ideality caused either by ineffective initial crosslinking (performed with low macromer concentration) [16] and/or enzymatic stiffening. However, it should be noted that this degree of network non-ideality (i.e., the exponent in the power law) is similar to other chemically crosslinked hydrogels [48]. Therefore, it can be concluded that the structure-function relationships of enzyme-stiffened hydrogels are similar to other chemically crosslinked hydrogels.
Gels 2019, 5, x FOR PEER REVIEW 10 of 16 network non-ideality caused either by ineffective initial crosslinking (performed with low macromer concentration) [16] and/or enzymatic stiffening. However, it should be noted that this degree of network non-ideality (i.e., the exponent in the power law) is similar to other chemically crosslinked hydrogels [48]. Therefore, it can be concluded that the structure-function relationships of enzymestiffened hydrogels are similar to other chemically crosslinked hydrogels. In summary, we have predicted and validated enzyme diffusion in gels with varying stiffness (Figures 2-4). We have also monitored Tyrase-mediated reactions and predicted the degree of gel stiffening from the perspective of product formation under varying enzyme and substrate concentrations (Figures 5 and 6). The modeling and experimental results have demonstrated that while a period of 3 h is sufficient for the enzyme to equilibrate the entire gel, it takes at least another 3 h for the reactions to be completed (Figure 7) [28]. Furthermore, the enzyme-stiffened hydrogels exhibit physical properties (e.g., swelling and modulus) similar to other chemically crosslinked hydrogels ( Figure 8). Collectively, the information acquired from this study should be highly useful in designing enzyme-mediated dynamic hydrogel system for fundamental material science and tissue engineering applications.

Macromer Preparation and Peptide Synthesis
PEG8NB and lithium phenyl-2,4,6-trimethylbenzoylphosphinate (LAP) were synthesized as described elsewhere [16,49,50]. Peptides (e.g., CYGGGYC) were synthesized via standard Fmoc coupling chemistry on an automated, microwave-assisted peptide synthesizer (Liberty 1, CEM, Matthews, NC, USA). The crude products were cleaved in a trifluoroacetic acid (TFA) cleavage cocktail composed of 7.6 mL trifluoroacetic acid (TFA), 0.2 mL triisopropylsilane (TIPS), 0.2 mL In summary, we have predicted and validated enzyme diffusion in gels with varying stiffness (Figures 2-4). We have also monitored Tyr ase -mediated reactions and predicted the degree of gel stiffening from the perspective of product formation under varying enzyme and substrate concentrations (Figures 5 and 6). The modeling and experimental results have demonstrated that while a period of 3 h is sufficient for the enzyme to equilibrate the entire gel, it takes at least another 3 h for the reactions to be completed (Figure 7) [28]. Furthermore, the enzyme-stiffened hydrogels exhibit physical properties (e.g., swelling and modulus) similar to other chemically crosslinked hydrogels ( Figure 8). Collectively, the information acquired from this study should be highly useful in designing enzyme-mediated dynamic hydrogel system for fundamental material science and tissue engineering applications.

Macromer Preparation and Peptide Synthesis
PEG8NB and lithium phenyl-2,4,6-trimethylbenzoylphosphinate (LAP) were synthesized as described elsewhere [16,49,50]. Peptides (e.g., CYGGGYC) were synthesized via standard Fmoc coupling chemistry on an automated, microwave-assisted peptide synthesizer (Liberty 1, CEM, Matthews, NC, USA). The crude products were cleaved in a trifluoroacetic acid (TFA) cleavage cocktail composed of 7.6 mL trifluoroacetic acid (TFA), 0.2 mL triisopropylsilane (TIPS), 0.2 mL distilled water, and 400 mg phenol. The cleaved and dried peptides were purified by reverse-phase HPLC (PerkinElmer Flexar system) using 95%/5% (v/v) water/acetonitrile (ACN) with a trace amount (0.1 vol %) of TFA as the starting mobile phase. A linear gradient of ACN was employed to separate the products through a semi-prep peptide C18 column (5 mL/min). The separated products were monitored with a UV/vis detector and the purified peptides were characterized with liquid chromatography coupled with mass spectrometry (1200 series LC/MS system, Santa Clara, CA, USA).

Modeling of Enzyme Diffusion into Hydrogels
We used Fick's second law of diffusion (Equation (10)) to estimate distribution of Tyr ase in hydrogel over space and time: Note that the diffusivity of the enzyme in hydrogel (D E gel ) negatively scales with hydrogel crosslinking density and positively correlates to hydrogel mesh size (ξ). As enzyme infiltrates the hydrogel, it catalyzes tyrosine residues into DOPA dimers, resulting in an increase in hydrogel crosslinking density (i.e., smaller mesh size) that may lead to lower diffusivity. Hence, D E gel may be dependent on the location of the enzyme in the hydrogel. In this mathematical model, we assume that enzyme diffusion proceeds bi-directionally from the top and bottom (i.e., along the x-axis) of the thin hydrogel with a thickness of h (i.e., neglecting diffusion from the gel edge). Initially, the entire hydrogel (i.e., 0 < x < h) is free of enzyme, which gives the following initial condition: we also assume that enzyme concentration in the solution remains unchanged (C E0 ), which yields the following boundary conditions: since enzyme diffusion proceeds bi-directionally, at any given time there is no flux at the center of the hydrogel (i.e., x = h/2), which yields the following boundary condition: As mentioned earlier, diffusivity of any solute in a crosslinked hydrogel is affected by gel mesh size. Hence, it is necessary to determine whether the additional DOPA dimer crosslinks significantly affects Tyr ase transport in the stiffened gel. In this regard, the Lustig-Peppas estimation of solute diffusivity in a highly swollen gel can be used to establish a correlation between Tyr ase diffusivity and hydrogel mesh size [35].
where D E sol is the diffusivity of enzyme in solution (5.05 × 10 −10 m 2 /s for Tyr ase ), R E is the hydrodynamic radius of the enzyme (4.5 nm for Tyr ase ) [51], and Y is the critical volume required for a successful translational movement of the substrate relative to the average free volume of a water molecule (1 for PEG-based gels) [35]. Note that this equation is used only when gel mesh size is larger than hydrodynamic radius of the soluble molecule (i.e., R E /ξ < 1) because no diffusion is possible when R E is larger than ξ. Clearly, when R E approaches ξ, D E gel becomes much smaller than D E solution .
Therefore, it is critical to determine the relative sizes of R E and ξ even for highly swollen PEG-based  (15) where the volumetric swelling ratio (Q) of the hydrogel is determined using mass swelling ratio (q) and the density (ρ) of PEG (1.087 g/cm 3 at 37 • C) and water (0.994 g/cm 3 at 37 • C): once Q is obtained, gel mesh size (ξ) can be derived using the following equation: where is the average bond length (1.47 Å) in the backbone of an ethylene glycol subunit (i.e., -CH 2 -CH 2 -O-), 3 is the total number of bonds in a PEG repeat subunit, C n is the Flory characteristic ratio (4 for PEG-based hydrogel), and M n is the number for average molecular weight of PEG8NB (20 kDa) [1]. To obtain the average molecular weight between crosslinks (M c ) within a step-growth hydrogel, the following equation can be used: where MW A and MW B represent molecular weights of the two macromer crosslinkers PEG8NB and CYGGGYC, respectively. f A and f B are the functionality of the macromer and peptide crosslinker (i.e., 8 and 2). Numerical model predictions were programmed and executed via Grapher (Arizona Software) using a spatial step size (∆x) of 1 µm and a temporal step size (∆t) of 1 s.

Characterization of Oxygen Consumption
During the oxidation of Tyr, Tyr ase exhibits four distinct oxidation states (deoxy-, oxy-, met-, and deact-Tyr ase ), and the oxygen molecule (O 2 ) is a key component initiating the overall catalytic processes [39]. The association of O 2 and Tyr ase forms oxy-Tyr ase that oxidizes Tyr (or other phenolic precursors) into DOPAquinone. Hence, the consumption of soluble O 2 can be used to gauge the extent of Tyr ase -catalyzed DOPA formation. To this end, we measured soluble O 2 contents in a 2 mL microtube which contained L-Tyr (10 mM) and Tyr ase (0.1, 1.5, and 3 µM) by a portable fiber optic oxygen probe and meter (PreSens Microx 4, Regensburg, Germany). The needle-type probe was extended into the center of the solution. Note that the O 2 quantifications were started 2 min prior to Tyr ase addition.

Tyr ase -Mediated Reaction Kinetics
The kinetics of Tyr ase -mediated reactions were modeled using standard Michaelis-Menten equation [52]: where C S and C P are the concentration of substrate and product, respectively; V S and V P are the velocity of substrate consumption and product formation; and V max is the maximum reaction velocity, which is equivalent to enzyme concentration (C E ) multiplied by the turnover number (k cat ).
In the event of substrate inhibition, a modified Michaelis-Menten equation was employed [53]: here, K i is the kinetic constant for substrate inhibition [38]. We utilized 3-Methyl-2-benzothiazolinone hydrazone (MBTH, 2 mM in pH 6.5 PBS) as an indicator to characterize enzymatic reactions between Tyr ase (0.6 µM) and various substrates (0-10 mM), including L-Tyr, L-DOPA, and CYGGGYC. MBTH is a hydrazone-based compound that complexes with DOPAquinone to produce a pink pigment, which can be quantified via measuring absorbance at 475 nm [40,41]. Note that cysteine residue also reacts with MBTH. Hence, prior to the MBTH assay, all peptides were conjugated with linear norbornene-functionalized PEG (PEGdNB) through thiol-norbornene photo-click reaction (1 mM LAP, 365 nm light at 5 mW/cm 2 , 2 min). The enzymatic reactions were allowed to proceed for 20 min and the increases in 475 nm absorbance were monitored with a microplate reader (Biotek, Synergy HTX, Winooski, VT, USA). At any given Cs, reaction velocity (V P ) can be obtained by calculating the slope of absorbance vs. the time curve. Next, V Max and Cs values were fitted to the general or modified Michaelis-Menten model (i.e., Equation (19) or (20)) to acquire the k cat and Michaelis constant (K M ). Additionally, to understand the correlation between C E and V P , we also treated the substrate (L-Tyr, 10 mM) with varying Tyr ase concentrations (0-3 µM). The reactions were monitored via the MBTH assay as described above.

Fabrication and Characterization of the Step-Growth PEG-Peptide Hydrogels
PEG-peptide hydrogels were prepared by coupling norbornene moieties of PEG8NB and the terminal cysteine moieties of CYGGGYC via thiol-norbornene photopolymerization. Hydrogels were prepared with different macromer contents but with the stoichiometric ratio of thiol to norbornene maintained at 1 to minimize network defects. In brief, aliquots of precursor solutions (45 µL) were deposited between two glass slides separated by silicon spacers (1 mm). Hydrogel crosslinking was initiated by 365 nm light (5 mW/cm 2 ) exposure for 2 min. Following the initial photopolymerization, gels were maintained in PBS at 37 • C for 24 h prior to characterization or Tyr ase -mediated stiffening. Shear moduli (strain-sweep mode) of the pre-and post-stiffened hydrogels were characterized using a digital rheometer (Bohlin CVO 100, Malvern Instruments, Malvern, Worcestershire, UK) equipped with an 8-mm diameter parallel plate geometry. Measurements of gel shear moduli were acquired from the average of the linear region on the modulus-strain curve (0.1-5% strain) with oscillation frequency set at 1 Hz.

Statistical Analysis
All experiments were repeated three times independently with four samples per condition. Experimental results were reported as Mean ± SEM with a sample size of at least three (n = 3). The data were analyzed by one-way or two-way ANOVA with GraphPad Prism 8 software. Single, double, and triple asterisks represent p < 0.05, 0.001, and 0.0001, respectively, and p < 0.05 is considered statistically significant.