High throughput computational evaluation of how scaffold architecture, material selection, and loading modality influence the cellular micromechanical environment in tissue engineering strategies.

BACKGROUND: In tissue engineering (TE) strategies, cell processes are regulated by mechanical stimuli. Although TE scaffolds have been developed to replicate tissue-level mechanical properties, it is intractable to experimentally measure and prescribe the cellular micromechanical environment (CME) generated within these constructs. Accordingly, this study aimed to fill this lack of understanding by modeling the CME in TE scaffolds using the finite element method.
METHODS: A repeating unit of composite fiber scaffold for annulus fibrosus (AF) repair with a fibrin hydrogel matrix was prescribed a series of loading, material, and architectural parameters. The distribution of CME in the scaffold was predicted and compared to proposed target mechanics based on anabolic responses of AF cells.
RESULTS: The multi-axial loading modality predicted the greatest percentage of cell volumes falling within the CME target envelope (%PTE) in the study (65 %PTE for 5.0% equibiaxial tensile strain with 50 kPa radial-direction compression; 7.6 %PTE without radial pressure). Additionally, the architectural scale had a moderate influence on the CME (maximum of 17 %PTE), with minimal change in the tissue-level properties of the scaffold. Scaffold materials and architectures had secondary influences on the predicted regeneration by modifying the tissue-level scaffold mechanics.
CONCLUSIONS: Scaffold loading modality was identified as the critical factor for TE the AF. Scaffold materials and architecture were also predicted to modulate the scaffold loading and, therefore, control the CME indirectly. This study facilitated an improved understanding of the relationship between tissue-level and cell-level mechanics to drive anabolic cell responses for tissue regeneration.

INTRODUCTION

Organ healing via regenerative medicine will afford revolutionary treatment for a myriad of diseases. The goal of tissue engineering (TE) is to drive native and introduced cells to produce a healthy, functional extracellular matrix (ECM) to repair and regenerate diseased native tissue. To enhance the regeneration potential, TE constructs are commonly laden with exogenous progenitor or stem cells. Consequently, the fate of these cells is paramount to establish long‐term biological function and mechanical integrity of the engineered tissue.

A major regulator of cell fate is mechanical loading. Localized stresses and strains have been shown to dictate cell viability, differentiation, and ECM deposition., , , As compared to other systems, the musculoskeletal system experiences a broad magnitude of mechanical loads. Consequently, cell fates in muscle, bone, articular cartilage, fibrocartilage, tendon, and ligament are all driven largely by mechanical cues. For example, models of bone fracture healing have used hydrostatic stress history and maximum principal strain history as mechanical measures to predict regeneration.,

This study focused on the specific example of regeneration of the annulus fibrosus (AF), a fibrocartilaginous component of the intervertebral disc (IVD). The biological and mechanical integrity of the AF is contingent on the production and maintenance of ECM by AF cells,, and diseased states of the IVD have been associated with a loss of tissue cellularity and dramatic changes to the organization and regeneration of the ECM., , Further, mechanical loading has been linked to inflammatory responses of AF cells, which may be critical for tissue homeostasis, or may invoke degenerative sequelae at supra‐physiological strains., , , The viability and ECM production of AF cells have been shown to depend on the magnitude and three‐dimensional combinations of mechanical loading. In vivo, the posterolateral AF experiences biaxial tensile strains of ~4% to 6% in the circumferential and axial directions and hydrostatic pressure generated by the adjacent nucleus pulposus.,

AF cells isolated from rabbits have demonstrated anabolic responses at maximum principal strains (ɛ 1) of 3% to 18%, and this response was maximized at 6% strain. At 1% strain, rabbit AF cells have demonstrated no significant changes in proteoglycan production, cell death, MMP‐1 expression, or MMP‐3 expression as compared to static loading. This remodeling window is supported in studies of human AF cells. Upregulation of catabolic factors associated with disc degeneration has been demonstrated at 20% strain. Decreased catabolic gene expression has been shown for human AF cells at 10% strain, and increased cell proliferation, collagen production, and glycosaminoglycan production has been reported at this strain magnitude. Accordingly, a maximum principal strain remodeling window of 3% to 18% was proposed for this study as a target for cell‐level loading to drive AF regeneration. Similarly, AF cells have exhibited anabolic responses for compressive hydrostatic strains (ɛ h < 0), , , , and an upper limit of 1 MPa compressive hydrostatic stress (σ h <  − 1 MPa) has been proposed for eliciting catabolic responses., Therefore, a hydrostatic stress remodeling window of 0 to 1 MPa was proposed for this study as a target for cell‐level loading to drive AF regeneration. This proposed CME target envelope based on strain and hydrostatic pressure is also in general agreement with previously reported micromechanical criteria for cartilage and fibrous tissue formation in fracture healings models.,

A ubiquitous strategy within the TE community is to fabricate composite constructs consisting of a biodegradable scaffold with a cell‐laden matrix. Such TE scaffolds have been engineered to replicate specific tissue‐level material properties of various musculoskeletal tissues,, , including AF., However, these scaffolds do not necessarily ensure that the mechanical loads induced at the cellular level are sufficient to drive cell survival, proliferation, and ECM formation. The relationship between tissue‐level loading and the cellular micromechanical environment (CME) is, therefore, essential to furthering our understanding of how best to design TE scaffolds. Yet, it is intractable to measure and prescribe the CME in cell‐laden matrices of TE scaffolds. The CME is three‐dimensional, heterogeneous, and dependent on scaffold loading, materials, and architecture; current experimental methods are not capable of accurately prescribing and/or measuring the CME. For example, in the aforementioned complementary experimental series, there is no physical method to know what CME is generated under global (ie, tissue level) loading of the TE scaffolds and whether or not that CME will be beneficial for regenerating the desired tissue. Optical strain measurement techniques, with image capture via high resolution digital camera or confocal microscopy, have been used with digital image correlation (DIC) to measure deformations in biological materials., , However, both DIC and confocal microscopy are not well suited for high‐throughput analyses. Digital photography techniques are limited in that they can only measure two‐dimensional surface strains and the resultant surface strains typically do not represent the complete deformation mapping within the scaffold. Confocal microscopy techniques may also be restricted by the opacity of the scaffold and hydrogel.

In addition to the experimental difficulties of measuring CME, TE experiments with dynamic mechanical loading may be prohibitively time consuming to characterize the relationship between tissue‐level and cellular‐level loading for a broad range of loading conditions. For example, in vitro investigation of TE scaffolds require sufficient culture time to elicit a measurable cell response, such as ECM production. Additionally, the need for complex apparatus to deliver precise multiaxial loading may limit a study group to successive cultures, which may take appreciable time to produce statistically powerful results. Subsequent study groups aimed to improve tissue regeneration will have a similarly protracted study duration. Therefore, in order to optimize the development timeline of tissue regeneration strategies, it is imperative that the most advantageous study groups are selected for experimental evaluation. However, there is currently no method to identify which particular scaffold design features and experimental conditions are most likely to drive improved tissue regeneration. For example, is the CME more sensitive to scaffold loading, materials, or architectural parameters?

In the absence of any feasible experimental methods, one possible tool to predict ECM formation in TE scaffolds is the finite element (FE) method. Cell fates in orthopedic tissues under mechanical loading have been modeled with FE in intervertebral disc, and bone fracture healing applications. The tissue‐level mechanics of TE scaffolds have also been studied using FE methods, , , and some of these models have been developed to predict mechano‐regulation of musculoskeletal regeneration., , However, there remains a need for a CME model that can: (a) be applied to all of the available volume that cells can occupy in heterogeneous TE scaffolds; (b) be applied parametrically to numerous candidate TE scaffold designs; (c) be broadly applied to a range of proposed target mechanics; and (d) be easily compared to in vitro cell cultures for validation. Therefore, in this study, a FE model of CME was developed to predict the regeneration potential of TE scaffolds. The influence of scaffold loading modalities, materials, and architectures on the theoretical healing potential were investigated. The results of this model were used to inform the design of a TE scaffold for AF regeneration in an ongoing experimental study.

MATERIALS AND METHODS

Scaffold base model

The development of a repeating unit cell model for evaluation of the CME has been previously published by our group. In brief, the unit cell is an idealized geometry of a three‐dimensional fiber deposition (3DF), angle‐ply scaffold which has previously demonstrated anisotropic material properties similar to the most relevant properties of native AF tissue (Figure 1A)., The unit cell model was parametrized based on scaffold architecture, materials, and loading. Base parameter values of all scaffold parameters are summarized in Table 1.

FIGURE 1

Scaffold model for the base geometry showing, A, the previously validated angle‐ply fiber scaffold; B, the double unit cell of the fiber scaffold with the definition of fiber angle (Θ); C, the final unit cell including the hydrogel infill showing the FE mesh and region of interest (ROI) for cellular micromechanical environment (CME) evaluation; and D, tri‐axial loading definitions of axial strain (ε a), circumferential strain (ε c), and radial pressure (σ r). The x‐, y‐, and z‐directions represent the axial, circumferential, and radial directions of the intervertebral disc (IVD), respectively

TABLE 1

Parameters and associated values for the scaffold base model

Category	Parameter	Symbol	Base model value	Ref.
Architecture	Fiber angle	Θ	34°	36
	Fiber spacing	S	1.0 mm	36
	Fiber diameter	D	0.3375 mm	36
	Layer height	–	0.6 × D	36
	Fiber contact radius	–	1.58 × D	36
Materials	Hydrogel elasticity	C 1	172 Pa	42
		C 2	383 Pa	42
	Hydrogel compressibility	D 1	3.41	42
		D 2	0.0806	42
	Fiber elastic modulus	E	265 MPa	36
	Fiber poisons ratio	ν	0.3	36
Loading	Axial strain	ɛ a	5.0%	‐
	Circumferential strain	ɛ c	5.0%	‐
	Radial pressure	σ r	0 MPa	‐

Note: Parameters are categorized as either architectural, material, or loading.

Fibers were prescribed linear elastic material properties of polycaprolactone (PCL) and the hydrogel infill was modeled with the compressible, second‐order reduced polynomial hyperelastic material properties of fibrin (Equation 1). where U is the strain energy potential, C 10 and C 20 are the fitted material elasticity constants, I1 is the first invariant of the strain tensor, D 1 and D 2 are the fitted material compressibility constants, and J is the Jacobian determinant of the deformation gradient tensor.

Tensile strains were applied to the model in the axial direction (x‐direction) and circumferential direction (y‐direction) to emulate the dominant in vivo loads experienced by the posterolateral AF., , , All nodes on the unit cell faces normal to the negative axial and circumferential directions were constrained against displacement in those respective directions; therefore, in‐plane sliding was allowable. Biaxial displacements were prescribed on the positive axial and circumferential faces to generate global strains. On the positive z‐direction face, all fiber nodes were constrained to equal z‐displacements and all fiber nodes on the negative z‐direction face were constrained against out‐of‐plane displacement. Previous results showed that constraining all nodes on the positive and negative z‐direction faces increased the region of interest (ROI) strain energy by 196% and 179%, respectively. To address this sensitivity of the ROI mechanics to the z‐direction boundary conditions, these two constraints were reviewed in the current study.

In the center of the unit cell, a ROI consisting of mesh elements was defined to contain all possible positions of cells within the hydrogel matrix with respect to the fiber architecture (Figure 1D). Using these ROI elements, a custom post‐processing script generated representative cell volumes of seeded progenitor cells in the hydrogel matrix (20 μm equivalent seed size) and evaluated the theoretical micromechanical environment of these cells.

CME evaluation

A CME post‐processing algorithm was developed to facilitate the evaluation of a constant three‐dimensional strain tensor for cell‐sized volumes in the ROI while maintaining the stability, accuracy, and efficiency of an FE model with larger and more complex elements. Specifically, the whole volume of the ROI was considered to characterize the CME in the scaffold. The deformation solution of the quadratic tetrahedral (C3D10H) elements in each model was reverse‐engineered to yield Green strain tensors for cell‐sized linear tetrahedral elements (C3D4). These linear tetrahedral elements had an effective seed size (20 μm) that is similar to the size of mature AF cells., From the cell‐volume strain tensor, the CME for all cell‐sized volumes in the ROI was categorized as either within (“satisfying”) or outside (“not satisfying”) the proposed target mechanics envelope. These target mechanics were derived from previously published 3D micromechanical criteria for anabolic responses of mature human AF cells. Specifically, the target mechanics envelope was based on maximum principal strain (3% < ɛ 1 < 18%), , , , and hydrostatic stress (−1 MPa < σ h < 0 MPa), , , , , as shown in Figure 2.

FIGURE 2

The proposed micromechanical target envelope for AF tissue regeneration based on hydrostatic stress and maximum principal strain criteria. The window of anabolic responses based on hydrostatic stress has a lower bound of −1 MPa (ie, hydrostatic compression) and an upper bound of 0 MPa., , , , , The window of anabolic responses based on maximum principal strain has a lower bound of 3% and an upper bound of 18%, , , ,

Parametric studies

Following validation, the loading, materials, and architecture of the base model were modified parametrically to investigate their relative influences on the predicted CME.

Scaffold loading

To evaluate the influence of biaxial loading on CME, the base model was prescribed the following biaxial strain conditions: (1) +6.0% and −6.0% axial strain with an array of circumferential strain from −6.0% to +6.0% in increments of 1.0% strain and (2) +6.0% and −6.0% circumferential strain with an array of axial strain from −6.0% to +6.0% in increments of 1.0% strain. These 48 combinations of biaxial strain were denoted as load array 1. The load ramps of these solutions were analyzed to yield a full series of biaxial loading conditions between −6.0% and +6.0% strain (eg, the solution for 6.0% equibiaxial strain contained the solution for 5.0% equibiaxial strain). Following preliminary results, two additional loading sets were considered in the study and used in the analyses of scaffold materials and architecture. Load array 2 was utilized to capture the most relevant range of load array 1 while minimizing the number of study points to reduce computational burden. Specifically, load array 2 was defined as the eight biaxial combinations of −5.0%, 0.0%, and +5.0% strain (excluding the unloaded condition). To capture the most pertinent scaffold loading modalities based on the results of load array 1, load array 3 was defined as the two conditions of equibiaxial (+5.0% axial and circumferential strain) and transverse‐constrained circumferential strain (+5.0% circumferential strain and 0.0% axial strain).

In addition to biaxial strain combinations, the influence of a compressive load in the radial direction was investigated by applying a pressure up to 1.0 MPa to the positive radial face following the full biaxial load. Due to numerical complexity, two cases of radial pressure were considered: prior to the biaxial strain and following the biaxial strain. This compressive load was considered for the base model under load array 3.

Scaffold materials

Eight conditions of modified material properties were evaluated under load array 2 (positive biaxial strains). Specifically, the following material property modifications were evaluated: (1, 2) upper and lower 95% confidence bounds of hydrogel compressibility coefficients (D 1 and D 2), (3, 4) upper and lower 95% confidence bounds of hydrogel elasticity coefficients (C 10 and C 20), (5, 6) 10‐fold increase and decrease of the hydrogel compressibility coefficients (C 10 and C 20), and (7, 8) increase and decrease of the fiber elastic modulus (E) by 20%.

Scaffold architecture

Modifications to four architectural parameters were explored. Scaffold fiber angles (Θ) ranging from 30° (increased biaxial asymmetry as compared to base model) to 45° (biaxial symmetry) were considered. The fiber spacing (S) in the scaffold was ranged from 0.6 to 1.4 mm, in increments of 0.1 mm and the fiber diameter (D) was studied from 0.20 to 0.45 mm (based on the range of common polymer fiber diameters produced via 3DF, , , ). Finally, the relative cell size was progressively increased to 10 times the original size. This generated effective architectural scale factors of 0.1 to 1.0 while maintaining the accuracy of the base model. Each scaffold architecture was prescribed the two loadings conditions of load array 3.

RESULTS

Base model

In the base model, the cell‐sized volumes within the ROI exhibited a distribution in hydrostatic strain and maximum principal strain (Figure 3). Specifically, the CME that satisfied the proposed target envelope was predicted for 7.2% of cell volumes (ie, a predicted target envelope, PTE, of 7.2% or 7.2 %PTE). A concentration of CME was observed for small positive hydrostatic strains (0%‐4%) and maximum principal strains ranging 5% to 15%. In equibiaxial 6.0% strain, the strain magnitudes experienced by some cell volumes exceeded 20%. The PTE varied between two subsets of cell volumes: (1) cells with direct contact with the PCL fibers and (2) cells with no direct contact with PCL fibers. These two subsets represented 88.4% and 11.6% of the total number of cells, respectively, and the PTE for these subsets were 6.1 %PTE and 15.65%PTE, respectively. At 5.0% equibiaxial loading, the base model had 8.1 %PTE. The two alternative boundary conditions with all nodes (hydrogel and fibers) constrained on the bottom and top faces resulted in increases in the PTE of 4.6 %PTE and 8.1 %PTE, respectively, as compared to the base model at 5.0% equibiaxial loading.

FIGURE 3

Unit cell solution for the base scaffold model with prescribed 6.0% equibiaxial strain: A, Mises strain contours of the full model; B, Mises strain contours of the region of interest (ROI); and C, distribution of cellular micromechanical environment (CME) as a function of hydrostatic strain and maximum principal strain for all ROI cell volumes. The red dashed box indicates the proposed micromechanical target envelope which contains 7.2% of the CME distribution for this model (7.2 %PTE; predicted target envelope)

Scaffold loading

The PTE was found to vary as a function of equibiaxial strain magnitude (Figure 4). Despite 7.2 %PTE for 6.0% equibiaxial strain, a peak of 9.5 %PTE was observed at 2.7% equibiaxial strain. The PTE rapidly and monotonically increased to this peak from 0.0 %PTE at ~0.7% equibiaxial strain and appeared to decrease approximately linearly following the peak. As compared to all cell volumes, the PTE for cells with no attachment to scaffold fibers (88.4% of cell volumes) was lower for all load fractions with a peak value of 8.9 %PTE. Cells with fiber contact (11.6% of cell volumes) exhibited a peak value of 15.6 %PTE at 5.0% equibiaxial strain.

FIGURE 4

Predicted target envelope (PTE) as a function of: A, loading magnitude in equibiaxial tension and B, biaxial strain ratios for the base scaffold model. In A, the peak and final value of PTE are identified with red circles for all cells volumes, cell volumes with no contact to fibers (unatt. cell volumes), and cell volumes in contact with fibers. In B, contours are shown in increments of 1.0 %PTE. Two local maxima are indicated with a white “X” and labeled with the corresponding %PTE and biaxial strain in parentheses

The base model prescribed with a biaxial loading array (load array 3) demonstrated a clear relationship between the biaxial load and the PTE (Figure 4B). The regions experiencing relatively low loading (less than ~1% strain in either direction) had no PTE (0.0 %PTE). Two local peaks in the PTE were observed, one in biaxial tension (9.8% for 2.8% axial and 4.2% circumferential strains) and one in the biaxial compression (10.5% for −3.6% axial and −5.4% circumferential strains). The biaxial compression peak was higher in magnitude and broader in loading range as compared to the biaxial tension peak. The average angle of the maximum principal strain directions from the loading plane were distinctly different between the two peaks (0.3° and 78° for biaxial tension and compression, respectively).

Numerical instabilities occurred when a radial pressure was applied to the model prior to biaxial strain. However, radial pressures of at least 50 kPa were resolvable following the application of biaxial strains of different magnitudes (Figure 5). In the absence of biaxial strain, the radial pressure of 37.5 kPa yielded 17.2 %PTE which was greater than any equibiaxial strain in the absence of radial pressure (maximum of 9.5 %PTE). At all magnitudes of radial pressure, a peak in the PTE was observed; for all non‐zero pressures, this peak was at 2% equibiaxial strain (maximum of 66.5 %PTE for 50 kPa radial pressure). Similarly, the influence of radial pressure was most pronounced near the peak in the PTE, as evidenced by a 40.5 %PTE increase from 0.0 kPa radial pressure to 12.5 kPa.

FIGURE 5

Predicted target envelope (PTE) for combined loading conditions of radial pressure and equibiaxial strain. The plotted lines shown the PTE for constant pressures in increments of 12.5 kPa and all non‐zero pressure lines are sampled increments of 1% equibiaxial strain. The %PTE is shown as a function of the load magnitude for constant radial pressures

Scaffold materials

Overall, the changes in PTE within the 95% confidence bounds in hydrogel elasticity and compressibility were less than 0.3 %PTE (Figure 6). Nonetheless, distinct trends were observed as a function of scaffold loading. The lower confidence bound of hydrogel elasticity and compressibility both demonstrated increased PTE in biaxial tension and decreased PTE in biaxial compression, with a gradual change in between. The inverse trend was observed for the corresponding upper confidence bounds. These trends are supported by the more extreme changes to the hydrogel properties.

FIGURE 6

Changes to the predicted target envelope (PTE) due to perturbations of material mechanical properties. The considered material perturbations were to the hydrogel elasticity constants (C) and hydrogel compressibility constants (D). The subscripts −95 and +95 represent the lower 95% and upper 95% confidence bounds, respectively. The subscripts −20 and +20 represent a 20% decrease and 20% increase of the base material values, respectively. C 0.1 and C 10 represent the tenfold increase and decrease in hydrogel elasticity, respectively. For each material condition, trends are shown as a function of biaxial loading

Similar to the upper 95% confidence interval, the PTE was increased in biaxial compression and decreased in biaxial tension for a tenfold increase in C 1 and C 2, and vice versa for a decrease in C 1 and C 2 (Figure 6). Decreasing the order of magnitude of the hydrogel elastic coefficients (C 1 and C 2) tenfold resulted in changes to the PTE within ±0.7 %PTE for all loading conditions. Conversely, a tenfold increase in C 1 and C 2 changed the PTE by ±0.9 to ±9.0 %PTE.

All conditions of changes in fiber elastic modulus changed the PTE by less than 0.02%. However, the changes in fiber elasticity were associated with proportional changes in stress on the unit cell For example, the 20% increase in fiber elastic modulus resulted in a 20% increase in the circumferential and axial stresses required to generate the same global strains.

Scaffold architecture

The PTE as a function of the four architectural parameters for 5% equibiaxial strain and 5% uniaxial tension (while constrained in the circumferential direction) are shown in Figure 7. Overall, the fiber spacing, fiber angle, fiber diameter, and architecture scale demonstrated a total range of 3.8 to 17.2 %PTE. For these load conditions, the attached fibers exhibited a greater PTE than the unattached fibers at full load in all considered architectures. However, for all variants of fiber angle and fiber diameter, the attached fibers only accounted for 9.4% to 12.7% of the total cell volumes in the ROI. The fraction of attached cells increased for decreasing fiber spacing (from 7.1% at 0.6 mm to 23.0% at 1.4 mm) and decreasing fiber scale (from 32.9% at 0.2 scale to 11.6% at 1.0 scale). These increases in attached fibers corresponded with an observable divergence of the PTE for all cells and unattached cells.

FIGURE 7

Predicted target envelope (PTE) in the scaffold model for four architectural parameters: fiber angle, fiber spacing, fiber diameter, and architecture scale. In the left column of plots, the %PTE is shown as a function of the architectural parameters. A selected value of each parameter is indicated with black arrows at the top and bottom of the plot. For each of these selected values, the %PTE is shown as a function of the strain magnitude in the right column of plots. In each plot, the PTE is shown for equibiaxial 5.0% strain (EQB, dark blue) and axially‐constrained circumferential 5.0% strain (TCC, light blue) as well as for: all cells (solid lines), cells not attached to fibers (dashed lines, denoted as unatt.), and cells attached to fibers (dotted lines, denoted as att.)

For the selected values of fiber angle, fiber spacing, and architecture scale, the axially‐constrained circumferential strain exhibited a greater PTE at full load as compared to equibiaxial strain. However, the loading profile data show that the equibiaxial strain has greater peaks in PTE during the load ramp. Many other trends and features were predicted as a function of the four architectural parameters, as detailed below.

Fiber angle

In equibiaxial tension, PTE appears to show no strong correlation with fiber angle, regardless of cell attachment (range of 5.4‐8.1 %PTE for all cell volumes). In axially‐constrained circumferential strain, the PTE generally decreased with increasing fiber angle for all conditions of cell attachment (range of 4.2‐9.2 %PTE for all cell volumes).

The loading profile data for a 20° fiber angle (Figure 7) showed a similar trend to the base model in equibiaxial strain (Figure 4), including a rapid increase in PTE (starting at a load fraction of ~0.1). Likewise, in axially‐constrained circumferential strain, a rapid increase in %PTE was shown, however, beginning at a greater load fraction (~0.3). The axially constrained circumferential strain data demonstrated a monotonic increase in PTE with increasing load fraction.

Fiber spacing

In equibiaxial strain, the PTE generation exhibited a peak as a function of fiber spacing (8.8 %PTE at 0.9 mm fiber spacing). This maximum had a pronounced reduction in the model's PTE in the direction of reduced fiber spacing (decreasing to 5.1 %PTE at 0.6 mm fiber spacing). In axially‐constrained circumferential strain, a monotonic decrease in PTE was observed with increasing fiber spacing (range of 3.8‐13.0 %PTE for all cell volumes). In both loading cases, similar trends were observed for attached and unattached cells.

The equibiaxial loading profile data for 0.6 mm fiber spacing demonstrated a peak in PTE of 15.3 %PTE at 0.30 load fraction. Near this peak, the unattached cells exhibited a higher PTE than the attached cells. A similar phenomenon was shown in axially‐constrained circumferential strain, which also appeared to plateau at ~13 %PTE (all cell volumes) for load fractions greater than 0.6.

Fiber diameter

In both loading conditions, variations in fiber diameter showed no distinct change in %PTE, regardless of cell attachment (range of 4.2‐9.2 %PTE for all cell volumes). The loading profile data for the 0.2 mm fiber diameter demonstrated a monotonic increase in PTE when equibiaxial strain was applied, and a plateau in PTE in uniaxial loading (with circumferential strain constraint). Both loading conditions provided prediction of less than 10 %PTE for all load fractions.

Architecture scale

The ROI strain energy for architecture scales of 0.2 and greater were within 1.0% of the base model in both loading conditions. However, the 0.1 scale ROI strain energy differed by at least 257% from the base model and was omitted. The circumferential and axial reaction forces varied proportionally with the respective constrained areas such that the applied stresses varied by less than 0.02% for all architecture scales.

In both equibiaxial strain and axially‐constrained circumferential strain, the PTE generally increased with reduced scale factor for all cells as well as the cells both attached and unattached to fibers. Loading profile data for the 0.2 scale demonstrated peaks of 17.2 and 15.7 %PTE for the equibiaxial strain and axially‐constrained circumferential strain conditions, respectively. In the axially‐constrained circumferential strain load ramp, the unattached cell volumes had a greater PTE than attached cell volumes up to a load ratio of 0.86.

DISCUSSION

In this study, a FE model was implemented to predict the CME within a TE scaffold. Specifically, a repeating unit cell of an angle‐ply laminate scaffold for AF regeneration was prescribed a variety of loads, materials, and architectures to assess the relative influences of these factors on the CME. The model with base parameters (optimized to match the tissue‐level properties of native AF) exhibited a distribution in maximum principal strain and hydrostatic strain (Figure 3C). Of this distribution, a fraction of the hydrogel volume fell within the proposed target envelope CME for AF regeneration. This level of PTE changed for all considered loads, materials, and architectures, however, the relative sensitivity of the PTE varied between these design factors.

The base model demonstrated a distribution of the complex, three‐dimensional mechanical state within the hydrogel matrix of the TE construct (Figure 3C). This heterogeneity can be attributed to the composite architecture. The fibrous scaffold is necessary to provide tissue‐level mechanical integrity to the TE construct. However, the fibers also transmit mechanical stimuli to progenitor cells seeded in the matrix and induce heterogeneity in the cell‐level mechanical environment. The results of this study predicted that the scaffold dramatically influences the CME; in the base model, some cell‐level strains exceeded 20% for just 6% global equibiaxial strain. Further, the transmission of mechanical loads to progenitor cells is highlighted by the consistently higher rate of PTE for cells attached to fibers as compared to cells with no fiber attachment (Figure 4A). This prediction of an enhanced cell response is consistent with the results of experimental cell cultures on fibrous scaffolds., ,

What is a sufficient level of PTE for AF repair?

Mature, healthy AF has a cell density of around 9000 cells/mm3., Based on the tetrahedral volumes in this study, only 1 %PTE is theoretically required to maintain this healthy cell volume. However, it is unlikely that such a low fraction of satisfactory CME (ie, CME that meets the requisite mechanics for ECM formation) would be sufficient for AF repair. First, a higher cell density may be required for the enhanced ECM production required for AF healing, as compared to the maintenance of healthy tissue. Second, seeded progenitor cells are dispersed within the hydrogel and do not occupy all available volumes. Therefore, the intersection between the distribution of satisfactory CME and the distribution of progenitor cells within the scaffold would likely lead to a lower fraction of cells with a satisfactory CME than predicted. For example, in the current study the satisfactory CME appeared to be concentrated around the scaffold fibers, however, the seeded cells may not be equally concentrated around these fibers. Third, a 1 %PTE indicates that 99% of cell volumes experience a CME outside of the target mechanics. The cells that occupy these regions may not contribute to the AF regeneration and could potentially produce some deleterious outcomes, including: apoptosis, cellular inactivity, altered cellular phenotypes, catabolic responses, and inflammatory cytokine release. Moreover, even if a small fraction of cells are apoptotic, this may induce apoptosis throughout the scaffold, regardless of CME.

The proposed CME target region for AF regeneration is also likely to influence the level of PTE. This target region was based on CME criteria from published data of uniaxial strain and hydrostatic stress experiments. These data resulted in discrete boundaries of the proposed AF target region. However, this is likely a simplification of the underlying continuous response of AF cells to three‐dimensional mechanical stimuli. Further, different cell types or variations in the cytoskeletal morphology may lead to deviations in the requisite CME to elicit and anabolic response. However, possible disparities in the specific anabolic range for cells with different cytoskeletal structures falls outside of the length scale (ie, subcellular modeling) consideration for this study. The CME criteria also did not include data from combined loading, which may alter the anabolic response of AF cells. However, as the influence of the micromechanical environment on AF remodeling is further explicated the literature, the CME target presented in this model may require revision. For instance, this study serves as a design tool for an ongoing experimental series of TE scaffolds seeded with AF cells and cultured with a prescribed mechanical loading regime. Results from experiments such as these may be leveraged to recursively evaluate the validity of the proposed remodeling windows.

Based on these factors, there is likely no distinct PTE threshold for AF repair. Accordingly, maximizing the level of PTE appears to be the most suitable approach to optimize the level of satisfactory CME. The highest PTE in this study resulted in only 33.5% of the hydrogel volume with an insufficient CME, which seems an acceptably low fraction of non‐compliant CME. If the proposed CME target correlates well with anabolic cell responses, the highest levels of PTE are predicted to maximize the probability of a successful regenerative response within the TE scaffold.

The PTE throughout the equibiaxial load profile demonstrated that the CME is highly dependent on the magnitude of scaffold loading (Figure 4A). The peak in PTE suggests that an optimal load magnitude exists. At lower loads, the PTE was dramatically reduced and converged toward zero. At higher loads, the PTE decreased more gradually. Given the stochastic distribution of in vivo loading, these results suggest that over‐prescribing mechanical stimuli may best ensure that the target mechanics are met. Further, the physiological loads in vivo are cyclic and, therefore, many CME conditions in the load profile would be experienced during each cycle.

The biaxial loading array also revealed an influence of the relative magnitudes of in‐plane loads on PTE (Figure 4B). Two distinct regions of biaxial loads presented elevated rates of PTE with similar peak values (one in biaxial tension and one in biaxial compression). The major difference between these two regions was the average orientations of the maximum principal strains. The proposed CME target envelope in this study aimed to promote anabolic responses in AF cells for enhanced ECM production, including collagen, which is a major component of AF. Maximum principal strains were used as a mechanical criteria for predicting the anabolic responses. However, the direction of this maximum principal strain has substantial implications for repair because it has been hypothesized to dictate the orientation of collagen fibers in organs, and TE scaffolds. Because collagen in healthy AF is highly oriented in the circumferential‐axial plane, these results suggest that the peak in PTE associated with biaxial tension in the circumferential‐axial plane (similar to in vivo loading of the AF) is more likely to restore the native collagen structure of the AF. Accordingly, subsequent studies on scaffold loading and architecture were focused on this biaxial tensile region to minimize computational time. Specifically, load array 3 was used to encapsulate the region of peak biaxial tensile strain.

The peak in tensile biaxial strain (2.8% axial strain and 4.2% circumferential strain) may also reflect the dominant constraints of AF tissue in vivo; AF is constrained by the relatively rigid vertebral bodies in the axial direction and by the relatively compliant adjacent AF tissue in the circumferential direction. Accordingly, it is possible that the biomimetic architecture is replicating the structural organization of the native AF tissue that regulates the anabolic responses of AF cells. This peak also occurred within the typical range of physiological strains experienced by the AF (up to 6%), further suggesting that the scaffold may be reproducing the healthy mechanoregulatory environment of the AF.

Although the dominant physiological loads experienced by the AF are in the circumferential‐axial plane, the AF also experiences a radial pressure from the nucleus pulposus. The nucleus pressure dissipates radially from the inner to outer annulus, such that the pressure typically ranges from of 0 to 1 MPa. Accordingly, this additional, third dimension of loading was considered in the study. Validation of the unit cell model indicated that the ROI CME may be sensitive to the radial direction boundary conditions. Accordingly, the specific values of applied radial pressure in this model may not translate to the true scaffold boundary conditions in vitro or in vivo. Regardless, the underlying trends give insight to the influence of radial pressure on the CME. It was evident that radial pressure alone could produce greater PTE than the peak in pure biaxial strains. This result can be attributed to the three‐dimensional mechanical equivalence of uniaxial compression and biaxial tension. In the absence of radial pressure, the CME distribution for biaxial tension showed a dense population of cell volumes with hydrostatic strains slightly more positive than the CME criteria (Figure 3). It is intuitive that superimposing a compressive load (ie, generation of more negative hydrostatic strains) would result in a dramatic increase in PTE; the combination of biaxial tension and radial pressure provided the greatest PTE (66.5 %PTE). Moreover, the influence of radial pressure appeared to be exacerbated for biaxial strain ratios with greater PTE.

The confidence intervals of the hydrogel mechanical properties yielded minimal changes to the PTE (Figure 6). These confidence intervals represent reasonably expected variation in the hydrogel properties due to fabrication or composition. This variation is unlikely to lead to substantial changes in the CME and resultant satisfactory mechanics. Conversely, the tenfold changes in magnitude of the hydrogel elasticity may represent the selection of a different hydrogel material. These material alterations did result in appreciable changes to the PTE. For example, increasing the hydrogel stiffness increased the PTE by up to 9.0 %PTE (from the baseline of 8.1 %PTE). This increase in PTE can be attributed to the stiffer hydrogel sharing a greater proportion of the global loading with the scaffold, however, with little impact on the global scaffold mechanics. As such, the results suggest that changing the hydrogel material may be an advantageous method to modify the CME without modifying the overall mechanics of the composite construct. Additionally, the studied variations in the matrix material properties may reflect mechanical changes associated with the temporal adaptation of the ECM within the scaffold. For example, collagen deposition would be expected to stiffen the matrix leading to local changes in the CME and, therefore, may alter the optimal scaffold loading parameters as the tissue matures.

Changes to the fiber elasticity had negligible effect on the PTE for a consistent scaffold strain. This intuitive result may be advantageous because TE processes can alter material properties, such as thermal degradation during additive manufacturing. However, the changes to fiber elasticity also altered the overall scaffold mechanics, as observed by proportional changes to the necessary stresses (ie, reaction forces) for a given deformation. Though changes to material properties lead to negligible changes to the CME for a given strain, in situ the scaffold deformation is driven by loading. Therefore, the scaffold materials are more likely to have an indirect influence the CME by changing the scaffold deformation for a given mechanical load.

The fiber angle, fiber spacing, and fiber diameter all had minimal direct influence on the CME. Although the PTE demonstrated some correlation with these parameters, the trends were sensitive to the ratio and magnitude of biaxial loading. The nature of these trends may be associated with the competing factors of: (a) fiber surface area per unit volume (ie, the fraction of cells in proximity to the scaffold fibers), (b) the total volume fraction of matrix in the construct, and (c) the deformation of each fiber segment. For example, for a given construct loading, increasing the fiber spacing resulted in increased fiber strain but a corresponding decrease in surface area per unit volume. Overall, changing the fiber angle, spacing, or diameter are not good candidates for enhancing the PTE generation.

In, contrast, reduction of the overall architecture scale (ie, maintaining a constant ratio of fiber spacing and fiber diameter) consistently increased the PTE. This increase can be attributed to an increase in fiber surface area per unit volume, while maintaining the fiber deformation. Accordingly, reducing the fiber scale appears to be a viable methods to tailor the CME in TE constructs. The scale factor with the greatest PTE (0.2 scale) corresponded to a fiber diameter of 67.5 μm. Extrapolation of the scale results suggests that even smaller fiber diameters may yield smaller architecture scales and further improved level of satisfactory CME. In order to create these fine architectures, technologies such as melt electrowriting can be used to precisely deposit fibers of PCL and other biodegradable polymers with diameters of 1 to 100 μm.

Which design factor is most critical for control of the CME?

The critical determining factor for controlling the CME in the TE scaffold was the multiaxial loading modality of the scaffold, as demonstrated by the combined biaxial strain and radial pressure. Reducing the overall architecture scale could also be leveraged to enhance the CME while maintaining the tissue‐level mechanics of the TE construct. The remaining material and architectural alterations (hydrogel elasticity, fiber elasticity, fiber angle, fiber spacing, and fiber diameter) may be utilized to tailor the global scaffold properties with minimal influence on the CME. The tissue level properties are essential for mechanical support to maintain the overall stiffness and range of motion of the spine. However, these tissue‐level loads will dictate the scaffold deformation and indirectly influence the CME. The proposed pathway of influence for scaffold materials, architecture, and loading modality is summarized in Figure 8.

FIGURE 8

Diagram of the proposed influence of tissue engineering (TE) scaffold design factors on the cellular micromechanical environment (CME). The scaffold materials and architecture strongly influence the scaffold loading modality, which strongly influences the CME. The scaffold materials and architecture only had a weak influence on the CME

Implications for the design of experiments and implants

The strong influence of scaffold loading on CME indicated the importance of delivering physiological loads to TE constructs to stimulate regeneration. Accordingly, in the context of in vivo implants, the attachment conditions of TE construct may be essential for the generation and maintenance of new tissue. Suturing is a simple and ubiquitous surgical technique which may afford mechanical integration of TE implants. However, sutures only provide discrete attachment points between the implant and the adjacent tissues that may not transmit loads optimally. A promising range of products for complete mechanical integration of an implant surface to an adjacent tissue are bioadhesives,, which may afford more consistent implant attachment as a continuum.

Similarly, in the context of in vitro tissue cultures, there is a need for advanced experimental apparatus to apply and measure the requisite mechanical loads. The resources necessary for such equipment may prohibit experimental groups from being cultured concurrently. This limitation exemplifies the need for predictive tools such as this to understand TE results and drive informed design changes. An example of this model‐aided design was demonstrated in the scaffold loading results; by visualizing the distribution in predicted CME, the radial pressure was identified as a potential method to increase the level of PTE, and subsequently verified. A similar approach may be used to drive experimental design. In the complementary series of TE scaffold cultures, there was no available tool to measure or predict the CME generated in the scaffold under global biaxial loading. The results of this study provided a rationale to select specific values of biaxial strain that would theoretically maximize the level of PTE in the scaffold.

The relationship of materials and architecture with the CME also indicates that the tissue‐level and cell‐level mechanics can be tailored relatively independently. Therefore, it seems promising that design constraints for both scaffold mechanics and CME can be met simultaneously. This result also suggests that variable architecture design could be implemented to afford greater control of the scaffold loading. For example, the scaffold architecture could be varied near stress concentrations in an implant geometry or varied in order to drive a gradient in tissue phenotype.

Study assumptions and limitations

The lack of experimental validation remains a limitation of this work. The need for complex experimental apparatus and statistically powerful study groups inhibits a thorough validation of all presented results. Nonetheless, pertinent scaffolds designs (ie, loading regimes, materials, and architectures) will be experimentally cultured to validate key results. Histological analyses can be conducted on experimental cultures to quantify the ECM formation as an analog of the model PTE. Further, histological images can be analyzed spatially to assess the influence of cell‐fiber attachment. As previously discussed, the extensive time and resources required for cell cultures highlights the importance of predictive models to understand the relationship between tissue‐level and cell‐level mechanics. Ultimately, this model is an advantageous tool to explain and interpret TE results, and inspire hypotheses for improved TE strategies.

Several assumptions were made in the model that idealize the TE scaffold for computational practicality. However, it is possible these assumptions result in model limitations that fail to capture the complex and variable behavior of physical scaffolds. The fiber scaffold model idealized fibers from fused deposition as perfectly cylindrical and, therefore, does not account for printing flaws, such as fiber sagging and fiber topology. Similarly, both the fibers and hydrogel were assumed as heterogeneous and isotropic materials, however, TE fabrication methods may induce some level of heterogeneity and anisotropy. The hydrogel was assumed to completely fill the scaffold with no voids or flaws and was assumed to fully bond with the PCL fibers at all material interfaces. It follows that progenitor cells seeded in the hydrogel are also assumed to perfectly attach to the fibers. The validity of these simplifications are dependent on the quality of the fabrication method, however, the model may not capture alterations to the CME due to imperfections in a fabricated scaffold. The material properties of the matrix were simplified as an isotropic, continuum solid, and the measures for the CME criteria (ie, hydrostatic strain and maximum principal strain) were a result of this simplification. Future work may also incorporate other mechanoregulatory factors, such as osmolarity or oxygen tension, into the model and CME evaluation. The presented work also predicted the ECM for the acute phase of the regenerative response. However, it is expected that the matrix material properties would temporally evolve due to tissue adaptation (ie, ECM synthesis), which may be of interest to further enhance the model.

CONCLUSION

This study demonstrated high‐throughput, computational analyses to predict the relationship between the tissue‐level and cell‐level mechanics of TE scaffolds with prescribed loading, materials, and architectures. The scaffold loading modality was identified as the most pertinent factor in TE of the AF. Scaffold materials and architecture were also predicted to control the scaffold loading, and therefore the CME indirectly. Overall, this study facilitated an improved understanding of the relationship between tissue‐level and cell‐level mechanics, which may be utilized to tailor the CME, drive anabolic cell responses, and promote tissue regeneration. The theoretical framework presented in this study is highly tailorable and can be adapted to alternative TE strategies or incorporated in larger scale biomechanical models. Ultimately, this tool provides a CME‐based rationale to predict which TE study conditions are most likely to leverage improved tissue regeneration.

CONFLICT OF INTEREST

The authors declare no conflicts of interest to declare.

AUTHOR CONTRIBUTIONS

Mitchell I. Page, Peter E. Linde, Christian M. Puttlitz: Fully involved in the study and preparation of the manuscript.