Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Cell Mol Bioeng. Author manuscript; available in PMC 2010 September 1.
Published in final edited form as:
Cell Mol Bioeng. 2009 September 1; 2(3): 306–319.
doi:  10.1007/s12195-009-0081-7
PMCID: PMC2782835

Pericellular Matrix Mechanics in the Anulus Fibrosus Predicted by a Three-Dimensional Finite Element Model and In Situ Morphology


Anulus fibrosus (AF) cells have been demonstrated to exhibit dramatic differences in morphology and biologic responses to different types of mechanical stimuli. AF cells may reside as single cell, paired or multiple cells in a contiguous pericellular matrix (PCM), whose structure and properties are expected to have a significant influence on the mechanical stimuli that these cells may experience during physiologic loading of the spine, as well as in tissue degeneration and regeneration. In this study, a computational model was developed to predict the micromechanical stimuli, such as stress and strain, fluid pressure and flow, of cells and their surrounding PCM in the AF tissue using three-dimensional (3D) finite element models based on in situ morphology. 3D solid geometries of cell-PCM regions were registered from serial confocal images obtained from mature rat AF tissues by custom codes. Distinct cell-matrix units were modeled with a custom 3D biphasic finite element code (COMSOL Multiphysics), and simulated to experience uni-axial tensile strain along the local collagen fiber direction. AF cells were predicted to experience higher volumetric strain with a strain amplification ratio (relative to that in the extracellular matrix) of ~ 3.1 – 3.8 at equilibrium, as compared to the PCM domains (1.3 – 1.9). The strain concentrations were generally found at the cell/PCM interface and stress concentration at the PCM/ECM interface. Increased numbers of cells within a contiguous PCM was associated with an apparent increase of strain levels and decreased rate of fluid pressurization in the cell, with magnitudes dependent on the cell size, shape and relative position inside the PCM. These studies provide spatio-temporal information on micromechanics of AF cells in understanding the mechanotransduction in the intervertebral disc.

Keywords: anulus fibrosus, intervertebral disc, cell, pericellular matrix, in situ morphology, mechanics


The intervertebral disc (IVD) is a heterogeneous tissue composed of the central nucleus pulposus and the surrounding anulus fibrosus (AF) that provides load support and flexibility for the spine. The AF tissue is maintained by a sparse population of cells that exhibit characteristics of both chondrocytes and fibroblasts, and have been referred to as “fibrochondrocytic cells”. These AF cells respond to different types of physical stimuli, such as compressive, tensile and shear stresses and strains, hydrostatic and osmotic pressures, and applied electric fields 19,20,48,49,58 with cellular responses that have consequences for altered matrix synthesis, proliferation and other biological phenomena. While in vitro studies are useful to study cellular responses to well-controlled physical stimuli, little is known about the magnitude and sense of mechanical stimuli at the cell-level, within the dense and highly organized extracellular matrix (ECM) of the IVD. Some studies have suggested the cell micromechanical environment to be highly complex, with a mismatch in the magnitude of strain at the tissue and cell level in the anulus fibrosus 10,11. More recent studies of direct strain measurements in other fibrocartilaginous tissues such as the meniscus suggest that strain magnitudes of a cell may be matched to the far-field ECM strains 56. These magnitudes of strain and corresponding stresses and related effects are potentially key regulators of diverse mechanobiologic responses that may occur in disc degeneration and regeneration 48,49.

AF cells are surrounded by a pericellular matrix (PCM) that is rich in type VI collagen 45,46,60 and reside as single cell, paired or multiple cells in a contiguous PCM in the AF tissue 13,29. Approximately half of AF cells in the outer periphery of the IVD have been found to reside in a contiguous PCM shared by two or more cells in the mature rat lumbar disc 13. This PCM is likely to be a key regulator of how mechanical loads are transduced from the ECM to the cell, as seen in articular cartilage 3,23; however, the mechanical consequences of different PCM structures or the presence of multiple cells within a single PCM in the AF remain unknown. Thus, research into the micromechanics of AF cells and their interactions with the surrounding PCM are important in determining the precise micromechanical stimuli experienced by cells that regulate the biological behaviors of AF cells.

Finite element modeling (FEM) has been used to predict cell mechanics where geometries, mechanical loading or mechanical interactions are complex and cannot be easily studied experimentally. Theoretical models of cell-matrix interactions have been developed for single cells in articular cartilage 4,14,24,2628,34,36,38,41,6163, the meniscus 57 and IVD 8,9, incorporating complex features of anisotropy, fluid-solid interactions, inhomogeneity and nonlinearity. In one biphasic model of AF cells developed by Baer and co-workers 8, cells were modeled as two-dimensional elliptical volumes, with results that suggested cell aspect ratio (i.e., elongation) and orientation relative to collagen lamella to be important determinants of the cell micromechanical environment. A previous study of three-dimensional (3D) in situ morphology of the PCM and cells in the rat IVD has shown that cell-matrix units (CMU, i.e., cell(s) and their associated PCM) in the outermost AF were more elongated and discoidal in shape than elliptical, however, with aspect ratios greater than 1 13. These findings indicate that PCM regions and cells in the AF possess fully 3D geometries that may not be well approximated using axi-symmetric spheroids. Furthermore, prior models of the AF cell mechanics have not incorporated a PCM region, which has since been demonstrated to significantly affect predicted cell stress and strain, and fluid pressure and flow profiles for chondrocytes 4,22,28,34. These observations suggest that a fully 3D FEM incorporated with 3D in situ morphology of the PCM and enclosed cell(s), and the material parameters of the PCM, as well as material anisotropy of the ECM, will be required to accurately assess the micromechanics of AF cells.

The objective of this study was to predict region-specific micromechanical stimuli by a theoretical model based on 3D in situ morphology of the PCM and enclosed cell(s) in the AF. The 3D solid geometries of the PCM and cell regions were constructed from serial confocal images of CMUs from mature rat AF tissues using custom codes (MATLAB and COMSOL Script). Three-dimensional finite element models of the extracellular matrix, PCM and cell sub-domains under idealized loading configurations were analyzed by a custom FEM code (COMSOL Multiphysics) that incorporates linearly elastic, biphasic behaviors and material anisotropy for the AF extracellular matrix. The predicted micromechanical environment, (e.g., stress, strain, fluid pressure and flow) experienced by inner or outer AF cells residing as single cells, paired, or multiple cells within a contiguous PCM were compared to reveal any regional or subgroup difference.


2.1 3D biphasic finite element method (FEM)

Different finite element formulations of the biphasic mixture theory 42 have been applied to simulate soft tissue mechanics 5,40,51,52,59. A 3D multi-phasic formulation of biphasic and triphasic mixture models 37,42 has recently been implemented in COMSOL Multiphysics (COMSOL, Inc., Burlington, MA) to analyze the mechanical, chemical and electrical signals and solute transport in the disc tissue under compression 30,64,65. Similarly, a 2D axi-symmetric u-v-p formulation of biphasic theory has also been implemented in COMSOL to model the dynamic responses of chondrons in cartilage 34. These models based on a u-p or u-v-p FEM formulation predict behaviors that are well-matched to those reported for earlier formulations of biphasic mixtures 5,40,51,52,59. In this study, a 3D mixed u-p formulation of the linear biphasic theory 40,59 was implemented in COMSOL, with the following coupled weak form in a given domain Ω with a boundary Γ


where u is the displacement vector (with components u,v,w); p is the pressure for the inviscid and incompressible fluid phase; w and q are admissible trial functions in the weighting space; k is a constant value for the hydraulic permeability; and t is the time. The effective Cauchy stress tensor for the linearly elastic solid phase s, is given by σEs=Cεs, where C is a fourth-order tensor of elastic coefficients; εs = ([big down triangle, open]u + [big down triangle, open]uT/ 2 is the infinitesimal strain tensor; t=σEsn is the prescribed surface traction on the boundary with an outward normal n; and Q=kpn is the fluid flux on the boundary.

The weak form formulations were implemented in COMSOL using “custom weak form” application mode with the structural mechanics module. The implemented 3D biphasic code was validated against the analytical solutions of viscoelastic responses of soft tissues in multiple experimental configurations including static or dynamic loading in confined or unconfined compression 6,42,53.

2.2 3D geometry registration and meshing

3D reconstructions and morphometric measurements of the PCM and cell regions were obtained in situ using fluorescence confocal microscopy of en bloc sections of the AF from the mature rat lumbar disc immunolabeled for type VI collagen 13. A total of ten distinct cell-matrix units (CMU) were chosen for the FEM in the inner (n=6) and outer AF (n=4) with aspect ratios around the mean values in each CMU subgroup (1 cell, 2 or 3+ cells). The reconstructed surfaces of the PCM and enclosed cell(s) in the CMU were registered as 3D solid geometry objects in COMSOL, respectively, by custom programs written in MATLAB (The Mathworks, Inc., Natick, MA) and COMSOL Script (COMSOL, Inc., Burlington, MA).

In detail, the reconstructed surface of the PCM in the original Cartesian coordinates was re-oriented through translation and rotation to be aligned with the principal axis of the CMU as the z-axis for FEM development, and with the center of the CMU as the origin of the new coordinates. The same transformation was operated on the reconstructed cell surface within the PCM for exact cell positioning. The transformed surface was re-sectioned to 20 ~ 70 z-slices per CMU or cell at a gradient of intervals with thicker slices in the middle of the CMU or cell (up to 2 μm) and thin slices on the tip of the CMU or cell (as low as 0.25 μm). The surface contour on each slice was fitted into an ellipse with 4 points defining its long and short axes. All fitted 2D surface contours (ellipses) were then stacked into a 3D solid object using a COMSOL Script function “LOFT” (Figure 1). The 3D solid objects of the PCM and cell(s), therefore, were imported into COMSOL and a cube representing the ECM with dimensions (a × a × h) of at least three times of the PCM dimensions was added to the geometry, in which a is the width of the transverse plane of the ECM, and h is the height of the ECM along the longitudinal direction (the principal axis of the CMU). A preliminary study confirmed that the choice of this ECM size was appropriate to predict cellular viscoelastic responses independent of geometry of the PCM and ECM sub-domains. Thus, a 3D solid geometry model with three sub-domains including the ECM, PCM and enclosed cell(s), was built based on their in situ 3D morphologies. These three sub-domains were exclusive to each other, but physically connected to each other (continuous displacement and pressure).

Figure 1
Geometry registration and meshing of the cell-matrix interaction model for a representative cell-matrix unit (CMU). (a) 3D reconstructed CMU; (b) 3D registered geometry of the CMU in tetrahedral element meshes, viewed in YZ plane; (c) Elliptical surface ...

The combined 3D solid objects were meshed using tetrahedral elements in COMSOL Multiphysics, with a linear shape function for pressure and quadratic shape function for displacement (Figure 1). The registered 3D geometries required a large number of tetrahedral elements to generate a good quality of meshes due to the curvature and roughness on the surfaces of the PCM and cells, as seen in examples of meshed CMUs from inner and outer AF regions (Figure 2). The meshes generally consisted of 18,000 ~ 35,000 tetrahedral elements and approximately 75,000 ~ 150,000 degrees of freedom, dependent on the CMU size and shape. The average minimum element quality index across all models was ~ 0.15, satisfying the requirement for tetrahedral elements in COMSOL (> 0.1).

Figure 2
Registered 3D solid geometries of the pericellular matrix and cell(s) in the annulus fibrosus (AF) in tetrahedral element meshes. Examples shown are from the inner AF (1, 3, 5) and outer AF (2, 4, 6), including 1 cell (1&2), 2 cells (3&4), ...

2.3 Material properties

The AF exhibits greater stiffness in the longitude,nal direction of fiber bundles than in directions perpendicular to the fibers (transverse plane) 1,16,17,50, necessitating the incorporation of material anisotropy for the ECM. Therefore, a transversely isotropic material was used to model the solid matrix of the AF with a dependence on five independent material constants: two Young's moduli in the longitudinal (EL) and transverse directions (ET) with respect to the local collagen fiber, Poisson's ratios νT and νLT, and one shear modulus (GLT)5. The solid phases of the PCM and cells were assumed to be isotropic with a dependence on only two independent material constants (Young's modulus E and Poisson's ratio ν). All three sub-domains were assumed to be linearly elastic, biphasic materials with a constant permeability.

The material constants were chosen from the literature for the ECM 17,21,32,35,43,55, PCM 2,4 and cells 25, as seen in Table 1.

Table 1
Material properties chosen for the extracellular matrix, pericellular matrix and cell in the annulus fibrosus

2.4 Boundary conditions

The finite element models were used to simulate a stress relaxation response under an idealized tensile deformation configuration along the longitudinal direction, for possible loading conditions 54 that correspond to: (Case 1) uni-axial tensile displacement instantaneously applied at the longitudinal surface (z = h / 2, step displacement corresponding to 5% strain across all models incorporated with a smoothed Heaviside function) of a 3D ECM volume along the longitudinal direction with fixed boundaries on all transverse surfaces; (Case 2) uni-axial tensile displacement instantaneously applied at longitudinal surface of a 3D ECM volume along the longitudinal direction with traction-free boundaries on all transverse surfaces. The displacements at the bottom surface (z = −h / 2) were fixed in all directions. It was assumed to be free draining (p = 0) on all surfaces, and the interfaces between different sub-domains were assumed to be fully bonded.

2.5 Solution procedure and post-processing

The models were solved using a generalized minimal residual method (GMRES) as an iterative linear system solver and a backward Euler difference approximation for time-dependent analysis. Simulations were performed on a Dell Inspiron, 3.2 GHz 32-bit computer with 2 GB RAM. The total simulated time was set as 20 seconds for all AF samples, which had been shown to be sufficient for different sub-domains to reach equilibrium for both cases in the preliminary study. The convergence was examined by refining the mesh and tightening the tolerance. The computational cost was 10,000 ~ 15,000 seconds across all models. Post processing results for stress, strain, fluid pressure and gradient were obtained using COMSOL post-processing features. Additional measures such as the average volumetric strain over the entire sub-domain, εν, and average fluid pressure, p, were calculated from the integrals of pressure and strain in each sub-domain, Ω, according to the following equations, respectively.


where V represented the volume of the sub-domain. Note that these averaging measurements provide valuable information of the overall response for each sub-domain with complex 3D geometry; however, their values in the ECM domain may be underestimated due to a boundary effect caused by averaging over all elements in the ECM that include the boundary elements on which the pressure or displacement was set to zero.

An effective von Mises stress at each spatial position in the sub-domain, σe, was defined as


as well as a deviatoric strain invariant, εdev,


where εm =(εxxyyzz/ 3. These measures were compared across models to detect any difference in kinematic and mechanical variables predicted to occur at ECM, PCM or cellular scales, as well as to detect differences between the inner and outer AF regions, amongst three CMU subgroups, and between the two computational cases (Case 1 and Case 2).


Case 1

The double inclusions of the PCM and enclosed cell(s) embedded in the ECM generally did not affect predictions of tissue mechanics in the majority of the ECM domain. Stress, strain, and fluid pressure were uniform in the relatively far-field of the ECM and inside the cell, as compared to highly non-uniform field effects observed in the PCM (Figures 3, ,5,5, ,66 and and8).8). At equilibrium, the average volumetric strain in the far-field ECM strain was εν=0.05; as expected, strain amplification effects over this value were observed in the PCM and cell domains (Figures 3 and and4)4) largely due to the more compliant material properties assigned to both the PCM and cell domains. The deviatoric strains in the ECM in the inner (data not shown) and outer AF (Figure 5) was uniform with a lower average value ~ 0.03. Finally, the effective von Mises stress was generally high in the ECM (~ 20 kPa) for all times.

Figure 3
Equilibrium volumetric strains at multiple length scales shown for loading conditions of Case 1 and Case 2. Strain distribution was generally uniform in the extracellular matrix (ECM) and cells, but non-uniform in the pericellular matrix (PCM). Strain ...
Figure 4
Temporal response of average volumetric strain at multiple length scales in Case 1. Comparisons between the pericellular matrix (PCM, a&b) and cell (c&d) domains in the 1 cell, 2 or 3+ cell cell-matrix unit subgroups showed that the average ...
Figure 5
Temporal responses of deviatoric strain at multiple length scales in the outer annulus fibrosus (AF) for loading described in Case 1. A gradient of deviatoric strain was shown at early times with the highest value observed in the cell domain. Strain concentrations ...
Figure 6
Temporal responses of the effective von Mises stress in the pericellular matrix of the inner annulus fibrosus for loading described in Case 1. The stress was highly heterogeneous in the pericellular matrix and relaxed over the time. Stress concentrations ...
Figure 8
Temporal responses of the pressure gradient predicted for the outer annulus fibrosus for loading conditions of Case 1. A pressure gradient, indicating an overall inward fluid flow after tensile deformation existed in all sub-domains at early times. The ...

Upon application of the step displacement, there was a transient cell volume loss predicted at early times as measured by values for εν (Figure 4). This was followed by a volume gain at later times in the cell domain due to the asymmetric cell positioning in the PCM and a net effect of spatially varying pressure gradients that drove fluid flow into or out of the cell at a local scale. The average volumetric strain at equilibrium in the PCM domain was modestly amplified over values for the ECM in a relatively narrow spatial band within the PCM domain, giving rise to a PCM strain amplification ratio of ~ 1.8. Much higher values for the strain amplification ratio were predicted for the cell domain with a dependence of volumetric strain on the cell number enclosed in one PCM (1 cell ~ 3.0 and 3+ cells ~ 3.8).

Average deviatoric strains were modestly higher in PCM (~ 0.04) and cell domains (~ 0.06) compared to ECM domains at early times (< 3s) (Figure 5), with predictions that the differences amongst domains decreased over time. The deviatoric strain concentration, however, was highest at the cell/PCM interface (up to ~ 0.15 in the inner AF and ~ 0.10 in the outer AF) with a distribution that depended on the shape and surface curvature of the individual cell surface (Figure 5). It is noteworthy that this strain concentration did not occur at, nor have an obvious dependence on geometry of the PCM/ECM interface. In contrast, the effective von Mises stress was highly heterogeneous in the PCM with stress concentrations that were seen mainly at the PCM/ECM interface (peak values of ~ 7 kPa in the inner AF (Figure 6) and ~ 8 kPa in the outer AF (data not shown)) at short times after loading.

The transient rise and decay of fluid pressure showed similar trends at ECM and PCM (Figure 7), with a fluid pressurization effect that occurred more rapidly in the outer AF across all models (~ 5s), as compared to the inner AF (~ 10s). This effect is likely due to a Young's modulus in the ECM of the outer AF that is nearly double that of the inner AF. The time to peak fluid pressure decreased with an increasing number of cells enclosed in one PCM, as compared to the 1 cell CMU subgroup. The fluid pressure gradient ([big down triangle, open]p ), that drives the fluid flow, was non-uniform in the narrow region of the PCM within a range of 107 – 109 N/m3 at early times after loading, consistent with the predictions that cell and PCM volumetric losses and gains occur over the short term (Figure 8).

Figure 7
Temporal responses of average fluid pressure in the pericellular matrix (PCM) and cell for loading conditions in Case 1. The average fluid pressures in the PCM (a&b) and cell (c&d) domains in the 1 cell, 2 or 3+ cell cell-matrix unit subgroups ...

Case 2

The mechanical predictions of the cell-PCM-ECM finite element models for the AF under uni-axial tensile displacements with traction-free transverse surfaces exhibited trends for stress, strain, fluid pressure and flow that were similar to that described for Case 1, except that the magnitudes of these measures were on different scales. The average volumetric strain at equilibrium in the ECM domain of the inner AF was reduced over that reported for Case 1 (~ 0.03) as expected, and corresponding values for average volumetric strain in the PCM domain were slightly higher (~ 0.04) with a strain amplification ratio of ~ 1.3. At the cellular level, the average volumetric strain was also higher in the 1 cell-containing CMU (~ 0.035) and the 3+ cell CMU (~ 0.08) with strain amplification ratios for the cells ranged from 1.2 to 2.7, lower than those in Case 1 (Figure 3). The sites for volumetric and deviatoric strain concentration were at the cell/PCM interface, as for Case 1, although significant spatial variations in volumetric strain distribution were observed for the laterally unconstrained case (Figure 3); positive values (volume gain) for volumetric strain were predicted mainly on the transverse faces of the cell/PCM interface and negative values (volume loss) were predicted for the longitudinal faces (Figure 3). The deviatoric strain and von Mises stress values and distributions were similar to that reported in Case 1 in PCM and cell domains (data not shown).

The time to peak fluid pressure averaged throughout the PCM or cell domain for laterally unconstrained loading was significantly longer than for Case 1 (Figure 10). In general, the time to peak pressure decreased with the increasing number of cells enclosed in one PCM, as compared to the 1 cell CMU subgroups, similarly as for Case 1. The pressure gradients also varied spatially in the PCM and cell domains (data not shown), depending on the local PCM geometry and cell position.


The mechanical interactions of cells, their PCM and ECM were predicted for AF cells using 3D, linearly elastic, anisotropic and biphasic finite element models based on their in situ morphologies for the first time. Model predictions for idealized loading configurations suggest that significant differences exist in mechanical factors at ECM, PCM, and cellular scales in both inner and outer AF regions, with a dependence on in situ geometry and assumed relative mechanical properties of the ECM, PCM and cells. In general, AF cells were predicted to experience volumetric strains that were significantly larger than those in the PCM and corresponding far-field ECM, with strain amplification ratios that varied from 1.2 up to 3.8, depending on the morphology, AF domain and loading conditions. The observation of strain amplification was expected for the PCM and cell domains, as the Young's modulus of the ECM domain was assumed to be 4 – 10 times higher than that of the PCM, and 1000 – 2000 times higher than that of the cell domain. The predicted values for volumetric strain were similar to that measured for chondrocytes within the PCM of articular cartilage 15.

The rate of fluid pressurization at the cell scale was generally slower than that in the PCM and ECM, even with the highest value for the assumed permeability. Predictions of the fluid flow and fluid pressurization will be particularly dependent on the selected material properties, particularly those related to permeability coefficients for the three different domains. Here, the assumed values for Poisson's ratio and permeability between the PCM and cell domains were associated with temporally varying pressure gradient predictions that suggested a pattern of short-term outward fluid flow from the cell, followed by a period of inward-directed fluid flow. While these flow patterns and pressurization profiles have potential to regulate cell signaling and solute transport, their predictions here are dependent on this choice of material properties and additional work would be required to validate this choice. Nevertheless, the pressure gradients were predicted to vary spatially in the narrow PCM region, resulting from its irregular in situ geometry and asymmetric cell positioning.

While the magnitudes for deviatoric strain in the PCM and cell were also amplified compared to ECM scales, this amplification effect was lower than for volumetric strains except at the cell/PCM interface where strain concentrations were observed. These findings for significant magnitudes of deviatoric strain at the cell/PCM interface point towards the potential for altered cell surface areas upon loading. In contrast, von Mises stress values were always predicted to be highest in the ECM and lowest in the cell, principally due to the spatially varying mechanical properties associated with these domains. Together, these model predictions reinforce the hypothesis that the PCM may play an important role in shielding the cells from instant fluid pressurization while enabling the transmission of kinematic variables of deviatoric and volumetric strain 24.

A major motivation for the current study was to evaluate the potential mechanical influence of multiple cells sharing one contiguous PCM within the AF, based on physically realistic cell and PCM morphologies. In general, increasing numbers of cells within one PCM was associated with higher volumetric and deviatoric strains within the cells and a slower time for pressure decay (i.e., sustained pressurization effect) under equivalent loading conditions. As an example, predicted average volumetric strain amplification ratio at equilibrium in the PCM of AF tissues ranged from 1.3 to 1.9 while this ratio at the cellular scale depended on the number of cells enclosed in one PCM; values were highest in cells residing in 3+ cell CMUs (up to 3.8) and lowest in 1 cell CMUs (~ 3). This prediction indicates that cells of physiologically realistic geometries and similar mechanical properties are mechanically interacting within the CMU with behaviors that differ from that of single cell-PCM units. In previous studies, investigators have reported a shift towards CMUs containing fewer cells with increased age in the AF 13. Taken together with the model predictions presented here, these findings suggest that the micromechanical stimuli experienced by AF cells may change with age, and thus, potentially contribute to the observed cell senescence with aging of the IVD.

Regional differences of cell-matrix interactions were predicted by region-specific models in the inner and outer AF, such as larger volumetric strains, higher effective PCM stresses, and lower peak deviatoric strains for cells in the outer AF as compared to inner AF regions. Many of these effects likely relate to the higher fluid pressurization rate in the outer AF region, which was almost double than that in the inner AF. These differences can arise from the adopted region-specific geometries (i.e., higher aspect ratio of PCM for cells in the outer AF) and region-specific material properties (higher anisotropic effects in the ECM in the outer AF, i.e., higher ratio of longitudinal and transverse stiffness) and were indeed, not surprising findings. The trend for these model predictions is in agreement with previous models of AF cell mechanics, which demonstrate that cell aspect ratio and ECM anisotropy strongly impacted the magnitudes of predicted stresses and strains in the microenvironment of AF cells 8.

A significant advantage of the FEM predictions for cell-matrix interactions presented here is that all 3D geometries of the PCM and cells were registered from their in situ morphologies, as compared to previous 2D axi-symmetric models of single cell mechanics in the AF 8,9. The use of directly-measured 3D geometries enabled the prediction of spatially heterogeneous strains and stresses near the interfaces of cell/PCM and PCM/ECM, and indicated the presence of significant magnitudes of deviatoric stresses and strains larger than that predicted for the axi-symmetric models. Moreover, the PCM containing paired or multiple cells with realistic spatial arrangement of individual cells within the matrix has furthered an understanding of potential effects of cells sharing a contiguous PCM on cell-matrix mechanics. A major effect not incorporated in these passive models of cell-matrix mechanics would be any electrical or chemical interactions amongst cells within one PCM, and between cell, PCM and ECM that were not physically modeled here. Such interactions have been successfully represented by triphasic models of cell-matrix interactions in studies of single chondrocytes and IVD cells 30,47,64,65 that suggest flow-related chemical gradients and electric field effects could regulate the kinematic variables of strain and also measures of stress as predicted here. These phenomena can be expected to have very significant effects on mechanotransduction of IVD cells through the generation of gradients of chemical or electrical potentials across the cell membrane. Incorporation of physicochemical effects into a COMSOL algorithm for cell-matrix modeling has already been achieved 64,65 and could be extended to study the in situ models of cell-matrix units in the AF as performed here.

Studies of the in situ strain fields within the IVD and other collagenous tissues suggest a complex mechanical state may exist at the length scale of a cell under even simplified loading configurations 7,10,11,56. In this study, simplified loading conditions were applied in the near cell region to uncouple effects of complex ECM loading from the complex effects generated by models of in situ cell matrix unit morphology obtained from normal rat disc cells. More accurate boundary conditions, representing kinematics in local normal and pathological disc tissues, could be implemented and subsequently transferred to the microscale. Challenges persist in understanding the local boundary conditions needed for application to the ECM in these cell-matrix models, however, as studies of collagenous tissues suggest that identifying cell-PCM adhesion sites is very challenging with effects such as discontinuities between cells and their PCM may exist at the local scale 7,10,11. In a recent experimental study of local strains within the fibrocartilaginous meniscus tissues, a nearly 1:1 coupling of ECM strain to local strain at the cellular length scales was measured for principal strain values, but not axial and transverse strain components 56. These observations suggest the existence of out-of-plane displacements and rotations upon loading of the ECM, even in simplified geometries. These observations indicate that the interface conditions between the ECM and PCM, and potentially the PCM and cell, may be significantly more complex than assumed in the present model. Alternate mechanical models for these interfaces, which incorporate discrete point contacts (i.e. adhesion sites between the cell and matrix) or active behavior such as time-varying cell-matrix bonding force, may lead to dramatically different model predictions for the micromechanical environment of the cell. While this model has addressed the transmission of mechanical stimuli from ECM to PCM and overall cellular scales, this model can be improved by adding more elements representing some subcellular structures such as a glycocalyx, cell membrane or cortex, cytoskeleton and nucleus, to further investigate mechanotransduction through these substructures.

In summary, cell-matrix mechanical interactions were predicted for AF cells using 3D, anisotropic, linearly elastic, and biphasic finite element models based on in situ morphology of the PCM and cells. Model predictions of region-specific micromechanics of the PCM and cells from the inner and outer AF regions were associated with differences in the geometries of the PCM and cells, zone of cell origin in the AF, spatial arrangement of multiple cells in one PCM, as well as material properties at multiple length scales. The incorporation of PCM regions in the FEM highly alters predictions for the local micromechanics of AF cells, as compared to previous models of AF cell mechanics 8,9, and suggests that the PCM may act to attenuate fluid pressurization and amplify strain in the cell as compared to single cell mechanical predictions3,22,34. These predicted micromechanical stimuli in the vicinity of AF cells may be sensed and transduced to ultra-structures inside the cells by different mechanisms, such as mechano- or osmo-sensitive ion channels on the cell membrane 44 or cell processes that may extend from the AF cell membrane to other cells or surrounding matrix 8,12,18,33. These physical signals and associated cellular responses are known to be important determinants of cell metabolic behaviors that govern cell survival and differentiation, matrix synthesis and degradation 31,39,48,49. The spatio-temporal responses of micromechanical stimuli in the vicinity of AF cells predicted by the geometrically accurate FEM models developed here provide useful information in understanding the magnitudes of anticipated micromechanical factors expected to affect mechanotransduction for cells of the IVD.

Figure 9
Average fluid pressure in the pericellular matrix (PCM) and cell for the loading conditions of Case 2. A complex pressurization pattern was observed at the pericellular (a&b) and cellular (c&d) scales for both inner (a&c) and outer ...


This study was supported by NIH grants AR47442, AR50245, and AG15768.


1. Acaroglu ER, Iatridis JC, Setton LA, et al. Degeneration and aging affect the tensile behavior of human lumbar anulus fibrosus. Spine. 1995;20(24):2690–2701. [PubMed]
2. Alexopoulos LG, Haider MA, Vail TP, et al. Alterations in the mechanical properties of the human chondrocyte pericellular matrix with osteoarthritis. Journal of Biomechanical Engineering. 2003;125(3):323–333. [PubMed]
3. Alexopoulos LG, Setton LA, Guilak F. The biomechanical role of the chondrocyte pericellular matrix in articular cartilage. Acta Biomater. 2005;1(3):317–325. [PubMed]
4. Alexopoulos LG, Williams GM, Upton ML, et al. Osteoarthritic changes in the biphasic mechanical properties of the chondrocyte pericellular matrix in articular cartilage. Journal of Biomechanics. 2005;38(3):509–517. [PubMed]
5. Almeida ES, Spilker RL. Finite element formulations for hyperelastic transversely isotropic biphasic soft tissues. Computer Methods in Applied Mechanics and Engineering. 1998;151(3–4):513–538.
6. Armstrong CG, Lai WM, Mow VC. An Analysis of the Unconfined Compression of Articular-Cartilage. Journal of Biomechanical Engineering-Transactions of the Asme. 1984;106(2):165–173. [PubMed]
7. Arnoczky SP, Lavagnino M, Whallon JH, et al. In situ cell nucleus deformation in tendons under tensile load; a morphological analysis using confocal laser microscopy. J Orthop Res. 2002;20(1):29–35. [PubMed]
8. Baer AE, Laursen TA, Guilak F, et al. The micromechanical environment of intervertebral disc cells determined by a finite deformation, anisotropic, and biphasic finite element model. J Biomech Eng. 2003;125(1):1–11. [PubMed]
9. Baer AE, Setton LA. The micromechanical environment of intervertebral disc cells: effect of matrix anisotropy and cell geometry predicted by a linear model. J Biomech Eng. 2000;122(3):245–251. [PubMed]
10. Bruehlmann SB, Hulme PA, Duncan NA. In situ intercellular mechanics of the bovine outer annulus fibrosus subjected to biaxial strains. J Biomech. 2004;37(2):223–231. [PubMed]
11. Bruehlmann SB, Matyas JR, Duncan NA. ISSLS prize winner: Collagen fibril sliding governs cell mechanics in the anulus fibrosus: an in situ confocal microscopy study of bovine discs. Spine. 2004;29(23):2612–2620. [PubMed]
12. Bruehlmann SB, Rattner JB, Matyas JR, et al. Regional variations in the cellular matrix of the annulus fibrosus of the intervertebral disc. J Anat. 2002;201(2):159–171. [PubMed]
13. Cao L, Guilak F, Setton LA. Three-dimensional morphology of the pericellular matrix of intervertebral disc cells in the rat. J Anat. 2007;211(4):444–452. [PubMed]
14. Chahine NO, Hung CT, Ateshian GA. In-situ measurements of chondrocyte deformation under transient loading. Eur Cell Mater. 2007;13:100–111. discussion 111. [PubMed]
15. Choi JB, Youn I, Cao L, et al. Zonal changes in the three-dimensional morphology of the chondron under compression: The relationship among cellular, pericellular, and extracellular deformation in articular cartilage. J Biomech. 2007;40:2596–2603. [PMC free article] [PubMed]
16. Ebara S, Iatridis JC, Setton LA, et al. Tensile properties of nondegenerate human lumbar anulus fibrosus. Spine. 1996;21(4):452–461. [PubMed]
17. Elliott DM, Setton LA. Anisotropic and inhomogeneous tensile behavior of the human anulus fibrosus: experimental measurement and material model predictions. J Biomech Eng. 2001;123(3):256–263. [PubMed]
18. Errington RJ, Puustjarvi K, White IR, et al. Characterisation of cytoplasm-filled processes in cells of the intervertebral disc. J Anat. 1998;192(Pt 3):369–378. [PubMed]
19. Ferguson SJ, Ito K, Nolte LP. Fluid flow and convective transport of solutes within the intervertebral disc. J Biomech. 2004;37(2):213–221. [PubMed]
20. Gu WY, Justiz MA, Yao H. Electrical conductivity of lumbar anulus fibrosis: effects of porosity and fixed charge density. Spine. 2002;27(21):2390–2395. [PubMed]
21. Gu WY, Mao XG, Foster RJ, et al. The anisotropic hydraulic permeability of human lumbar anulus fibrosus. Influence of age, degeneration, direction, and water content. Spine. 1999;24(23):2449–2455. [PubMed]
22. Guilak F. The deformation behavior and viscoelastic properties of chondrocytes in articular cartilage. Biorheology. 2000;37(1–2):27–44. [PubMed]
23. Guilak F, Alexopoulos LG, Upton ML, et al. The pericellular matrix as a transducer of biomechanical and biochemical signals in articular cartilage. Ann N Y Acad Sci. 2006;1068:498–512. [PubMed]
24. Guilak F, Mow VC. The mechanical environment of the chondrocyte: a biphasic finite element model of cell-matrix interactions in articular cartilage. J Biomech. 2000;33(12):1663–1673. [PubMed]
25. Guilak F, Ting-Beall HP, Baer AE, et al. Viscoelastic properties of intervertebral disc cells. Identification of two biomechanically distinct cell populations. Spine. 1999;24(23):2475–2483. [PubMed]
26. Haider MA. A radial biphasic model for local cell-matrix mechanics in articular cartilage. SIAM J. Appl. Math. 2004;64(5):1588–1608.
27. Haider MA, Guilak F. Application of a three-dimensional poroelastic BEM to modeling the biphasic mechanics of cell-matix interactions in articular cartilage, Comput. Methods Appl. Mech. Engrg. 2007;196:2999–3010. [PMC free article] [PubMed]
28. Haider MA, Schugart RC, Setton LA, et al. A mechano-chemical model for the passive swelling response of an isolated chondron under osmotic loading. Biomech Model Mechanobiol. 2006;5(2–3):160–171. [PubMed]
29. Hastreiter D, Ozuna RM, Spector M. Regional variations in certain cellular characteristics in human lumbar intervertebral discs, including the presence of alpha-smooth muscle actin. J Orthop Res. 2001;19(4):597–604. [PubMed]
30. Iatridis JC, Laible JP, Krag MH. Influence of fixed charge density magnitude and distribution on the intervertebral disc: applications of a poroelastic and chemical electric (PEACE) model. J Biomech Eng. 2003;125(1):12–24. [PubMed]
31. Iatridis JC, MacLean JJ, Roughley PJ, et al. Effects of mechanical loading on intervertebral disc metabolism in vivo. J Bone Joint Surg Am. 2006;88(Suppl 2):41–46. [PubMed]
32. Iatridis JC, Setton LA, Foster RJ, et al. Degeneration affects the anisotropic and nonlinear behaviors of human anulus fibrosus in compression. J Biomech. 1998;31(6):535–544. [PubMed]
33. Johnson WE, Roberts S. Human intervertebral disc cell morphology and cytoskeletal composition: a preliminary study of regional variations in health and disease. J Anat. 2003;203(6):605–612. [PubMed]
34. Kim E, Guilak F, Haider MA. The dynamic mechanical environment of the chondrocyte: a biphasic finite element model of cell-matrix interactions under cyclic compressive loading. J Biomech Eng. 2008;130(6):061009. [PMC free article] [PubMed]
35. Klisch SM, Lotz JC. A special theory of biphasic mixtures and experimental results for human annulus fibrosus tested in confined compression. J Biomech Eng. 2000;122(2):180–188. [PubMed]
36. Korhonen RK, Julkunen P, Rieppo J, et al. Collagen network of articular cartilage modulates fluid flow and mechanical stresses in chondrocyte. Biomech Model Mechanobiol. 2006;5(2–3):150–159. [PubMed]
37. Lai WM, Hou JS, Mow VC. A triphasic theory for the swelling and deformation behaviors of articular cartilage. J Biomech Eng. 1991;113(3):245–258. [PubMed]
38. Likhitpanichkul M, Guo XE, Mow VC. The effect of matrix tension-compression nonlinearity and fixed negative charges on chondrocyte responses in cartilage. Mol Cell Biomech. 2005;2(4):191–204. [PubMed]
39. Lotz JC, Hsieh AH, Walsh AL, et al. Mechanobiology of the intervertebral disc. Biochem Soc Trans. 2002;30(Pt 6):853–858. [PubMed]
40. Meng XN, LeRoux MA, Laursen TA, et al. A nonlinear finite element formulation for axisymmetric torsion of biphasic materials. International Journal of Solids and Structures. 2002;39(4):879–895.
41. Michalek AJ, Iatridis JC. A numerical study to determine pericellular matrix modulus and evaluate its effects on the micromechanical environment of chondrocytes. J Biomech. 2007;40(6):1405–1409. [PubMed]
42. Mow VC, Kuei SC, Lai WM, et al. Biphasic creep and stress-relaxation of articular-cartilage in compression - theory and experiments. Journal of Biomechanical Engineering-Transactions of the Asme. 1980;102(1):73–84. [PubMed]
43. Perie D, Korda D, Iatridis JC. Confined compression experiments on bovine nucleus pulposus and annulus fibrosus: sensitivity of the experiment in the determination of compressive modulus and hydraulic permeability. J Biomech. 2005;38(11):2164–2171. [PubMed]
44. Pritchard S, Guilak F. The role of F-actin in hypo-osmotically induced cell volume change and calcium signaling in anulus fibrosus cells. Ann Biomed Eng. 2004;32(1):103–111. [PubMed]
45. Roberts S, Ayad S, Menage PJ. Immunolocalisation of type VI collagen in the intervertebral disc. Ann Rheum Dis. 1991;50(11):787–791. [PMC free article] [PubMed]
46. Roberts S, Menage J, Duance V, et al. 1991 Volvo Award in basic sciences. Collagen types around the cells of the intervertebral disc and cartilage end plate: an immunolocalization study. Spine. 1991;16(9):1030–1038. [PubMed]
47. Sengers BG, Oomens CW, Baaijens FP. An integrated finite-element approach to mechanics, transport and biosynthesis in tissue engineering. J Biomech Eng. 2004;126(1):82–91. [PubMed]
48. Setton LA, Chen J. Cell mechanics and mechanobiology in the intervertebral disc. Spine. 2004;29(23):2710–2723. [PubMed]
49. Setton LA, Chen J. Mechanobiology of the intervertebral disc and relevance to disc degeneration. J Bone Joint Surg Am. 2006;88(Suppl 2):52–57. [PubMed]
50. Skaggs DL, Weidenbaum M, Iatridis JC, et al. Regional variation in tensile properties and biochemical composition of the human lumbar anulus fibrosus. Spine. 1994;19(12):1310–1319. [PubMed]
51. Spilker RL, Maxian TA. A Mixed-Penalty Finite-Element Formulation of the Linear Biphasic Theory for Soft-Tissues. International Journal for Numerical Methods in Engineering. 1990;30(5):1063–1082.
52. Spilker RL, Suh JK, Mow VC. Effects of Friction on the Unconfined Compressive Response of Articular-Cartilage - a Finite-Element Analysis. Journal of Biomechanical Engineering-Transactions of the Asme. 1990;112(2):138–146. [PubMed]
53. Suh JK, Li ZF, Woo SLY. Dynamic Behavior of a Biphasic Cartilage Model under Cyclic Compressive Loading. Journal of Biomechanics. 1995;28(4):357–364. [PubMed]
54. Tsantrizos A, Ito K, Aebi M, et al. Internal strains in healthy and degenerated lumbar intervertebral discs. Spine. 2005;30(19):2129–2137. [PubMed]
55. Umehara S, Tadano S, Abumi K, et al. Effects of degeneration on the elastic modulus distribution in the lumbar intervertebral disc. Spine. 1996;21(7):811–819. discussion 820. [PubMed]
56. Upton ML, Gilchrist CL, Guilak F, et al. Transfer of macro-scale tissue strain to micro-scale cell regions in the deformed meniscus. Biophys J. 2008 [PubMed]
57. Upton ML, Guilak F, Laursen TA, et al. Finite element modeling predictions of region-specific cell-matrix mechanics in the meniscus. Biomech Model Mechanobiol. 2006;5(2–3):140–149. [PubMed]
58. Urban JP. The role of the physicochemical environment in determining disc cell behaviour. Biochem Soc Trans. 2002;30(Pt 6):858–864. [PubMed]
59. Wayne JS, Woo SLY, Kwan MK. Application of the U-P Finite-Element Method to the Study of Articular-Cartilage. Journal of Biomechanical Engineering-Transactions of the Asme. 1991;113(4):397–403. [PubMed]
60. Wu JJ, Eyre DR, Slayter HS. Type VI collagen of the intervertebral disc. Biochemical and electron-microscopic characterization of the native protein. Biochem J. 1987;248(2):373–381. [PubMed]
61. Wu JZ, Herzog W. Analysis of the mechanical behavior of chondrocytes in unconfined compression tests for cyclic loading. J Biomech. 2006;39(4):603–616. [PubMed]
62. Wu JZ, Herzog W. Finite element simulation of location- and time-dependent mechanical behavior of chondrocytes in unconfined compression tests. Ann Biomed Eng. 2000;28(3):318–330. [PubMed]
63. Wu JZ, Herzog W, Epstein M. Modelling of location- and time-dependent deformation of chondrocytes during cartilage loading. J Biomech. 1999;32(6):563–572. [PubMed]
64. Yao H, Gu WY. Physical signals and solute transport in human intervertebral disc during compressive stress relaxation: 3D finite element analysis. Biorheology. 2006;43(3–4):323–335. [PubMed]
65. Yao H, Gu WY. Three-dimensional inhomogeneous triphasic finite-element analysis of physical signals and solute transport in human intervertebral disc under axial compression. J Biomech. 2007;40(9):2071–2077. [PMC free article] [PubMed]