|Home | About | Journals | Submit | Contact Us | Français|
Recently a cartilage growth finite element model (CGFEM) was developed to solve non-homogeneous and time-dependent growth boundary value problems [1, 2]. The CGFEM allows distinct stress constitutive equations and growth laws for the major components of the solid matrix, collagens and proteoglycans. The objective of the current work was to simulate in vitro growth of articular cartilage explants in a steady-state permeation bioreactor in order to obtain results that aid experimental design. The steady-state permeation protocol induces different types of mechanical stimuli: when the specimen is initially homogeneous it directly induces homogeneous permeation velocities and indirectly induces non-homogeneous solid matrix shear stresses; consequently, the steady-state permeation protocol is a good candidate for exploring two competing hypotheses for the growth laws. The analysis protocols were implemented through the alternating interaction of the two CGFEM components: poroelastic FEA using ABAQUS and a finite element growth routine using MATLAB. The CGFEM simulated 12 days of growth for immature bovine articular cartilage explants subjected to two competing hypotheses for the growth laws: one that is triggered by permeation velocity and the other by maximum shear stress. The results provide predictions for geometric, biomechanical, and biochemical parameters of grown tissue specimens that may be experimentally measured and, consequently, suggest key biomechanical measures to analyze as pilot experiments are performed. The combined approach of CGFEM analysis and pilot experiments may lead to the refinement of actual experimental protocols and a better understanding of in vitro growth of articular cartilage.
The avascularity and other intrinsic properties of articular cartilage (AC) preclude its ability to heal following damage or degeneration . Possible repair strategies include the replacement of damaged AC with tissue engineered constructs grown in vitro  or with graft tissue [5, 6]. The long-term goal of this work is to develop continuum mechanics models of AC growth that may be used with experimental approaches to improve AC repair strategies. AC often is modeled as poroelastic material with solid matrix (SM) and fluid constituents . The SM is composed of cells (chondrocytes) surrounded by collagen (COL) fibers, that resist tensile and shear loads, and proteoglycans (PGs), that resist compressive loads due to the fixed charge density of their glycosaminoglycan (GAG) components. The SM behaves as a fiber-reinforced composite with anisotropic and strain-dependent mechanical properties that exhibit strong tension-compression asymmetry .
In recent years there has been growing interest in the development of tissue growth models, including models of transport and metabolic processes and/or matrix accumulation in native AC tissue  and engineered constructs [10-12], continuum growth mechanics with small strains for native AC tissue  and engineered constructs , and general continuum mechanics growth models with large strains [15-18]. In this paper, growth of native AC tissue is studied which presents challenges different than those that model tissue engineered constructs, as the matrix is more mature and, consequently, modeling large deformations with more complex stress-strain equations becomes more relevant. Although the transport and metabolic processes are not considered in detail here, these different types of models may be merged as more data becomes available as discussed previously .
A cartilage growth finite element model (CGFEM) was developed to solve non-homogeneous and time-dependent growth boundary-value problems (BVPs) [1, 2]. That approach integrated the poroelastic FEA solver of ABAQUS with custom MATLAB routines that solved the incremental growth BVP  for each finite element. Although the incremental growth BVP is limited to homogeneous growth, non-homogenous growth of the total specimen is allowed because each finite element can experience different homogeneous growth. The CGFEM may be used to provide an a priori prediction of in vitro growth experiments that can be used to guide experimental design and reduce the number of time-consuming and costly pilot experiments.
AC exhibits depth-dependent biochemical and biomechanical properties that are correlated with one another . In order to enhance the depth-dependence of biomechanical properties of immature tissue grown for cartilage repair, it may be necessary to apply mechanical loading in vitro that produces depth-dependent mechanical stimuli. The CGFEM was used to study the evolution of depth-dependent properties using a confined compression protocol ; however, growth in a permeation bioreactor is more common.
The steady-state permeation BVP was solved using large deformations of a biphasic mixture with an isotropic SM in  and the dynamic permeation BVP was solved using infinitesimal deformations of a triphasic mixture with an isotropic SM in . The results of  predicted that the permeation velocity is homogeneous while the SM strain, stress, and fluid pore pressure are non-homogeneous. The steady-state permeation protocol induces different types of mechanical stimuli: when the specimen is initially homogeneous it directly induces homogeneous permeation velocities and indirectly induces non-homogeneous SM shear stresses; consequently, the steady-state permeation protocol is a good candidate for exploring two competing hypotheses for the growth laws. If the growth laws depend explicitly on permeation velocity as suggested for AC explants in [24, 25], then homogeneous growth may occur. Indeed, GAG and COL synthesis and accumulation has been positively correlated with permeation velocity for tissue engineered constructs grown in permeation bioreactors [26, 27]. However, if the growth laws depend explicitly on SM shear stress as suggested in  and indirectly induced by permeation, then non-homogeneous growth may occur.
Nearly all of the articular cartilage growth experiments in a permeation bioreactor involve tissue engineered constructs; one study reporting long-term culture in a permeation bioreactor with explanted human cartilage  does not discuss several of the unique difficulties posed by using more mature explanted tissue. For explanted tissue that will have a substantially lower permeability than tissue engineered constructs, these difficulties include the large pressure gradients and associated matrix compaction that may affect cell viability, and the need for substantial pre-compression strains to prevent platen-specimen separation during growth which would be difficult to quantify experimentally.
The objective of this work was to use the CGFEM with a finite strain anisotropic SM based on a PG-COL stress balance to predict how immature AC explants will evolve during growth in a steady state permeation bioreactor, thereby providing extensions to earlier analyses of articular cartilage subjected to permeation [22, 23]. The results reported here correspond to the first stage corresponding to the long-term goal of growing and shaping constructs to have the desired geometric and biomechanical properties needed for more successful clinical outcomes. Subsequent stages involve using these results to design experiments including biomechanical and biochemical testing before and after growth, followed by iterations on the analyses and experiments.
The specific aims were to: 1) modify the CGFEM to incorporate new PG and COL stress constitutive equations for finite deformations; 2) solve the permeation boundary-value problem using poroelastic FEA; 3) solve the total specimen growth BVP for two competing hypotheses for the growth laws; and 4) simulate post-growth experiments needed to quantify the biomechanical properties and to estimate the growth laws.
Consider a continuous body B that occupies a reference configuration Κ0 (B) at a time t0. At time t, B occupies a new configuration Κ(B). A point P on body B has a reference position vector X at t0 and then later a position vector x at t. The motion of B is therefore described by x = χ(X,t). The deformation gradient tensor that describes this motion is F = χ(X,t)/X. The Jacobian J = detF describes the volumetric change in B. The (absolute) velocity vector v is given by the material time derivative of x = χ(X,t).
The theory in this section is equivalent to [22, 30]. The solid and fluid constituents are denoted with the superscripts s and f, respectively, whereas an arbitrary constituent is denoted with the superscript α. Constituent volume fractions ϕα are defined as ϕα = ρα / ραT where ρα is the apparent density (mass/tissue volume) and ραT is the true density (mass/constituent volume). The intrinsic incompressibility constraint relates ϕs and ϕf:
Using (1), the continuity equation for biphasic mixtures reduces to
where vα is the constituent absolute velocity vector and div is the divergence operator. The permeation velocity (i.e. relative or effective pore pressure fluid velocity) is
which reduces to vf during stead-steady permeation because vs = 0. The total mixture Cauchy stress T is decomposed as
where Tα are constituent partial Cauchy stresses defined per unit tissue area. The void ratio e is
so that the strain-dependent permeability k (units m2/Pa·s) from  can be expressed in terms of void ratio as
where k0 and e0 are the initial permeability and void ratio, respectively, and M is a non-dimensional material constant.
The theory in this section is from [8, 31-33] and describes the incremental growth BVP. The analysis is limited to pre- and post-growth equilibrium configurations (Κr and Κg, respectively) of an unloaded homogeneous SM element undergoing homogeneous growth (Fig. 1) because this solution is obtained for each finite element in the CGFEM (see section 2.4). The SM element is modeled as a mixture of two growing elastic materials, PG and COL, with distinct mechanical properties and growth laws. The superscripts p and c will be used to designate the PG and COL constituents, respectively. The deformation gradient tensor F describes the evolution of the SM stress-free configuration due to growth. The immobility constraint states that all PG and COL molecules are bound to the SM, so that their deformation gradient tensors Fp and Fc equal F. Also, Fp and Fc are decomposed into growth tensors that describe mass deposition (or removal) and elastic growth tensors that ensure continuity of the growing SM element. Using the immobility constraint one obtains
A crucial assumption, which has been elaborated upon in [8, 20, 31, 34], leads to experimental prescriptions for and : the constituent density and stress functions are independent of . Consequently, the continuity equations become decoupled into elastic and growth continuity equations
respectively, where the mass growth functions cα are the net rate of mass deposition per unit current mass (s-1), ρα and are the apparent densities in Κg and Κr, respectively, , , and det is the determinant operator. When growth (i.e. cα) is homogeneous and constant during a time step Δt, (8)2 becomes
where and Δmα are the original mass and change in mass of a constituent, respectively. Consequently, cα and can be computed from experimental mass data.
Since the SM element is homogeneous and unloaded in Κr and Κg, the traction-free boundary condition on all surfaces and the equations of motion at equilibrium (i.e. divTs = 0) are satisfied if
where Ts satisfies the stress balance hypothesis
and Tp and Tc are the PG and COL partial Cauchy stresses, respectively.
General stress constitutive equations are defined in terms of total elastic tensors and relative to distinct PG and COL reference configurations and :
The stress constitutive theory used in this paper was developed in [35, 36] and is presented in the Appendix. It is emphasized that the manner in which evolve during a computational growth solution, presented below, is not trivial.
The theory in this section is from [2, 20]. A fundamental assumption in the CGFEM is that the time scale for the mechanical effects due to mass deposition (i.e. days) is several orders of magnitude greater than the time scale for the mechanical effects due to applied biomechanical loading (i.e. seconds) (Fig. 2). Consequently, the CGFEM has two computational components that work interactively in an incremental fashion to solve the total specimen growth boundary value problem: a total specimen poroelastic finite element analysis (FEA) component using ABAQUS and a finite element growth routine (EGR) component using MATLAB (Fig. 2).
Initially, the stress constitutive equations for the PG and COL relative to and , respectively, are defined (see Eqn. (12)) and implemented using the UMAT subroutine in ABAQUS. At the beginning of a general increment n of growth, the biomechanical factors affecting growth have been determined for every SM finite element using the poroelastic FEA component and applied biomechanical loads as described in section 2.6. Then, each finite element from the total specimen configuration nΚR is individually unloaded in the EGR to the SM element stress-free configuration nΚr through the deformation nFu (Fig. 2). Initially, the element is assumed to be stress-free so that 1Fu = I for the initial increment. For future increments the EGR computes nFu from the history of the past increments. Following unloading, the incremental growth BVP corresponding to the theory presented in section 2.2 is solved for each element. The element undergoes a deformation nδF from the stress-free configuration nΚr to the grown stress-free configuration nΚg:
where nδF and are unknown while are prescribed by the growth laws. It is emphasized that (13) has not been linearized; the incremental nature of the solution is achieved only by approximating as shown in section 2.5. The EGR solves (13) for the elastic growth tensors in terms of the other quantities:
Next, the constituent stresses in configuration nΚg are expressed using their stress functions relative to nΚr. Noting that the total elastic tensors are updated as shown below in Eqns. (16-17), the EGR uses Eqns. (10-12) to solve for the incremental growth deformation nδF and, consequently, using (14).
Once the CGFEM determines the grown, yet still stress-free, SM configuration nΚg, the element is deformed to the stressed post-growth configuration nΚR through the deformation [(nδF)(nFu)]-1. The pre- and post-growth configurations nΚR have equivalent material point locations, but different SM stress tensors and stress constitutive equations.
The CGFEM carries out the EGR for each element of the entire specimen. Although each element grows homogeneously, the growth of elements throughout the specimen may be non-homogeneous. When this occurs, the total specimen is not in equilibrium in the post-growth configuration nΚR. In the final step to solving the total specimen growth BVP, the CGFEM uses the elastic FEA component to solve for the non-homogeneous deformation nFr that ensures total specimen equilibrium in the grown configuration nΚG that will be residually stressed if non-homogeneous growth occurs. The EGR executes subsequent increments the same as the first with the exception of the initial unloading deformation.
As stated before, the EGR determines the unloading for the next increment n+1 (n>1) using the deformations from previous increments according to:
because the grown SM stress-free configuration nΚg is the stress-free configuration of the SM element needed for increment n+1 (i.e. n+1Κr) and the mapping in (15) is the reverse mapping of the SM element from nΚG to nΚg.
Finally, the stress constitutive equations need to be updated because only the elastic deformations during increment n change the constitutive stresses. For increment n=1, recall that each element in 1ΚR is assumed stress-free; thus, 1Fu = I and 1ΚR=1Κr is stress-free. The COL reference configuration is assumed to have zero COL stress (i.e. Tc =0) while the PG reference configuration is the same as 1Κr (i.e. ) and the PG has an initial swelling stress (i.e. Tp ≠ 0). Thus, a coupled set of nonlinear algebraic equations is solved to determine the initial COL swelling strains (i.e. ) that satisfy Ts =0. Then, using a theorem from , the constituent stress equations defined for the elastic deformation relative to 1Κr are derived from the general forms (12) as
After a general increment, the deformations and constitutive stress equations are updated according to
The unstimulated growth rates for cα are based on calculations in  using experimental data of calf patellofemoral groove (PFG) cartilage from 0.4-0.9 mm deep grown without mechanical stimulation for 13 days in a solution of DMEM, 20% fetal bovine serum (FBS), and 100 μg/ml ascorbate . Thus, the PG and COL unstimulated growth rates used here are 2.68%/day and 0.43%/day, respectively (Table 1).
Two different dose-independent growth triggers are used. Steady-state permeation directly induces a homogeneous permeation velocity; for the permeation trigger growth is stimulated by the magnitude of the permeation velocity at the element centroid. The trigger magnitude is 0.25μm/s which when exceeded increases PG synthesis by 50% based on observations in . In that study, experimental and theoretical results for calf PFG AC explants subjected to dynamic unconfined compression suggested that biosynthesis was stimulated in an on-off nature when the threshold fluid velocity value was exceeded and the authors noted that this observation was generally consistent with several previous results [38, 39]. COL growth was not measured in that study; thus, measurements of hydroxyproline increase in chondrocyte-seeded constructs cultured in a perfusion bioreactor  and the conversion factor for hydroxyproline to COL mass  were used to estimate a stimulated rate increase for COL that is approximately half the stimulated increase in PG and, consequently, 25% for this study. The stimulated growth rates are 4.02%/day and 0.53%/day for PG and COL, respectively (Table 1). Since the tissue is subjected to a permeation velocity double the magnitude of the trigger value, it is expected that stimulated growth will occur in every element rendering homogeneous growth; thus, the homogeneous initial specimen will remain homogeneous throughout the growth simulation.
Steady-state permeation indirectly induces a non-homogeneous SM shear stress; for the shear trigger growth is stimulated by the magnitude of the maximum SM shear (i.e. Tresca) stress at the element centroid. The shear trigger is based on a study quantifying cartilage growth in simple shear  in which calf PFG cartilage was subjected to an offset 10% compression with superposed shear strains of 1% and 3%. The results of that study suggested a trend between tissue shear and PG and COL stimulation and, since there was not a significant difference between stimulated rates at 1% and 3% shear strain levels in that study, here a trigger with an “on-off” nature is assumed. Based on those results, in this study stimulated PG and COL rate increases of 25% and 38%, respectively, are used. The stimulated growth rates are 3.35%/day and 0.59%/day for PG and COL, respectively (Table 1). However, the trigger value needed to be computed for the present study. The Tresca trigger value was estimated to be 0.07 MPa by simulating the experiment in  using ABAQUS and the constitutive models used in this study.
Five specimens were simulated. The first specimen immediately undergoes the post-growth measurement process on day 0 as a verification step. The two growth triggers are applied to specimens grown in the permeation bioreactor with growth terminated on day 6 or day 12 before post-growth testing. All specimens have identical day 0 biomechanical properties corresponding to calf PFG cartilage harvested from a depth of ~ 0.65 - 0.90 mm. Using true densities of the constituents  and constituent contents for this tissue source , the initial volume fractions were calculated (Table 2). The PG material constants were determined from this data using an intra- and extra-fibrillar water compartmental model (see Eqn. A2)  (Table 3). The COL bimodular parameters (see Eqn. A6) were obtained from a fit to experimental data performed in  and the shear modulus (see Eqn. A4) was measured previously for a deeper region  (Table 3). Permeability coefficients k0 = 2.46 × 10-15 m2/Pa-s and M=8.88 were chosen from experimental measurements for this tissue source .
The analysis was developed in parallel with the design of pilot experiments. A cylindrical disc (3.2 mm diameter, 0.29 mm thick) of calf PFG cartilage is harvested at a depth of ~ 0.75mm with the test axis in the 3-dir. (i.e. perpendicular to the articular surface). Using an aspect ratio (diameter/height) of 2.5, convergence studies suggested that a model 12 elements thick was needed; the final mesh used 7812 poroelastic 8-node linear brick elements. The time increment Δt was chosen to be 1 day based on convergence studies here and in an earlier study .
The first analysis step is to apply an initial compression to prevent separation of the tissue from the porous platens during permeation; preliminary FEA indicated pre-compression to 0.25mm (14%) is sufficient. In the apparatus this is accomplished by placing the disc into the bioreactor chamber between porous platens (Fig. 3A), which compress the specimen as the outer chamber tightens against the inner chamber (Fig. 3B). One quarter of the disc was modeled and displacement boundary conditions were applied to the resulting faces on the symmetry planes, the outer radial surface, and the drain side (Fig. 4). The supply side surface was displaced from z = 0.29mm to z = 0.25mm in one static analysis increment in ABAQUS. These boundary conditions are not changed until the analysis step that releases the disc from the confined compression (CC) chamber following the termination of growth.
Recall that the CGFEM analysis is split into FEA and EGR components. For each increment n, the FEA is performed in three steps. First, a static step solves for a new tissue equilibrium stress state after previous growth, if any. Second, a poroelastic step ramps the permeation velocity from zero to 0.50 μm/sec over a 500 second period. Third, a poroelastic step allows the tissue to reach equilibrium (i.e. steady-state permeation). During these steps, the displacement boundary conditions are the same as for the initial compression and the effects of growth are neglected. In the poroelastic steps, additional boundary conditions include specifying the permeation velocity and the pore pressure on the drain side, which is set to zero to simulate atmospheric pressure. These additional boundary conditions differ from those used in previous analyses of steady-state  and dynamic  permeation analyses, in which the pressure differential was specified and the permeation velocity was solved for; it is important to note that continuous fluid pressure at the boundaries are implied in both our and those previous studies.
Steady-state permeation output parameters of interest include the SM stress (to determine the growth laws for the shear trigger), SM strain (to monitor the presence of high compression strains that could compromise cell viability), and pore pressure (to monitor pressure gradients to aid hardware design) at the element centroids. Using the FEA results and the assumed growth laws of section 2.5, the EGR then solves the incremental growth BVP outlined in section 2.4 to determine the growth of each tissue element, the new SM stress constitutive equation relative to the stress-free grown element, and the total specimen post-growth configuration. This post-growth configuration represents a new structure for the total specimen for the following increment. The analysis is repeated for increment n+1 until growth is terminated.
Upon termination of growth, four post-growth analysis steps are conducted that simulate experiments that may quantify the evolved biomechanics due to growth. In the first post-growth step the grown specimen is transferred from the permeation bioreactor to a CC chamber for mechanical testing. The transfer step requires that the bioreactor and CC chamber be aligned so that the specimen does not expand in the radial direction if it has developed residual stress (Fig. 5A). The corresponding FEA step models this transfer by changing the boundary conditions on the specimen and allowing the model to reach equilibrium: in particular, the radial displacement boundary conditions are maintained while allowing the specimen to expand axially. The results of this simulation are used to compute the average grown height from the node coordinates (Fig. 5B) and compression offsets for the subsequent CC analysis.
In the second post-growth step the CC analysis is performed according to established protocols . In that protocol the specimen is placed between two porous platens and sequentially ramped to 15%, 30%, and 45% strains (Fig. 5C), allowed to reach equilibrium, and subjected to a series of oscillatory strains. This protocol allows the measurement of the strain-dependent secant aggregate modulus (HA) and permeability coefficients (k0, M) . The corresponding FEA step models the CC analysis by applying appropriate boundary conditions on the specimen and is simplified by considering only the equilibrium points at each strain offset and ignoring time-dependent poroelastic effects. The element mesh is analyzed in three layers of equal numbers of elements thick: supply side (S), middle (M), and drain side (D) slices. In MATLAB, the FEA results are post-processed: the stress predictions are used to calculate HA at the three strain offsets for each slice as well as the total specimen. Also, average void ratios for each slice and the total specimen are used to calculate the permeability at each strain offset with Eqn. (6) from which (k0, M) are re-calculated relative to the grown configuration.
In the third post-growth experiment the specimen is released from the CC chamber and undergoes free-swelling to equilibrium followed by geometric measurements of height, diameter, edge lift, and wet weight (WW) (Fig. 6A). The FEA is performed as a static step: all of the restraining boundary conditions are released allowing the disc to expand to a new free-swelling height and diameter, while allowing for curvature. The nodal coordinates and new biochemical composition of each finite element are post-processed in MATLAB to calculate total specimen averaged values of PG and COL volume fractions, tissue WW, free-swelling diameter and height, and edge lift. The edge lift is defined as the gap between a flat surface and the edge of the disc that may be measured optically and may be used to calculate the curvature.
In the fourth post-growth experiment the disc is cut into the S, M and D slices, each slice is weighed wet to obtain WW, and assayed to determine PG (measured as glycosaminoglycans (GAG)) and COL masses (Fig. 6B). To perform the analysis, the FEA results are passed to MATLAB, which is used to quantify the individual slices that are obtained from the disc slicing procedure. As before, three layers of equal numbers of elements thick (S, M, and D slices) are used. The GAG and COL masses and tissue WW are averaged for each slice.
Simulation of the permeation experiment confirmed that the supply side tissue face remained in contact with the porous platen. For every day of growth for both permeation and shear triggers, the 14% initial pre-compression was enough to maintain a small compressive stress at the surface. The evolution of the steady-state permeation solution during growth was similar for both the permeation and shear triggers; only the results for the shear trigger are presented (Fig. 7). Substantial pore pressure gradients develop on the drain side of the specimen (Fig. 7A) as the stresses induced by permeation during the ramp period compress the drain side of the tissue and, consequently, lower the permeability in this region. The pressure gradient further increases on the drain side during growth as the permeability further decreases due to mass deposition leading to decreased porosity. The pore pressure drop across the specimen on day 0 was 0.29 MPa (41 psi) and on day 12 was 0.54 MPa (78 psi) for the shear trigger, occurring mostly on the drain side (z=0) of the specimen. Although the radial, circumferential, and axial stresses all increase during growth (results not shown) the Tresca shear stress profiles do not change substantially (Fig. 7B).
The following results are limited to day 12 values because these results suggest that the actual experiments should be conducted for a minimum of 12 days to produce measurable changes. The simulations predicted homogeneous biochemical and mechanical properties for the permeation trigger and non-homogeneous properties for the shear trigger, as expected. In the former case, growth was stimulated in every tissue element because the permeation trigger is homogeneous in steady state permeation; in the latter case, stimulated growth occurred only on the drain side because of the spatial profile of the shear trigger. More specifically, for the shear stress trigger there is a stationary “growth” boundary such that only the three rows of elements adjacent to the drain side experience growth as the shear trigger value is exceeded, in contrast to the moving growth boundary simulated during confined compression loading with a permeation growth trigger .
Simulations of the CC experiments predicted that the HA increases and the k0 decreases (Fig. 8) during growth. For the shear trigger, the drain side slice had the highest HA while the middle and supply side HA values were equal. In each case, the minimum HA occurred at 30% strain; this results typifies the stress-softening phenomenon that has been measured for this tissue site . Not shown are the permeability coefficient M values of 8.69 and 8.74 for the permeation and shear triggers, respectively, and M values of 8.78, 8.78, and 8.67 for the shear-S, shear-M, and shear-D slices, respectively.
Simulations of the release experiments predicted increases in free-swelling height and diameter after growth (Fig. 9). The permeation trigger exhibited larger increases in height and diameter than the shear trigger. The curvature of the shear trigger specimen increased, as measured by edge lift, due to the increase in heterogeneous properties during growth. In particular, the edges curved towards the supply side.
GAG and COL masses for both cases were plotted normalized by initial tissue wet weight (WWI) to represent content and by final tissue wet weight (WWF) to represent concentration (Fig. 10). For both triggers, the changes in GAG content (GAG/WWI) were greater than the corresponding changes in COL content (COL/WWI), as expected since the prescribed growth rate for PG was greater. For the shear trigger the final GAG and COL contents were 6.6% and 2.6% higher in the drain side slice than the supply side slice. Although COL mass was deposited, the overall free-swelling volume of the tissue increased enough to cause the COL concentration (COL/WWF) to decrease. The dominant growth of the PG constituent contributed the most to the increase in tissue WW so the increase in GAG concentration (GAG/WWF) was positive. Contours of PG (Fig. 11) and COL (Fig. 12) volume fractions normalized by day 0 values showed similar trends.
This study used the CGFEM to simulate growth of AC explants in a steady-state permeation bioreactor as well as post-growth experiments needed to quantify the evolving biomechanical properties and to test competing hypotheses for the growth laws. Commonly, tissue engineering experiments which are both time consuming and expensive have been conducted in the absence of simulations that attempt to predict their outcomes; a contribution of this work is to illustrate how tissue growth models can be used to refine experimental protocols before the experiments are conducted. This work was accomplished with experimental design as a foremost concern, which helped to guide the FEA protocols. In turn, the CGFEM results can now be used to refine the experimental protocols. One conclusion is that post-growth experiments may involve four steps: 1) direct transfer to a CC chamber without relieving residual stress; 2) CC testing and measurement of strain-dependent aggregate modulus and permeability; 3) release from the CC chamber and geometric measurements; and 4) biochemistry measurements from multiple slices.
Several lessons were learned during the steady-state permeation analysis that affects experimental design. Permeation velocities above 0.5 μm/s cause large SM strains, a low permeability, and, consequently, numerical convergence problems (discussed below). These high strains and their associated stresses are undesirable in real experiments and may lead to cell death . Pilot experiments should be conducted with permeation velocities in the range of 0.25-0.50 μm/s while minimizing the initial compression in the SSP bioreactor needed to avoid platen-tissue separation.
Two growth triggers were considered with equal unstimulated growth rates but distinct stimulated growth rates. The geometry changes in the specimens may imply the dominance of one growth trigger mechanism over another. The heterogeneous results for the shear trigger caused the disc edges to curve up by 47 μm (Figs. 9, ,1111--12).12). Experimental measurements that produce a significant change in edge lift and curvature would indicate the presence of heterogeneous growth, possibly ruling out permeation velocity as the dominant trigger because this stimulus is homogeneous throughout the specimen (unless the properties are non-homogeneous on day 0).
The results suggest that a significant increase in PG mass, but not COL mass, will be possible to measure, given the variability in PG and COL contents for this tissue source [37, 48, 49]. PG mass increased markedly (61% and 40% for permeation and shear triggers; GAG/WWI results of Fig. 10) during growth while COL mass increased considerably less (7% and 6% for permeation and shear triggers; COL/WWI results of Fig. 10). To detect such an increase in COL mass will require many specimens. However, these predictions are based on growth laws estimated from studies with bovine AC explants but with different loading protocols [25, 40]; thus, these conclusions are not definitive.
Measuring the geometry of the specimen, as opposed to biochemical and biomechanical measures) may be the best indicator of non-homogeneous growth, as the edge lift change of 47 μm that is predicted to arise in the shear trigger specimen will likely be measurable using optical methods similar to those employed for measuring Poisson's ratios in unconfined compression experiments . Although PG mass varied substantially between the day 0 and day 12 specimens, the PG mass contents among the three shear trigger specimen slices were only slightly different (~7%). Thus, it may be difficult to detect non-homogeneous PG growth in the real experiments. Also, the results suggest that it may be possible to detect significant differences for the day 12 CC modulus and permeability with respect to day 0 but not with respect to the shear trigger slice location (Fig. 8). However, monitoring the trends in non-homogeneous biochemical and biomechanical property evolution may lead to alterations in the experimental design. For example, a power analysis of pilot experiment results along with additional parameter studies with the CGFEM may determine the number of specimens and days of growth needed to detect significant differences arising due to non-homogeneous growth.
Several limitations are related to simplifying assumptions made due to a lack of relevant experimental data. For example, for chondrocytes seeded in agarose culture it has been shown that sulfate and proline incorporation decrease while GAG content increases from ~ 5 mg/ml to ~15 mg/ml during the first 40 days of culture , suggesting that growth rates may be inversely proportional to matrix content. However, the GAG content for our tissue source are > 25 mg/ml (Williamson et al. 2001); consequently, due to the different maturity levels of these studies, it is unclear if the growth rates from studies with tissue engineered constructs can be used here. Since we are not aware of studies with calf AC from our tissue site and depth that have been grown in a permeation bioreactor, our approach is to maintain a relatively simple form in the experimental design stage while planning experiments of different growth durations to quantify the temporal change in the growth functions.
Another limitation regards the assumption of triggers with an “on-off” nature, because the previous studies that support this hypothesis [25, 40] used different loading conditions (i.e. dynamic compression and shear) than the one used here (i.e. steady-state permeation). In order to further test this hypothesis for the shear trigger, experimental results as outlined here can be used to quantify differential growth rates which can then be compared to spatial profiles of shear stress. In order to further test this hypothesis for the permeation trigger, different permeation velocities can be applied; however, this is a concern because higher permeation velocities may trigger cell death via excessive drag-induced matrix compression as discussed above.
Since measurements of geometry changes appear to be a good indicator of heterogeneous growth as discussed above, a limitation of the present study is that the specimens were assumed to be cylindrical in the pre-growth configuration and that curvature develops due to heterogeneous growth. This approach neglects both the possibility of initial curvature and other deviations from a perfectly cylindrical geometry in the pre-growth configuration as well as heterogeneous growth that may result from damage caused by the harvesting procedures. In order to address the first of these two concerns, pilot experiments should include the more accurate representation of the initial geometry via optical methods similar to those employed for measuring Poisson's ratios in unconfined compression experiments  followed by specimen-specific FEA and care should be used when interpreting results in light of possible damage effects.
In an extended simulation of the permeation trigger, the steady state permeation analysis failed to converge on day 14 of the growth process. This phenomenon does not imply a limitation of mixture theory; rather, as discussed above, results as the pore pressure gradients and SM strains increase markedly near the drain side during the growth process. Future efforts should attempt to resolve this limitation.
For the permeation trigger, the predicted increases in GAG content and concentration were 61%, and 37%, respectively (recall that GAG content/concentration is GAG mass normalized by initial/final wet weight, respectively). In several studies with similar specimens cultured under free-swelling conditions [37, 44, 49, 51] GAG content increased markedly while GAG concentration did not increase. This phenomenon has been hypothesized to be partly due to the decrease in collagen-specific pyridinoline crosslink concentration that is consistent among these studies which produces a weaker collagen network that allows substantial tissue swelling [8, 44]. Furthermore, these observed decreases in pyridinoline concentration appears to be related to the observed decreases in both tensile and compressive mechanical properties [8, 44] due to the weakened collagen network. Consequently, post-growth biochemical measurements should include pyridinoline so that, if needed, collagen remodeling can be implemented in later CGFEM analyses. Since it is not known if the collagen network will remodel in the permeation bioreactor, this feature was not implemented in this study.
The majority of compression studies with AC tissue do not report the phenomenon of “stress-softening” which the current model captures: in compression, the modulus first decreases before starting to increase again at some critical strain level. The stress-softening phenomenon may be attributed to the relaxation at small compressive strains of a COL network that is initially in tension to restrain the GAG swelling pressure; it has been detected in several laboratories [21, 44, 52-54]. Any model that attempts to incorporate the “stress-balance” hypothesis of [41, 55] should predict this type of behavior; other than the model used here, similar models have been developed by other authors [52, 54, 56].
In the model used here, it is assumed that all of the PG and COL molecules are “bound” to the extracellular matrix and, consequently, the mass growth functions refer to the net rate of mass deposition into this “bound” state. Since the goal of this study is to study the evolution of continuum-level structure-function properties, the present model does not take into account the permeation of mobile molecules, either newly synthesized or newly degraded, or the fact that the net rate of deposition may be expressed in terms of other factors such as the rate of synthesis of precursor molecules or the rate of conversion from mobile to immobile molecules. However, the general growth mixture theory from which the CGFEM is derived has the flexibility to model these phenomena, as discussed previously [19, 20].
The continued refinement of the CGFEM may lead to a better understanding of the phenomena underlying cartilage growth, degeneration, and repair. For example, difficulties associated with osteochondral graft repair strategies include the development of a smooth convex joint surface, lateral integration with surrounding tissue, concern about transplantation from low to high weight-bearing sites, and mismatch between donor and repair site thickness [5, 57, 58]. In the context of repairing damaged AC, the CGFEM may help to define experimental protocols that may produce a tissue implant with desired non-homogeneous biochemical, biomechanical and geometric properties to help improve integration and survival upon implantation.
This work was supported by the National Science Foundation (Grant No. CMS-0556025) and the National Institutes of Health (NIBIB Grant No. 1R15EB006004-01).
The theory in this section is from [35, 36]. Recall from (12) that general stress constitutive equations are defined in terms of total elastic tensors and relative to distinct PG and COL reference configurations and :
For the PG constituent, an isotropic polyconvex strain energy function is chosen that leads to
where (α1, α2) are material constants determined using an intra- and extra-fibrillar water compartmental model  with representative experimental data  for calf PFG cartilage. For the COL constituent, an anisotropic convex strain energy function is chosen that can be decomposed into isotropic and anisotropic terms:
The isotropic term corresponds to a neo-Hookean material
where μ is the shear modulus, I is the identity tensor, and T is the transpose operator. The anisotropic term corresponds to the bimodular fiber-reinforced strain energy function developed in . Three primary fiber families are orientated along mutually orthogonal basis vectors E1, E2, and E3 where E1 is aligned with the local split-line direction and E3 is normal to the articular surface. Strong interaction terms are generated by six secondary fiber families that lie as symmetric pairs in each of the three planes (1-2, 1-3, 2-3) defined by the basis vectors and at angles of ϕ±12, ϕ±13, and ϕ±23 to E1, E2, and E3 directions: E±12 = cosϕ±12E1 ± sinϕ±12E2, E±13 = cosϕ±13E1 ± sinϕ±13E3, and E±23 ± cosϕ±23E2 ± sinϕ±23E3 . These fiber families generate nine structural tensors
where is the tensor dyadic product. Consequently, the anisotropic term can be expressed as
where λ1, λ2, etc. are the stretches along fiber directions, γ1, γ2, and γ3 are material constants along the primary fibers, δ is the material constant along the secondary fibers, and all material constants are active only when their corresponding fibers are in tension (i.e., λ1 > 1, λ2 > 1, etc.).
Ficklin TP, Davol A, Klisch SM. Simulating the growth of articular cartilage explants in a permeation bioreactor to aid in experimental protocol design. Journal of Biomechanical Engineering, 131:041008, 2009.
In memory of Walter D. Pilkey (1936-2007), Professor of Mechanical Engineering at the University of Virginia, who was the FEA course instructor, MS thesis advisor, and mountain climbing partner of SM Klisch.