Mechanical Behavior of 3D Printed Poly(ethylene glycol) Diacrylate Hydrogels in Hydrated Conditions Investigated Using Atomic Force Microscopy

Three-dimensional (3D) printed hydrogels fabricated using light processing techniques are poised to replace conventional processing methods used in tissue engineering and organ-on-chip devices. An intrinsic potential problem remains related to structural heterogeneity translated in the degree of cross-linking of the printed layers. Poly(ethylene glycol) diacrylate (PEGDA) hydrogels were used to fabricate both 3D printed multilayer and control monolithic samples, which were then analyzed using atomic force microscopy (AFM) to assess their nanomechanical properties. The fabrication of the hydrogel samples involved layer-by-layer (LbL) projection lithography and bulk cross-linking processes. We evaluated the nanomechanical properties of both hydrogel types in a hydrated environment using the elastic modulus (E) as a measure to gain insight into their mechanical properties. We observed that E increases by 4-fold from 2.8 to 11.9 kPa transitioning from bottom to the top of a single printed layer in a multilayer sample. Such variations could not be seen in control monolithic sample. The variation within the printed layers is ascribed to heterogeneities caused by the photo-cross-linking process. This behavior was rationalized by spatial variation of the polymer cross-link density related to variations of light absorption within the layers attributed to spatial decay of light intensity during the photo-cross-linking process. More importantly, we observed a significant 44% increase in E, from 9.1 to 13.1 kPa, as the indentation advanced from the bottom to the top of the multilayer sample. This finding implies that mechanical heterogeneity is present throughout the entire structure, rather than being limited to each layer individually. These findings are critical for design, fabrication, and application engineers intending to use 3D printed multilayer PEGDA hydrogels for in vitro tissue engineering and organ-on-chip devices.


Analysis of force-distance (P-h) curves and contact point identification
The samples were initially analyzed by extracting force-distance (P-h) curves from NanoScope Analysis software (version 1.85) and JPK SPM Data Processing 8.0.13, and the elastic modulus (E) was estimated using both Hertz model and Oliver and Phar method by analyzing the extend (loading) and retract (unloading) curves separately. For using the Hertz model, to calculate the elastic modulus (E), a number of assumptions were made: (i) indentations are frictionless and for small strains, (ii) indentation depth is less than 1/3 of the probe radius, (iii) deformation is linearly elastic, (iv) for small fast indentations the viscoelasticity can be negligible. 1,2 JPK SPM software was used for batch processing of the (P-h) curves by following the steps: (i) smoothen the data using the Gaussian function, (ii) subtract the baseline offset by fitting range of 75% to 95%, (iii) correct for any tilt, (iv) identify contact point by adjusting the x-offset where the contact point is assumed to be the first point at which the recorded force is in a positive value, (v) Hertz model for spherical indenter was used to the extend curve.

S-4
The indentation depth (δ) or the tip sample separation is measured from the identified contact point to the maximum displacement. The PEGDA Poisson's ratio is assumed to be 0. 5. 4 Analysis of (P-h) curves of AFM has always been challenging especially in locating the contact point to calculate the indentation depth which is a key parameter to measure E. 2,3,5 For hard samples, it is a straightforward step where it can be easily determined by using the extend curve and looking at the clear transition point. However, for soft samples, deformation due to long range forces happens before tip-sample contact, making the contact transition difficult to define and to differentiate between the cantilever movement and sample indentation depth. This can in turn can lead to large errors and deriving over/underestimated E for the sample. 6 In this study, we used a simple procedure where we defined the contact point as a point at which the recorded force crosses the y-axis for the first time. It should be noted that a few steps must be followed before identifying the contact point which is explained in the Methods and Materials section. Such an approach, even though not ideal, if used to analyze all (P-h) curves consistently can provide relatively good comparative results to understand the variation in the mechanical properties as a function of the build (printed) direction.
In Oliver and Pharr method the projected contact area of the indenter is determined based on the indentation depth. For fitting the (p-h) curves an external MATLAB script is used to fit the model to the retract curve where E is calculated following the retraction of the tip from the sample. S-7

Comparison of elastic modulus maps based on two different fitting ranges (Hertz model) and a fitting based on Oliver and Phar method
It is a common problem with soft materials to establish a contact point and to fit the unloading slope in the (P-h) curve ( Figure S2a). To combat this issue, two different distance ranges, from 0.0 to 1.5 µm and from 0.5 to 2.0 µm, were chosen to fit the data for the slope to obtain E and a comparison of these results are shown in Figure S4a and b. These results showed that even though the magnitude of E increased when the higher range was used, the overall pattern and trend did not change. Therefore, it can be suggested that the results are comparable and point to a distinct and clear phenomenon concerning the change in E in a single layer of the multilayer sample. The lower fitting range was used for the remaining measurements. In terms of the difference between Hertz model ( Figure S4a and b) and Oliver and Phar method ( Figure S4c) it can be seen that the results are reporting a similar pattern; however, the magnitude of E has increased by 32% and 92% for both monolithic and multilayer samples, respectively when Oliver and Phar method was used.
To measure E using the Oliver and Phar method, the approach involves identifying a suitable multiplier that can be utilized in the script to determine the correct contact point from the differential comparison between the first point and the subsequent points on the extension curve.
This process involves examining each individual curve to ensure that the fitting is correct and the contact point has been accurately identified. However, this task is tedious and requires significant data processing time. As a result, we chose to use the Hertz model as the primary model for the study based on the assumptions outlined in the supporting information.

Explanations of effect of surface topography on E measurements
In this study, two important arguments are addressed about the outer surface topography of the sample and the contact point. First, looking at Figure 2aii, it could be argued that the peak to valley (trough) depth in the multilayer sample is higher and therefore, the measured E map may be influenced by how contact is established. In the Hertz model, eq S1 and S2, an ideal contact is assumed when the AFM tip is normal to the material's surface and the contact area follows the shape of the tip. Figure S6 postulates scenarios of the type of contact between the AFM tip and the sample's surface with red arrows highlighting the first contact point between the two objects.
Given the size of the AFM tip (radius of 3.46 µm) and the printed layer thickness (20 µm), Figure   S6ai shows how the contact may appear at the two points if the tip was approaching the interface between two consecutive layers, the trough. Conversely, if the tip is away from the interface, the contact scenario may differ and may occur at a single point ( Figure S6aii). Two additional scenarios may also occur for the peak area of the printed layer where the tip may land exactly normal to the surface, the ideal situation, or on the side of the tip where the tip may slip on the surface rather than creating a well-established contact with the sample. Looking at Figure S5, where the exact position of the contact can be determined from the white dots on the apparent topography image, the scenarios mentioned above are plausible and cannot be ruled out.
S-11 Figure S6. Schematic of different scenarios of contact between the AFM tip and the outer surface of the 3D multilayer PEGDA hydrogel at the interface. The red arrows demonstrate the contact points changing as the AFM tip moves from one area to another on the surface.
Secondly, since the surface topography is considered rough for the AFM measurements, it could be argued that there may be variations in terms of how much material is being indented beneath the surface which may influence the estimate of E. Since the indentation depth into the sample surface is well over 1 µm, the effect of surface topography on the measured E may be considered negligible. S-12

Attempting to eliminate the outer surface topography
Attempting to eliminate the influence of the outer surface topography of the printed multilayer sample by reducing surface roughness, a cross-section of the multilayer sample was cut with 100 µm thick stainless steel blade in its hydrated state ( Figure S7a). However, the surface topography and roughness increased with peak to valley (trough) difference reaching maximum of 5.3 and 6.7 µm in both multilayer and monolithic samples, respectively ( Figure S7b), rendering it impossible to perform any sensible modulus measurements ( Figure S7c). Another way to improve surface topography was to freeze and cryo-section the sample as it may improve the surface flatness. 7 However, such approach will completely change the hydrogel properties and its real-life condition and will lead to E measurements that are effectively not representative of the real hydrated system. 8,9 S-13 S-14

Measuring light absorbance for prepolymer solution with fixed QY concentration
A prepolymer solution containing 9 mg/mL QY, 5 mg/mL LAP, and 200 mg/mL PEGDA was prepared. The solution was poured onto a coverslip at the bottom of a vat, and the UV light intensity was adjusted to 20.51 mW/cm 2 . After being exposed to UV light for 3 s, the excess solution was pipetted out from around the formed layer, and the surface of the layer was gently blotted with a paper towel. The thickness of the formed layer (known as cross-linking depth (Cd)) was measured and recorded as 46 ± 2 µm. The intensity of the UV light passing through the formed layer was measured and found to be 3.06 mW/cm 2 . Using eq S3 the light absorbance was calculated: 10 where (I0) is the initial intensity at the bottom of the vat and I(L) is the intensity of the light as it passes through the sample of thickness (L), the light absorbance was calculated to be 0.80.
A combination between eq S3 and Beer-lambert law eq S4 can be used to calculate the extinction coefficient (ε): 10 where (C) is the concentration of the photoabsorber in the prepolymer solution and L is the constant path light (which is equivalent to Cd in this case). The extinction coefficient was calculated to be 5,474.6 M -1 cm -1 . S-15

Measuring cross-linking depth as light passes through a cross-linked single layer hydrogel with 20 μm layer thickness
Initially, single layer hydrogel with a thickness of 20 µm and a QY photoabsorber concentration of 9 mg/mL in the prepolymer solution was printed using an intensity of 20 mW/cm 2 .The intensity of the UV light passing through the resulting cross-linked layer was measured to be 3.29 ± 0.1 mW/cm 2 (as depicted in Figure S9). Next, a glass slide was inserted into the printing vat and the UV light intensity was adjusted to 3.29 mW/cm 2 . The same prepolymer solution used to print the 20 µm layer thickness was poured onto the glass slide and exposed to the adjusted UV light for 3 s (as shown in Figure S9). The resulting cross-linked layer was measured using a micrometer and found to have a thickness of 3 ± 2 µm. This suggests that only the top 3 µm of each layer in the multilayer sample were double-exposed (as illustrated in Figure S9), which is consistent with our previous findings and supports the formation of a highly cross-linked planar sub-layer. 11 Based on our measurements, it appears that the intensity of the UV light passing through the prepolymer solution is insufficient to penetrate and double expose the entire previously formed

S-17
layer (as shown in Figure S10). The QY photoabsorber present in the prepolymer solution significantly absorbs light as it propagates through the solution during photo-cross-linking. This causes an exponential light intensity attenuation, 10,12 resulting in only a small portion of the top section of the previously formed layer near the interface being double exposed. Figure S10. Schematic representation of the projection lithography system depicting the layer formation process (a-c) and (d) area within the surface of the previously formed layer which is potentially double exposed to UV light just below the interlayer interface.

S-18
10. Optical microscopic images of outer surface topography of multilayer hydrogels Figure S11. Optical microscopic image of the printed multilayer hydrogel sample outer surface topography of a (a) cylindrical, (b) cuboid shape samples. S-19

Alternative options to LAP photoinitiator
There is a constant search for alternatives to metal-based photoinitiators in order to examine whether they can replace LAP and lead to a more homogeneous structure. If Irgacure 2959 was soluble in water instead of methanol, it could have been a viable option for tissue engineering studies. 1 > However, a more suitable photoinitiator with stable radicals such as TEMPO could be used to achieve well-controlled free-radical cross-linking. TEMPO's stable free-radical property is due to its steric bulk, which impedes the reaction of other free radicals, providing control over polymerization kinetics. 14 While TEMPO has been shown to be cytotoxic and unsuitable for tissue engineering, alternative options could be explored to examine if they may result in more homogenous structures.