In Silico Tissue Engineering: A Coupled Agent-Based Finite Element Approach

Introduction

Tissue engineering and regenerative medicine, despite all challenges, has been one of the fastest growing fields over the past 20 years. The dominant approach in the field has been to develop functional tissue-engineered constructs by mimicking the form of the native tissue, ¹ whereby research has been strongly focused on developing new biomaterials, ² improving manufacturing techniques,^3–5 and improving the techniques in cellular biology. ⁶ However, hidden in the “form follows function” dictum, are biological contributions to the function, which are the caveats to reverse engineering the tissue construct, since these missing parts of the puzzle expand the solution space exponentially. Fortunately, in silico tissue engineering merits potential for improving our understanding of the tissue construct by narrowing down the solution space and by optimizing the manufacturing processes,^1,7 specifically in the case of TEVGs, where the net of mechanobiological events, involved in regulating the growth and remodeling of the tissue, stretches across several temporal and spatial length scales and conceptualizing these processes is a daunting task.

TEVGs are emerging as an alternative to saphenous vein grafting and autologous blood vessels.^8,9 However, the application of TEVGs, specifically in the case of small-caliber arteries (<6 mm), remains limited due to low patency rates^9,10 associated with poor mechanical strength and compliance mismatch.^11–13 Achieving the mechanical properties of the native tissue has been a challenge since the development of the first TEVG by Weinberg and Bell. ¹⁴ Niklason and Langer ¹⁵ were among the very first groups that tried to address this issue by developing a novel bioreactor that allowed culturing TEVGs under conditions of pulsatile flow. Since then, various groups have tried to improve the mechanical strength of TEVGs by application of various loading regimens and/or manufacturing techniques.^16–18

Agent-based modeling (ABM), finite element analysis (FEA), and coupled AB-FEA have been widely used to study mechanobiological interactions in various diseased and healthy states, such as Duchenne muscular dystrophy,^19,20 volumetric muscle loss, ²¹ neovascularization in porous stents, ²² and vascular tissue engineering. ²³ ABM is a computational technique for simulation of complex adaptive systems, ²⁴ in which the overall response of the system is defined by the culmination of interactions and behaviors of individual agents (e.g., cells), and their response to external stimuli. These interactions and behaviors are defined using “rules,” which can be in the form of mathematical equations or logical expressions.

In this study, we extend a coupled AB-FEA computational modeling framework, previously developed by our group, ²⁵ to investigate the effect of various loading scenarios (uniaxial, biaxial, and equibiaxial), boundary conditions, stretch magnitudes, and growth factor concentrations on growth and remodeling of TEVGs. We set up the model based on the work of Huang et al. ¹⁶ The AB-FEA model reported the production and alignment of collagen fibers, proliferation of the smooth muscle cells (SMCs), and concentrations of platelet-derived growth factor (PDGF), transforming growth factor-β (TGF-β), and vascular endothelial growth factor (VEGF) over time. For validation of the model, we compared the collagen content, SMC density, and fiber orientation with Huang et al. ¹⁶

Methods

We used a modified version of a coupled AB-FEA computational framework, which was previously established by our group, ²⁵ to study the growth and remodeling of TEVGs under various loading scenarios. The coupled AB-FEA framework utilizes REPAST Simphony libraries ²⁴ for the ABM simulations and the commercial FEA package ANSYS (ANSYS, Inc., Canonsburg, PA). We set up the model based on the experimental conditions reported by Huang et al. ¹⁶ (inner radius = 3.2 mm, thickness = 200 μm, and initial SMC population = 2.625e5). To reduce the computational time, we simulated a 1.6 mm segment of the graft. Huang et al. used a polyglycolic acid (PGA) mesh as the scaffold. The PGA meshes were wrapped around a distensible silicone tube and a pulsatile pressure of 270 mmHg/−30 mmHg was applied to the inner surface of the tube. This loading regimen led to a circumferential strain of 1.5%–2% on the surface of the tube. In this study, we did not include the silicon tube in the FEA simulations and instead applied the reported circumferential strain to the luminal surface of the tube. The PGA mesh was also excluded from the model as reviewing prior work ²⁶ showed it did not contribute to the mechanics and is biocompatible, nontoxic, and FDA approved for human implantation. ²⁷ We investigated the effect of various loading modes (uniaxial, biaxial, or equibiaxial) and combinations of growth factors, present in the serum, on the growth and remodeling of TEVGs, with specific focus on formation of collagen fibers and their orientation. All the figures were generated in MATLAB (MathWorks, MA).

Agent-Based Model

The ABM was developed in the REPAST Simphony environment (lattice based) and includes two superclasses of agents: patches and cells. The lattice size (patch dimension) was set to 100 × 100 × 400 μm (X, Y, and Z). The cell agents interact with a virtual representation of the spatially dependent soluble factors and the physical environment. In this study, as our focus was on SMC-populated scaffolds, the only cell type present in the model were SMCs. Structurally significant extracellular matrix (ECM) proteins were collagen fibers and gelatin; elastin was excluded due to the limited quantifiable information on production of elastin by SMCs. In this model, gelatin is the product of collagen degradation and it can influence the growth and remodeling by modulating the mass balance and affecting the mechanics (gelatin's stiffness^28–30 is lower than collagen and higher than the cells, please refer to Table 2 to compare the material constants). Diffusion of the growth factors, cytokines, and matrix metalloproteinase (MMP) was handled by numerically solving the transient partial differential equations using a forward in time-centered in space algorithm in each ABM time step. The ABM ran the MATLAB Java engine in the background to perform these calculations. Each time step is a “tick” and represents 6 h of time.

Collagen fiber orientation

Cyclic mechanical stretch induced the production of collagen by SMCs.^47–49 The new collagen fibers were deposited parallel to the SMC long axis. ³⁸ We assumed that collagen fibers retain their orientation during loading. In each time step and patch, the resultant orientation of the collagen fibers was determined by calculating the weighted average of fiber orientations deposited during previous time steps. The newly deposited collagen, in each time step, was used as the weight in the weighted average. In each patch, there were two families of fibers (Fig. 1).

OPEN IN VIEWER

FIG. 1.

Schematic illustrating fiber angle definitions. (a) Fiber angle α is the angle between the collagen fiber A and the radial direction. Fiber angle β is the angle between the component of the collagen fiber A in the Z-θ plane, A_Zθ, and the θ axis. (b) Resultant fiber orientation, A_S, was determined by calculating the weighted average of the fiber orientations deposited in each time step (A₁ and A₂).

Finite Element Analysis

The interaction module, within the ABM, managed running the FEA simulations and transferring the data back and forth between the ABM and FEA. In each interaction, the FE models were created in batch mode, using an ANSYS Parametric Design Language script that built the FE model based on the point cloud (describing the geometry), point-wise mass fractions of contents (collagen, gelatin, and the cells), and fiber orientations available from the ABM's most recent tick. More details are provided in prior work. ²⁵

We used a mass fraction-based version of the anisotropic hyperelastic strain energy density (SED) function, originally proposed by Holzapfel et al., ⁶⁴ and assumed the arterial tissue to be incompressible. Noteworthy, this SED function does not account for the active mechanical properties of SMCs. Equations (2.1–3) depict the SED functions (i.e., W, w_i, and w_c) implemented in the FEA simulations, where, φ_i and ψ_c are mass fractions of the isotropic (cells and gelatin) and anisotropic contents (the two collagen fiber families), respectively. C_i is the initial shear modulus of the isotropic content i (which can be cells, elastin, or gelatin), κ₁ and κ₂ are constants, A and B are unit vectors that describe directions of collagen fiber families, and I₁ is the first invariant of the right Cauchy-Green tensor, C.

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split} W = \mathop \sum \limits_{i = 1}^2 { \O _i} \cdot {w_i} + { \psi _c} \cdot {w_c} = \mathop \sum \limits_{i = 1}^2 { \O _i} \cdot {{{c_i}} \over 2} \left( {{I_1} - 3} \right) \\+ { \psi _c} \cdot \left( {{{{ \kappa _1}} \over {2{ \kappa _2}}} \left( {{e^{{ \kappa _2}{{ \left( {{{ \bar I}_4} - 1} \right) }^2}}} - 1} \right) + {{{ \kappa _1}} \over {2{ \kappa _2}}} \left( {{e^{{ \kappa _2}{{ \left( {{{ \bar I}_6} - 1} \right) }^2}}} - 1} \right) } \right)\end{split} \tag{{\rm Equation} \ 2.1} \end{align*} \end{document}

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \bar I_4} = A \cdot \bar C \cdot A \tag{{\rm Equation} \ 2.2} \end{align*} \end{document}

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \bar I_6} = B \cdot \bar C \cdot B \tag{{\rm Equation} \ 2.3} \end{align*} \end{document}

We adapted the ANSYS AHYPER material model (Equations 2.4–7), which allows for modeling soft fiber-reinforced biological tissues, to mimic the two-fiber model suggested by Holzapfel et al. ⁶⁴ The AHYPER model describes the SED as the summation of volumetric (W_v) and isochoric (W_d) SEDs, where C is the right Cauchy Green tensor, a_i and b_j are material constants of the isotropic part of the SED, c₁ and c₂ are the material constants for the first fiber, e₁ and e₂ are the material constants for the second fiber, A and B are the unit vectors for each fiber direction, I₁ is the first invariant, and I₂ is the second invariant.

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {W_s} = {W_v} \left( J \right) + {W_d} \left( { \bar C , \;A \otimes A , \;B \otimes B} \right) \tag{{\rm Equation} \ 2.4} \end{align*} \end{document}

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { W_v } \left( J \right) = \frac { 1 } { d } \cdot { \left( { J - 1 } \right) ^2 } \tag { { \rm Equation } \ 2.5 } \end{align*} \end{document}

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split} {W_d} \left( { \bar C , \;A \otimes A , \;B \otimes B} \right) = \mathop \sum \limits_{i = 1}^3 {a_i}{ \left( {{{ \bar I}_1} - 3} \right) ^i} + \mathop \sum \limits_{j = 1}^3 {b_j}{ \left( {{{ \bar I}_2} - 3} \right) ^j} \\\quad\quad\quad\quad\quad+ {{{c_1}} \over {2{c_2}}} \left( {{e^{{c_2}{{ \left( {{{ \bar I}_4} - 1} \right) }^2}}} - 1} \right) + {{{e_1}} \over {2{e_2}}} \left( {{e^{{e_2}{{ \left( {{{ \bar I}_6} - 1} \right) }^2}}} - 1} \right)\end{split} \tag{{\rm Equation} \ 2.6} \end{align*} \end{document}

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \bar I_4} = A \cdot \bar C \cdot A \tag{{\rm Equation} \ 2.7} \end{align*} \end{document}

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \bar I_6} = B \cdot \bar C \cdot B \tag{{\rm Equation} \ 2.8} \end{align*} \end{document}

The material constants were obtained from experiments performed by Seyedsalehi et al. ⁶⁵ for collagen fibers, and Norris and McGraw ²⁹ and Qiu et al. ⁶⁶ for gelatin and SMCs, respectively (Table 2). Each patch, within the ABM, was associated with a unique set of material properties in the FEA, and we did not use a material table or clustering algorithm to combine similar materials. The boundary conditions were chosen to mimic the loading scenarios reported by Huang et al. ¹⁶

Model Setup

We constructed the model to match the geometry, initial values, and boundary conditions of the TEVG developed by Huang et al. ¹⁶ and investigated the effect of various (a) loading modes (uniaxial, biaxial, or equibiaxial), (b) boundary conditions, (c) stretch magnitudes, and (d) growth factor concentrations on the matrix content, SMC alignment, and collagen fiber orientation of the TEVG during a 12-week period. The graft was fixed at one axial end and, based on the loading scenario, was stretched (1.04 or 1.08), fixed, or free to move at the other end. A circumferential stretch (1.02, 1.04, or 1.08) was applied to the lumen. TGF-β (1.26 ng/mL) and PDGF-AA (1 ng/mL) were added to the model to represent culture media, on the abluminal side of the graft, as these growth factor concentrations are approximately equivalent to their concentrations in 10% serum. We then compared the trends in the predicted collagen content, fiber alignment, and SMC density with the values reported by Huang et al. ¹⁶ at 12 weeks.

Results

Role of loading mode

Below, we show the role the loading mode has on the composition, collagen fiber angles, and synthesis of soluble factors within the vascular graft. We set up the AB-FEA model based on initial conditions reported by Huang et al. ¹⁶ The vascular grafts were simulated under static (fixed axial ends), uniaxial (fixed axial ends, circumferential stretch = 1.02), biaxial (axial stretch = 1.08 and circumferential stretch = 1.02), and equibiaxial conditions for 12 weeks. The first three loading modes were based on conditions reported by Huang et al. ¹⁶ Equibiaxial loading was not reported on by Huang et al.

Loading mode and ECM composition

Biaxial loading (axial stretch = 1.08 and circumferential stretch = 1.02) had the highest impact on collagen content. At week 12, the collagen content was 6.49e3 ng (197.81 ng/mm³) under static, 7.51e3 ng (230.16 ng/mm³) under uniaxial, and 9.92e3 ng (301.55 ng/mm³) under biaxial conditions. Experimentally, under the same conditions, Huang et al. ¹⁶ also observed collagen content was highest under biaxial loading (static: 54% ± 7%, uniaxial: 48% ± 5%, and biaxial: 67% ± 3% of dry weights), although collagen was not significantly different between any loading condition. At week 12, SMC density was similar under static, uniaxial, and biaxial conditions (4514 Cells/patch, 4518 Cells/patch, and 4526 Cells/patch, respectively). Experimentally, at 12 weeks, Huang et al. also found SMC density was not significantly different between loading conditions (i.e., 970 ± 135 nuclei/mm² for the static, 895 ± 35 nuclei/mm² for the uniaxial, and 1256 ± 43 nuclei/mm² for the biaxial case). For direct validation of the theoretical results, we cannot perform a one-to-one comparison of the measured values since Huang et al. ¹⁶ reported collagen content as percentage of dry weight and the SMC density in units of nuclei/mm² (and not nuclei/mm³). Regardless, the AB-FEA model captured similar results to the experimental TEVG.

Loading mode and collagen fiber orientation

We categorized the fibers based on their radial positions. For each radial position, we plotted the histogram of the fiber angles [α: the angle between the Z-θ plane and radial direction, and β: the angle between Z and θ axes (Fig. 1)]. Figures 2 to 5 depict the orientations of collagen fibers at the start of the simulation and after 12 weeks of uniaxial, biaxial, and equibiaxial mechanical stimulation. Figures 6 and 7 depict the distributions of maximum principal stress and strain for the biaxial loading case.

OPEN IN VIEWER

FIG. 2.

Initial collagen fiber orientations illustrating random assignment before loading. (a) Vector plot depicts the orientations of the two fiber families at each patch in the space (only every other vector in the XY plane is shown for visualization in this and similar figures). Histogram of (b) α and (c) β collagen fiber angles grouped by radial position (±50 μm).

OPEN IN VIEWER

FIG. 3.

Under uniaxial loading (λ_θ=1.02, axially fixed), collagen fiber orientations are mainly perpendicular to the direction of stretch (β = ∼85°). (a) Vector plot of collagen fiber directions. Histograms illustrate (b) fiber orientations are primarily helical (α = 85°) and (c) aligned with the z-axis of the TEVG (β = ∼85°), since SMCs preferred to avoid the 1.02 circumferential strain and align and deposit collagen axially. SMC, smooth muscle cell; TEVG, tissue-engineered vascular graft.

OPEN IN VIEWER

FIG. 4.

Under biaxial loading (λ_Z=1.08, λ_θ=1.02), collagen fiber orientations vary from circumferential (β = ∼20°) to axial as we move from the lumen to the outer boundary of the TEVG. (a) Vector plots of collagen fiber directions. Histograms illustrate (b) fibers are mostly helical (α = 85°) and (c) fibers in the inner layers, aligned dominant β angle between 15° and 35°. In the outer layers, majority of fibers aligned in the axial direction.

OPEN IN VIEWER

FIG. 5.

Equibiaxial loading (λ_Z=1.04, λ_θ=1.04) reversed the (a) distribution of the collagen fiber orientation throughout the thickness of the TEVGs (compared to biaxial loading). Histograms illustrate (b) fibers are mostly helical (α = 85°) and (c) fibers in the inner layers were primarily oriented at β = 55°. In the middle layers, three distinct fiber orientations were observed (β = 35°, 45°, and 55°). In the outer layer, there was no dominant orientation.

OPEN IN VIEWER

FIG. 6.

Average maximum principal stress reaches the steady state (maximum) after 4 weeks of biaxial loading (λ_Z=1.08 and λ_θ=1.02). (a) Distribution of maximum principal stress at week 12. (b) Normalized average maximum principal stress throughout simulation. (c) Histogram of the distribution of maximum principal stress at 12 weeks.

OPEN IN VIEWER

FIG. 7.

Average maximum principal strain approaches steady state (minimum) after 4 weeks of biaxial loading (λ_Z=1.08 and λ_θ=1.02). (a) Distribution of maximum principal strain at week 12. (b) Normalized average maximum principal strain throughout simulation. (c) Histogram of the distribution of the maximum principal strain values for all patches at 12 weeks.

In the absence of mechanical stimuli, whether it was at the start of the simulations or under static conditions, collagen fibers were randomly aligned (Fig. 2). Under uniaxial circumferential loading (Fig. 3), the majority of the collagen fibers aligned perpendicular to the direction of strain (β = ∼85°, Fig. 3c) and avoided the radial direction (α = ∼85°, Fig. 3b). A fiber's unit vector does not have any component in the radial direction when α = 90°; therefore, the fibers with α values close to 90° demonstrate a helical structure, which is consistent with vessels observed in vivo. In contrast, under biaxial loading, collagen fiber orientation varied from circumferential (β = 25°, Fig. 4c) in the lumen to mostly axial in the outer layers. These predictions are similar to what is observed in vivo^67,68 and the orientations measured by Weidenhamer and Tranquillo ⁶⁹ Moreover, qualitative results from all simulated cases (static, uniaxial, and biaxial) agree with the general trend observed experimentally by Huang et al. ¹⁶ Specifically, Huang et al. ¹⁶ also observed at 12 weeks under static loading conditions random orientation of collagen fibers throughout the artery wall; under uniaxial conditions, collagen fibers axially oriented in the inner layers and circumferentially aligned in the outer, and under biaxial conditions, the reverse, circumferentially aligned in the inner and axially oriented in the outer layers. Uniquely, under equibiaxial loading, the physiological trend that was observed in the biaxial case reversed. The β angle was highest in the inner layers and gradually reduced throughout the thickness. In the middle layers, three distinct fiber orientations were observed (β = 35°, 45°, and 55°) and interestingly, there was no preferred direction in the outer layer.

Role of Axial Boundary Conditions

To investigate the effect of axial displacement on collagen deposition and orientation, we compared the axially fixed and free TEVGs under uniaxial circumferential loads. The grafts were either fixed at both ends (0 displacement constraints in the Z direction) or fixed only at one axial end and simulated under 1.02 and 1.04 circumferential stretch for 12 weeks. In contrast with the axially fixed TEVGs (Fig. 3), collagen fibers aligned in the circumferential direction in the axially free TEVGs (Fig. 8). This phenomenon, which was also observed in vitro by Weidenhamer and Tranquillo, ⁶⁹ can be explained by the preference of cells to avoid stretch and compaction.^70,71 Moreover, when axial displacement was allowed, fiber orientation was not a function of stretch magnitude (no difference between 1.02 and 1.04 cases), which is in contrast with the axially fixed cases (Supplementary Figs. S1 and S2). There was no difference in collagen mass and SMC population between the axially fixed and free cases (Table 3, Supplementary Figs. S3 and S4).

OPEN IN VIEWER

FIG. 8.

Under circumferential loading and allowed axial displacement, collagen fibers align in the circumferential direction. TEVGs were only fixed in one axial end and stretched circumferentially (2%) for 12 simulated weeks. (a) Vector plot of fiber directions. In contrast with TEVGs fixed at both axial ends (Fig. 3), collagen fibers aligned in the circumferential direction. Histograms illustrate (b) fibers are mostly helical (α = 85°) and (c) all fibers were aligned at β = 5°.

Role of Stretch Magnitude

Increasing the stretch magnitude (circumferential, axial, or both) increased the total collagen mass, but did not affect the SMC population (Table 3). Axial stretch (λ_Z) was unstrained, 1, 1.04, or 1.08 and the circumferential stretch (λ_θ) was set to 1.02, 1.04, or 1.08. When both axial ends were fixed, increasing the circumferential stretch (λ_θ) increased the percentage of fibers that were aligned circumferentially in the outer layers; the majority of fibers, however, were still aligned along the axial directions (please refer to Supplementary Fig. S5). In the biaxial (and equibiaxial) loading modes, the fiber alignment did not depend on the stretch magnitude, but rather upon the ratio of axial and circumferential stretches (please refer to Fig. 5 and Supplementary Fig. S6 to compare the 1.04 and 1.08 equibiaxial loading). SMC population was not affected by the stretch magnitude.

Growth Factor Concentration

We investigated the effect of various combinations of growth factor (TGF-β and PDGF-AA) concentrations in the culture media on growth and remodeling of the biaxially loaded TEVGs. These combinations are shown in Table 4; rows 1–3 represent serum concentrations of 1%, 5%, or 10%. The grafts were simulated for 12 weeks under biaxial loading (λ_Z = 1.08, λ_θ = 1.02). Decreasing serum concentration from 10% (TGF-β = 1.26 ng/mL and PDGF-AA = 1 ng/mL, row 3) to 5% (row 2) or 1% (row 1), significantly slowed down the growth of the TEVGs. At week 12, the TEVGs with 10% serum reached a wall thickness of ∼1 mm, while the TEVG thickness for the 5% and 1% serum was 0.441 mm and 0.179 mm, respectively. Collagen fiber orientation was similar under 5% and 10% serum throughout the thickness. However, fibers were more randomly aligned in the Z-θ plane under 1% serum.

Discussion

Despite recent advances, developing a TEVG that mimics the performance of the native tissue is still a challenging task, specifically with regard to the long term mechanical performance, as no single design has yet been able to achieve—and retain—the desired mechanical strength and compliance. ² This, in part, is due to the intrinsic complexity of the mechanobiological interactions involved in the growth and remodeling of the vascular tissue, which pose a challenge to engineering the grafts. In this study, we presented a multiscale computational modeling framework that can facilitate the development of vascular grafts and studied the effect of various loading modes (uniaxial, biaxial, and equibiaxial), boundary conditions, stretch magnitudes, and growth factor concentrations on ECM composition and orientation of collagen fibers in TEVGs.

Biaxial loading—and not equibiaxial or uniaxial—with dominant stretch in the axial direction led to physiological fiber alignment throughout the thickness of the TEVG. Holzapfel reported the fiber orientations in the intima, media, and adventitia of the abdominal aorta to be 18.8°, 37.8°, and 58.9°, respectively. ⁶⁸ Our model predicted a similar trend for the biaxial loading mode (λ_Z = 1.08 and λ_θ = 1.02) and showed that under uniaxial circumferential loading, axial boundary conditions play a key role in determining the SMC orientation and therefore, collagen fiber alignment (Figs. 3, 8, and 9). This trend was in accord with the observations of Foolen et al. ⁷² They reported that after 3 days of uniaxial cyclic stretch, stress fibers aligned in the direction of the constraints. Weidenhamer and Tranquillo ⁶⁹ also reported a similar alignment of collagen fibers in TEVGs. In that study, under uniaxial circumferential stretch, collagen fiber orientation was dependent on the axial boundary conditions. When the scaffold was fixed in the axial direction, fibers aligned in the axial direction, and in contrast, when scaffolds were allowed to compact in the axial direction, in response to increased circumferential stretch, fibers aligned in the circumferential direction. In addition, de Jonge et al. reported that cells and collagen fibers realign to the direction of uniaxial constraints. ⁷³

OPEN IN VIEWER

FIG. 9.

Under circumferential loading and allowed axial displacement, SMCs align in the circumferential direction. TEVGs were only fixed in one axial end and stretched circumferentially (2%) for 12 simulated weeks. (a) Vector plot of SMC orientation in each patch in the XY plane. Histograms of (b) the α college fiber angle suggests that SMCs avoided the radial direction and (c) the β fiber angle suggests that most SMCs align circumferentially throughout the TEVG thickness.

Mechanical loading increased the collagen content of simulated vascular grafts in a stretch magnitude-dependent manner (Table 3); the SMC number, however, was not affected by mechanical stretch. In this study, we did not implement any ABM rule that directly affects—or describes—the proliferation rate of SMCs as a function of stretch, as there was no consensus in the literature on the inhibitory or augmenting effect of stretch on proliferation rate. While Chapman et al., ³³ Sumpio and Banes, ⁷⁴ and Hipper and Isenberg ⁷⁵ have reported an inhibitory effect for stretch, Mills et al., ⁷⁶ Wilson et al., ⁷⁷ and Iwasaki et al. ⁷⁸ have reported the SMC proliferation rate to increase in response to cyclic stretch. Wilson et al. ⁷⁷ attributed the stretch-induced SMC proliferation to increased expression—and autocrine action—of PDGF-BB by de novo SMCs; however, other researchers have reported that vascular SMCs do not produce PDGF-BB.^31,32

In the 12-week simulation using biaxial loading (λ_Z=1.08, λ_θ=1.02), a decrease in the provided serum concentrations of TGF-β and PDGF-AA (10% serum to 1% or 5% serum) substantially affected the population of SMCs (59.65% decrease for 5% serum and 96.90% decrease for 1% serum), the total collagen mass (60.75% decrease for 5% serum and 95.86% decrease for 1% serum), and the average thickness (55.16% decrease for 5% serum and 81.87% decrease for 1% serum), and consequently slowed down the growth and remodeling of the graft. Increasing TGF-β concentration (from 1.26 ng/mL to 12.6, 25.2, or 50.4 ng/mL) alone, while PDGF-AA concentration was kept constant at 1 ng/mL, either did not affect (TGF-β = 12.6 ng/mL) or slightly increased (TGF-β = 25.2 and 50.4 ng/mL) the total collagen mass. The population of SMCs slightly decreased (TGF-β = 12.6 and 25.2 ng/mL) and then increased (TGF-β = 50.4 ng/mL). Interestingly, the collagen mass per cell was highest at the lowest concentration of TGF-β (Table 4, row 1). In this case, the average thickness of the graft, at tick 336, was 179.54 μm. Further investigation showed that the maximum principal strain, which stimulates the production of TGF-β by SMCs, drops after the graft thickness reaches ∼200 μm (Fig. 6b, day 28); therefore, SMCs that were located beyond that thickness had access to less TGF-β and produced less collagen. This observation suggests that it might be necessary to, gradually, increase the stretch magnitude, as the thickness of the graft increases (implementing a time-dependent loading profile) to achieve the desired mechanical properties.

In conclusion, growth and remodeling of SMC-populated TEVGs without topographic cues under various loading conditions were presented using our coupled AB-FEA modeling framework. Through investigation of the evolution of collagen fiber orientations, collagen content, and SMC count during 12 weeks of culture under various loading scenarios, it is evident that biaxial (axially dominant) loading is necessary to achieve the physiological distribution of fiber orientations observed in vivo.