Thin Bio-Artificial Tissues in Plane Stress: The Relationship between Cell and Tissue Strain, and an Improved Constitutive Model

INTRODUCTION

Bio-artificial tissues under development for the replacement of injured or diseased tissue in the human body must be not only biologically compatible, but also mechanically compatible. For this reason, they must reproduce the mechanical behavior of the healthy tissues they replace. The aim of this work is to develop an improved set of constitutive equations that describe how a class of bio-artificial tissues behaves mechanically, based upon an understanding of the mechanical properties of the tissue's constituents.

The improved constitutive model presented in this work is an extension of the Zahalak model (Zahalak et al., 2000). The primary limitation of Zahalak's model is its assumption that, irrespective of their relative mechanical properties, cells deform in registry with the tissue. We show that although this approximation is reasonable for tissues with extremely high cell concentrations and low elastic mismatch between cells and extracellular matrix, it needs to be addressed for other tissues. We address this by including a correction factor called the "strain factor". When a tissue is strained uniaxially, the strain factor is the ratio between the average strain along a cell's axis and the average tissue strain resolved in that direction.

With the adjustment we present, the Zahalak model applies to a broad range of tissues. However, in developing models for characterizing the strain factor, we refine our attention to a class of bio-artificial tissue constructs consisting of relatively stiff, elongated (length/width ratio on the order of 5-40) fibroblasts cultured in a relatively compliant reconstituted collagen matrix. These tissue constructs are far more compliant than most living tissues, due to the high compliance of reconstituted collagen; in such constructs, the matrix is much more compliant than fibroblasts (Wakatsuki et al., 2000; Zahalak et al., 2000).

We consider in these models only the short-term and very long-term response of these constructs, and thus treat the constructs as an incrementally linear elastic collagen matrix (e.g., Parry, 1988; Roeder et al., 2002) populated by perfectly bonded linear elastic cells. Although the instantaneous elastic moduli of collagen are strain dependent (Pryse et al., 2003; Ozerdem and Tozeren, 1995; Pins et al., 1997), the approximation of linearity is appropriate for small strain increments applied monotonically (Fung, 1981). A discourse on the limitations of modeling biological tissues with linear kinematics and Hooke's law is presented by Prager (1969). We further refine our attention to tissue constructs in which the cells have remodeled the matrix into a very thin membrane whose thickness is in the order of the cell diameter (see, for example, Wakatsuki et al., 2000).

The primary model we use for characterizing the strain factor is a scaling model. We validate the scaling model through comparison to numerical simulations and an exact solution for a special case, and then calibrate the model's single free parameter using Monte Carlo simulations involving several different types of computational analyses. The scaling model needs to account for very strong interactions between neighboring and overlapping cells, which can form tightly linked networks. To account for these effects, we incorporate a first-order statistical model, and validate the resulting "percolation model" with Monte Carlo simulations.

Background

From a mechanics viewpoint, this work is related to earlier work on composite materials reinforced with short fibers. A rich foundation of studies on short-fiber composites considers the problem from a great number of perspectives (e.g., Fukuda and Kawata, 1974; Chou, 1992; Budiansky and Cui, 1995; Tucker and Liang, 1999). This work makes use of homogenization procedures, which estimate overall mechanical properties based upon the solution of model problems, and unit cell approaches, which estimate the overall mechanical properties based on the response of idealized microstructural representations.

Several assumptions and modeling techniques used in this article have been employed in prior studies of the mechanical environment of cells, which have focused largely on relatively compliant, roughly spherical cells in cartilage. Baer and Setton (2000) treated the cells and matrix as linear to study the short-term and long-term mechanical environment of such cells. Wu et al. (1999) and Wu and Herzog (2000) apply both unit cell analysis and linear elastic homogenization theory to these tissues. More complicated constitutive models such as biphasic theory (Mow and Ratcliffe, 1997) have been applied to this problem by Bachrach et al. (1995) and Guilak and Mow (2000).

This study differs from earlier work in its focus on tissues containing oblate cells in a relatively compliant matrix. The following reviews the analytical foundation of the specific models used in this work.

Zahalak model

METHODS

This section describes the numerical and analytical models used to evaluate strain factors in two dimensions, and the update to the Zahalak constitutive model. The following section describes the numerical models, the idealizations of cells, and the limitations of the numerical models. The section "Analytical predictions for strain factors" describes the analytical and statistical models developed to predict strain factors and establish percolation thresholds. The section "Extension of the Zahalak model to incorporate strain factors" describes the way that strain factors are incorporated to update Zahalak's model.

Numerical models

Numerical analyses served three purposes: 1), to validate the concept of a strain factor; 2), to validate the scaling model over a broad range of material parameters and cell concentrations; and 3), to find the scaling model's one free parameter through Monte Carlo simulations. The analyses employed both a commercial finite element (FE) package (ADINA v.7.5.2), and a specialized FE code written with MATLAB.

The tissues considered were thin membranes subjected to uniaxial inplane stretching, which were modeled with plane stress conditions (Saada, 1993). Analyses all required a planar mathematical discretization of a region containing a prescribed number of cells with prescribed orientations (Figs. 1 and 2). Since some random distributions of cells required extensive statistical analyses, many different FE meshes were needed.

Boundary conditions simulated a tissue that was infinitely long in the direction of the applied stretch, and constrained from contracting in the direction perpendicular to this applied stretch. As depicted in Fig. 1, the top and bottom edges of all FE meshes were restrained from moving vertically, the left edges were restrained from moving horizontally, and the right edges were constrained to remain vertical while displacing. All edges were free of shear tractions. Analogous periodic boundary conditions were used in simulations involving applied shear strains.

The cells and matrix were parametrically assigned linear elastic, isotropic material properties. The matrix was assigned a Poisson's ratio of v = 0.49, as were the cells in simulations requiring this.

In the following sections, we describe the models for cells and their limitations; the model for cell distributions; the FE discretizations; and the specially written FE code we developed for large analyses.

Sparse cell populations

Models in which the cells occupied a 2D region with elastic modulus E^sub c^ were studied for comparison (Fig. 1 b). In these meshes, the cell was modeled with 2D elements like those of the matrix.

Convergence studies. The primary challenge in attaining convergence stemmed from the singularities in the elastic problems considered here: the elasticity solutions for ribbon cells and rectangular cells embedded in an elastic matrix predict infinite matrix strains at the end of the cell. FE analyses inherently smooth end effects over an area that is on the order of the element size.

We refined FE meshes until further refinement showed no effect on the strain factor. Convergence studies indicated accuracy on the order of a few percent; error decreased as the ratio of the cell and matrix moduli approached 1. All strain factors calculated using the FE method are upper bounds. FE models overpredict stiffness when the discretization and interpolation schemes used in the solution cannot precisely replicate the analytical displacement field. Elements used for the cells performed better than those used for the matrix because of the relatively uniform strain fields in the cells. The net result was a larger underprediction of strain in the matrix than in the cell, resulting in an overprediction of the strain factor.

Dense cell populations