Development of a Finite Element Model to Predict the Cellular Micromechanical Environment in Tissue Engineering Scaffolds

Cell fate in tissue engineering (TE) strategies is paramount to regenerate healthy, functional organs. The mechanical loads experienced by cells play an important role in cell fate. However, in TE scaffolds with a cell-laden hydrogel matrix, it is prohibitively complex to prescribe and measure this cellular micromechanical environment (CME). Accordingly, this study aimed to develop a finite element (FE) model of a TE scaffold unit cell that can be subsequently implemented to predict the CME and cell fates under prescribed loading. The compressible hyperelastic mechanics of a fibrin hydrogel were characterized by fitting unconfined compression and confined compression experimental data. This material model was implemented in a unit cell FE model of a TE scaffold. The FE mesh and boundary conditions were evaluated with respect to the mechanical response of a region of interest (ROI). A compressible second-order reduced polynomial hyperelastic model gave the best fit to the experimental data (C10 = 1.72×10-4, C20 = 3.83×10-4, D1 = 3.41, D2 = 8.06×10-2). A mesh with seed sizes of 40 μm and 60 μm in the ROI and non-ROI regions, respectively, yielded a converged model in 54 minutes. The in-plane boundary conditions demonstrated minimal influence on ROI mechanics for a 2-by-2 unit cell. However, the out-of-plane boundary conditions did exhibit an appreciable influence on ROI mechanics for a two bilayer unit cell. Overall, the developed unit cell model facilitates the modeling of the mechanical state of a cell-laden hydrogel within a TE scaffold under prescribed loading. This model will be utilized to characterize the CME in future studies, and 3D micromechanical criteria may be applied to predict cell fate in these scaffolds.


Material testing of fibrin
Hydrogel specimens were divided into two groups for unconfined compression (n = 9) and confined 13 compression (n = 7) testing. Prior to testing, the diameter and height of each specimen was measured with 14 digital calipers (CD-6"PSX, Mitutoyo, Kawasaki, Japan  22 In ABAQUS (Dassault Systèmes, Vélizy-Villacoublay, France), stress-strain data for confined 23 compression tests (data up to 5.0% strain) were paired with unconfined compression tests (data up to 33% strain) and imported as uniaxial compression and volumetric compression data, respectively. The data were 1 then fit to polynomial (first order to third order), Ogden (first order to third order), reduced polynomial 2 (first order to third order), and Van der Waals compressible hyperelastic models. The most appropriate 3 model for finite element implementation was selected based on model stability and fit to experimental data.

4
The averaged material models (elasticity and compressibility independently) were compared to the 5 experimental data resampled 0.1% strain increments with lack-of-fit F-tests (α = 0.05 for all tests).  20 Validation of the unit cell size was conducted prior to mesh refinement because preliminary data 21 indicated that mesh convergence was dependent on the unit cell size. All unit cell validation was performed 22 with a uniform quadratic tetrahedral (C3D10) mesh with 40 µm seed size, which was subsequently 23 validated. The selected unit cell size was then prescribed a series of mesh sizes to evaluate convergence of 24 the mesh and ROI mechanics. Uniform meshes of 40 µm, 35 µm, and 30 µm seed size were considered.

Mesh refinement
1 maintained for CME evaluation. The remaining unit cell was seeded with increasing size from 40 µm to 2 250 µm. For all meshes, the element growth rate within the ROI was set to 1.0 and the element growth rate 3 outside of the ROI was set to 5%. The only hyperelastic models with extensive elastic and volumetric stability were the first-order 16 reduced polynomial (neo-Hookean; all data fits were stable) and second-order reduced polynomial (the 17 elastic data from one sample was unstable and omitted, resulting in n = 6 for this model). Of these two 18 models, both exhibited no statistically significant evidence for lack of fit to the volumetric data (p = 0.99 19 for both). Similarly, for the elastic data, there was no statistically significant evidence for lack of fit of the 20 second-order model to the experimental data (p = 0.99). However, the first-order model had statistically 21 significant evidence of lack of fit to the experimental data (p < 0.001). Based on the stability and accuracy 22 measures, the second-order reduced polynomial model was selected for subsequent analyses ((1; Figure 3).
= 10 ( 1 -3) + 20 ( 1 -3) 2 + 1 ( -1) 2 + 2 ( -1) 4 (1) 1 where U is the strain energy potential, C 10 and C 20 are the fitted material elasticity constants, I 1 is the first 2 invariant of the strain tensor, D 1 and D 2 are the fitted material compressibility constants, and J is the 3 Jacobian determinant of the deformation gradient tensor. Average fitted coefficients for the second-order 4 reduced polynomial model were obtained by averaging the coefficients from all samples ( Table 1).

14
The results of unit cell validation demonstrated that cell sizes greater than 1.5-by-1.5 had ROI strain 15 energy within 0.01% of the 3 by 3 unit cell (Figure 4). The ROI strain energies of the 1-by-1 unit cells with 16 centered and offset ROI were 38% and 65% greater, respectively, than the 3-by-3 unit cell. Due to the observed influence of geometric boundary conditions on model stability in preliminary studies, the 2-by-2 1 unit cell was selected as the most suitable size for this study. All three of the reduced boundary conditions 2 on the in-plane (x-and y-direction) faces resulted in changes to the ROI strain energy of 0.5% or less.
3 4 Figure 4. Validation of unit cell size, mesh refinement, and unit cell thickness for the scaffold 5 CME model. For each validation, the convergence of strain energy was weighed against CPU 6 time. In the unit cell validation, a centered ROI was shown to be similar to an ROI offset by 1/8 7 of a unit cell in order to improve boundary conditions. 8 9 The mesh refinement yielded total strain energies within 1% of the 30 µm mesh for all considered 10 meshes. Similarly, the ROI strain energy was within 1% of the 30 µm mesh for all meshes and exhibited 11 apparent convergence for meshes finer than 100 µm. Therefore, the uniform 40 µm mesh was deemed 12 acceptable as compared to finer uniform meshes and the 40 µm ROI mesh was also accepted for all coarser 13 meshes. To select the most appropriate mesh with a 40 µm ROI size, the ROI strain energy and 14 computational time were considered. The 60 µm mesh (0.9 hours CPU time) had ROI strain energy within 15 0.5% of the uniform 30 µm mesh (10.6 hours CPU time). Therefore, the mesh with 60 µm general size and 16 40 µm ROI size was selected for all subsequent analyses.

17
The ROI SE demonstrated a monotonic increase as a function of bilayer count; the one-, two-, and 18 three-bilayer scaffolds yielded SE of 68%, 87%, and 93%, respectively, of the four bilayer scaffold SE.

19
Based on these results, the two bilayer (N = 2) was selected for all subsequent model validation. The fully 20 constrained boundary condition on both out-of-plane (z-direction) faces failed to converge. When 21 considered independently, the fully constrained boundary condition on the top and bottom faces resulted in 22 increases to the ROI strain energy of 179% and 196%, respectively. All of the modified boundary conditions 23 in the out-of-plane direction (z-direction) yielded changes in the total model strain energy of 0.03% or less.