- Email: [email protected]
Fractional-order elastic models of cartilage: A multi-scale approach
Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Fractional-order elastic models of cartilage: A multi-scale approach Richard L. Magin a,*, Thomas J. Royston b a b
University of Illinois at Chicago, 851 S. Morgan St., Chicago, IL 60607, USA University of Illinois at Chicago, 842 West Taylor St., Chicago, IL 60607, USA
a r t i c l e
i n f o
Article history: Received 22 April 2009 Received in revised form 4 May 2009 Accepted 4 May 2009 Available online 12 May 2009 PACS: 43.80.Cs 46.35.+z 87.61.c 45.10.Hj Keywords: Mathematical models Mechanical stress Medical systems Stress Tissues Hysteresis
a b s t r a c t The objective of this research is to develop new quantitative methods to describe the elastic properties (e.g., shear modulus, viscosity) of biological tissues such as cartilage. Cartilage is a connective tissue that provides the lining for most of the joints in the body. Tissue histology of cartilage reveals a multi-scale architecture that spans a wide range from individual collagen and proteoglycan molecules to families of twisted macromolecular fibers and fibrils, and finally to a network of cells and extracellular matrix that form layers in the connective tissue. The principal cells in cartilage are chondrocytes that function at the microscopic scale by creating nano-scale networks of proteins whose biomechanical properties are ultimately expressed at the macroscopic scale in the tissue’s viscoelasticity. The challenge for the bioengineer is to develop multi-scale modeling tools that predict the three-dimensional macro-scale mechanical performance of cartilage from micro-scale models. Magnetic resonance imaging (MRI) and MR elastography (MRE) provide a basis for developing such models based on the nondestructive biomechanical assessment of cartilage in vitro and in vivo. This approach, for example, uses MRI to visualize developing proto-cartilage structure, MRE to characterize the shear modulus of such structures, and fractional calculus to describe the dynamic behavior. Such models can be extended using hysteresis modeling to account for the non-linear nature of the tissue. These techniques extend the existing computational methods to predict stiffness and strength, to assess short versus long term load response, and to measure static versus dynamic response to mechanical loads over a wide range of frequencies (50–1500 Hz). In the future, such methods can perhaps be used to help identify early changes in regenerative connective tissue at the microscopic scale and to enable more effective diagnostic monitoring of the onset of disease. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction For the repair and regeneration of load-bearing tissues, such as bone, spinal disks, or articular cartilage, knowledge of the mechanical strength and stiffness is critical. Such biomechanical properties depend, in the case of cartilage, on the distribution and micro-scale architecture of the component elastomer and protein adhesion molecules (collagen, elastin, proteoglycan and hyaluronan) that form the extracellular matrix (ECM). Chondrocytes, the cells of cartilage that mold the ECM, function at the microscopic scale by creating nano-scale networks of elastomers, while the tissue biomechanics of cartilage is expressed at the macroscopic scale through the free movement and shock absorbing properties of the articulating joints of the body [1]. Consequently, in designing knee or hip joint replacements from either synthetic or natural cell-based materials,
* Corresponding author. Tel.: +1 312 996 2331. E-mail addresses: [email protected] (R.L. Magin), [email protected] (T.J. Royston). 1007-5704/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2009.05.008
658
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
the challenge for the bioengineer is to develop visualization and assessment tools that predict the three-dimensional (3D) macro-scale mechanical performance from noninvasive micro- and nano-scale observations and measurements. In this paper, we describe novel multi-scale viscoelastic models for the design and in vivo analysis of engineered chondrogenic tissue using magnetic resonance imaging (MRI) and magnetic resonance elastography (MRE). MRI and MRE are 3D imaging methods that generate spatial maps of tissue relaxation times, diffusion coefficients, magnetization transfer ratios and estimates of the mechanical shear moduli. In the past we have applied MRI and MRE to monitor the growth of tissue engineered bone, fat and cartilage [2–4]. These previous studies demonstrated a reduction in T1 and T2 relaxation times during the stages of osteogenesis [4] and a decrease in tissue stiffness, but not T2 with adipogenesis in tissue-engineered fat [5]. This initial work highlighted the basic capabilities of microscopic MRI to measure the tissue T1 and T2 relaxation times and MRE to determine material stiffness. Our long term research goal – outlined in Fig. 1 – is to extend these methods to assess the structure and function of other complex multi-scale tissue engineered biomaterials at the nano-, micro- and macro-scales. However, in order to accurately interpret the microscopic basis for macroscopic measurements, an improved understanding of engineered tissue mechanics is needed. Conventional linear, fixed element macro-scale models (e.g., Voigt, Maxwell, etc.) are insufficient for describing many phenomena observed in biomaterials [1]. New models that account for behavior arising from micro-scale phenomena, including memory, hysteresis and fractional order dynamics (fractional calculus) need to be formulated to computationally predict macro-scale behavior. Fractional calculus models represent a relatively simple way to describe dynamics in complex, porous, or composite systems, such as cartilage [6–10]. Hence, they are finding increasing use in many areas of science and engineering from the nano- to the macro- scale. There is a multi-scale generalization inherent in the definition of the fractional derivative that accurately represents interactions occurring over wide dynamic ranges of space or time [6,7]. Thus, there is no need to segment or compartmentalize systems into thousands of subsystems or subunits – a system generalization that often creates more complexity than can be experimentally evaluated. Also, experiments and measurements conducted on complex materials can shed light on the meaning of fractional order operations, when the results match the intermediate order system dynamics predicted by fractional calculus [8,9]. Using fractional calculus, we can begin to unravel the contributions to the physics that follows from the fundamental model dynamics, the geometry, the interaction between system components, and the physical barriers [10]. Fractional order models often work well, particularly for dielectrics and viscoelastic materials over extended ranges of time and frequency [8,9]. In heat transfer and electrochemistry, for example, the half order fractional integral is the natural integral operator connecting the applied gradients (thermal or material) with the diffusion of ions of heat [6,10]. Thus, a focus of this current work is to investigate dynamic processes in cartilage using MRI and MRE and to apply fractional order derivative models where they are useful and appropriate. Articular cartilage is a thin heterogeneous tissue (with a typical thickness less than 2 mm) that coats the skeletal joint surfaces. Its physical and chemical integrity are essential for maintaining normal joint movements. Studies of the anatomy of articular cartilage show an organized structure, consisting of four distinct zones: (i) superficial zone with collagen fibrils (macromolecular arrangements of collagen fibers) oriented along the cartilage surface (5–10% of thickness in human), (ii) transitional zone with randomly oriented collagen (40–45% of thickness), (iii) radial zone with radially oriented collagen (40–45% of thickness) and (iv) calcified zone acting as a transition region from soft hyaline cartilage to subchondral bone (5–10% of thickness). Each zone has different biomechanical properties, which largely reflect the composition of the ECM [11]. All zones; however, are avascular, lack sensory neurons, and contain only a few mature chondrocytes, whose low met-
Fig. 1. A multi-scale research paradigm for cartilage tissue engineering with emphasis on the structure and composition of the extracellular matrix (ECM).
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
659
abolic activity precludes rapid proliferation [12]. In addition, the proteoglycans of cartilage inhibit cell adhesion. Thus, in the adult, articular cartilage has a very limited capability for self-repair. 2. Cartilage tissue engineering Cartilage tissue engineering begins by introducing cells to the injured site, usually embedded in a porous polymer matrix or scaffold, and proceeds with the administration of specific tissue growth and differentiation factors [13]. The ultimate goal is for the implant to take over the role of ECM production and thus to restore the function of damaged cartilage. Natural and synthetic scaffolds are available in various physical forms, but so far naturally derived collagen, the major structural component of cartilage, has received the most attention. This is largely due to the fact that cartilage can be recognized by cellular enzymes and be remodeled as needed [14]. Both chondrocytes and mesenchymal stem cells are used as cell sources in cartilage tissue engineering. The advantage of using chondrocytes is that the cells do not have to differentiate before they can secrete cartilage-specific ECM components (e,g., collagen, proteoglycans), and hence are able to produce cartilage-like tissue in a short time (<2 months) [15]. However, chondrocytes have the limitations of low expandability and complexity in their interactions with the dense connective tissue membrane (periosteum) that lines the surface of all bones [16]. Also, in clinical situations, a large number of donor chondrocytes are not always available. In contrast, stem cells are relatively easy to harvest, readily expandable from small quantities in vitro in cell culture, and reliably differentiate into chondrogenic cells [17]. Even though cartilage is a relatively uncomplicated tissue – given its low cell density, limited cell diversity and the absence of vascular structure or nerve supply – so far tissue engineering techniques have failed to produce a mature tissue with all the structural and mechanical properties needed for it to function in vivo at the joint surface [18,19]. Therefore, there is a need to provide the tissue engineer with additional information about the structure and composition of the ECM in engineered cartilage; information that can be used to optimize the selection of materials and administration of growth factors. 3. Magnetic resonance elastography The stiffness of biological tissue is a key regulator of cell differentiation that modulates embryonic development and tissue regeneration. Mechanical properties such as stiffness (shear modulus and viscosity) can be measured by using MRI to visualize acoustic shear wave motion in the tissue [20]. In this technique (commonly referred to as magnetic resonance elastography, MRE) the mechanical vibrations in the tissue are viewed in a phase contrast MR image that encodes the spatial and temporal pattern of strains associated with the propagation of the acoustic waves. These image data are combined with a viscoelastic model of the tissue to construct a map of the ‘‘shear stiffness” throughout the tissue. Mechanical shear waves, typically with amplitudes of less than 100 lm, and with frequencies of 100–1000 Hz, are generated in tissue by a piezoelectric or a speaker coil oscillator directly coupled to the region of interest. By using multiple phase offsets and motion encoding gradients, the shear wave images are reconstructed, from which the local values of the tissue viscoelastic properties can be calculated [20]. MRE studies using clinical 1.5 T MRI systems have established that changes in the mechanical properties of tissues are associated with developing disease: malignant tumors appear to be stiffer than benign tumors; fibrosis and cirrhosis tend to increase liver stiffness; and articular cartilage softens in developing osteoarthritis. Recent micro-MRE studies at 11.7 T have achieved higher resolution: generating stiffness maps of small biological specimens with a higher in-plane resolution (tens of microns, depending on the MR system resolution) and extending measurements of the shear modulus up to 100 kPa [2]. A necessary step in constructing the stiffness map of a tissue from the MRE data is the assumption of a linear viscoelastic model for the tissue. In the simplest case a single spring (Hookean) or a two element parallel spring-dashpot (Voigt) constitutive model are used. In addition to assuming linearity, many techniques have been proposed to solve the inverse problem of estimating the tissue elastic modulus, E, and viscosity, f, from MRE measurements of wavelength and attenuation coefficient. Some techniques use algorithms that make an assumption of tissue incompressibility, and then use a spatial filter to remove displacements caused by compression wave motion. Other MRE reconstruction techniques are based on an algebraic inversion of the differential equation (AIDE) for the medium (e.g., [21]). AIDE techniques directly estimate elastic modulus and attenuation with data smoothing to minimize noise problems. AIDE techniques may also assume incompressibility to reduce complexity and noise susceptibility of the inversion. Other techniques use finite element (FE) analysis iteratively within the material identification algorithm or apply other model-based reconstructions (e.g., [22]). In our work we have to date used the two element Voigt model, but one of the goals of the current work is to employ new multi-scale constitutive models of cartilage to better understand how microscopic MRE measurements can be related to macro-scale biomechanics. 4. Cartilage tissue multi-scale models In order to correctly interpret MRE measurements of cartilage in terms of its micro-scale tissue structure a better understanding of multi-scale tissue mechanics is needed. Recent work [23] has shown that the macro-scale mechanics of cartilage does not follow directly from the mechanics of the individual fibrils, but that the tissue’s strength strongly depends on the arrangement of the fibrils at the micro-scale. Cartilage tissue also shows hysteresis at the macro-scale, which is not evident in isolated fibril mechanics, but rather is indicative of a friction-like sliding motion of the tissue components [1]. Finally, the
660
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
mechanical properties of cartilage are time and history dependent, that is, the observed strain, e(t), will be different when the rate at which stress, r(t), is applied changes and when different stress paths are followed to arrive at the same final load stress. While the dynamics of individual collagen fibrils at the micro-scale may be reasonably modeled using conventional viscoelastic approaches incorporating Hookean (r(t) = E e(t)) and Newtonian fluid elements (r(t) = gde(t)/dt) for elasticity (E) and viscosity (g), respectively (e.g., Voigt, Maxwell, Kelvin models, etc.), extending this to the macro-scale with many fibrils arranged in complex orientations, where the orientation also affects the dynamics, requires some form of spatial averaging and approximation. Considering the problem to be linear (neglecting hysteresis), this is precisely the type of problem for which fractional calculus techniques have shown great promise [24,25]. If, in addition, non-linear hysteretic behavior can be incorporated into the fractional calculus model, then we can develop simplified mathematical expressions, of hysteresis models for cartilage that have a range of applications. Current work in applying fractional calculus to viscoelastic theory starts with the idea that it is the order of the derivative of the strain that – to a first approximation – characterizes the material’s behavior (assuming a one-dimensional stress– strain relation, r(t) = E*dae(t)/dta). In this model, the order of the derivative is zero for a Hookean solid (E* = E) and to one for a Newtonian fluid (E* = g). Viscoelastic materials occupy the intermediate range with a fractional order ‘‘a” between zero and one. Thus, one can build a multi-component fractional equivalent of the ‘‘standard” linear solid (or fluid) by replacing one or more springs and dashpots with ‘‘Springpots”. There are a number of papers (and one monograph) that outline this approach [24,6]. The rationale for the fractional calculus approach is that it encompasses a wider range of observed viscoelastic data (in the time or frequency domain) in a simpler (fewer elements), yet linear model. Many complex, composite and biomaterials exhibit fractional order behavior through power law decays in their creep or relaxation response and their dynamic elastic modulus. Alternatively, quasi-linear elastic and non-linear elastic models has been used to analyze these cases. In our work to ‘‘span” the nano-, micro- and meso-scales we choose to apply the fractional calculus models. Such fractional order constitutive models can also be applied to our elastography experiments. Note in particular, Coimbra and coworkers [26–28], however, go further; they generalize the fractional derivative parameter ‘‘a” so that it becomes a function of time: a(t). Such time varying derivative operators have only recently been introduced (e.g., [29]) and are little studied. There are many issues still to be resolved; but, this model of linear viscoelastic materials provides a bridge between the different fractional calculus models, and it can perhaps also include strain rate and strain path dependence in the formalism. Briefly, Coimbra and coworkers show that the so-called variable order definition of the fractional derivative is simply a properly weighted sum of a large number of fractional order derivative terms. This extends the idea that a single fractional order derivative is the sum of properly weighted integer order spring and dashpot elements [29]. In addition, the variable order fractional derivative model of stress/strain is postulated to be directly related to the rate of change of the long range order of the molecules within a viscoelastic material. In the case of simple compression of a carbon–epoxy composite, Ramirez and Coimbra [28] demonstrate that a(t*) is proportional to an ‘‘entropy” like function t*ln(t*) where t* is a dimensionless time for compression experiments spanning a wide range of strain rates. For example, in this paper, they present a normalized stress versus strain plot as a single curve that extends over eight orders of magnitude in strain rate using only one parameter. Thus, the variable order fractional derivative yields a constitutive model that on the one hand connects with the underlying molecular disorder (disorder that evolves with compression of the composite), and on the other hand with a macroscopic description of mechanical state through an extended viscoelastic model of the physical specimens under study. In our work, we apply these approaches to model complex biomaterials, which exhibit complex viscoelastic behavior. The hypothesis is that fractional and perhaps variable order calculus will provide a greater understanding of the molecular events that occur as cartilage grows, ages and responds to injury. Such understanding can lead to a simpler description of the underlying multi-scale processes that occur when composite materials, such as cartilage, are stressed or strained. Now, returning to our consideration of hysteresis modeling, recently Royston [30] established the relationship of an Iwan hysteresis model consisting of series Coulomb friction and elastic slide components arranged in parallel, to several other models of hysteretic behavior that have been used to simulate a wide range of physical phenomena, most significantly the classical Preisach model. These hysteresis models all possess the necessary and sufficient properties of ‘‘wiping-out” and ‘‘congruency”. Wiping-out means that the current state is only dependent on past extrema; congruency means that the response to a cyclic input is identical in form (forms the same shape hysteresis curve) regardless of previous extrema, which only serve to shift the entire output loop. It was then possible to extend analytical tools developed for the Preisach model to all of the models with which it is related, and vice versa. Such tools include experimental identification, inversion and analysis of dynamic energy flow and dissipation. This prior work provides a framework upon which there can be extensions in several directions relevant to biological tissues. For example, while Royston [30] focused on one-dimensional behavior, coupled multi-dimensional anisotropic behavior could be considered. Given the developed links between the elasto-slide and Preisach model formulations, a logical approach may be to extend techniques developed for multi-dimensional vector-based Preisach models to the anisotropic multi-scale problem [31]. In a more general sense, since hysteresis in and of itself, is actually a manifestation of the presence of multiple scales in a system [32], it is appropriate that it be included as a part of a multi-scale modeling approach to elastic biomaterials.
661
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
Fractional calculus models provide a relatively simple way to describe the biomechanics of complex, porous, or composite systems. There is a multi-scale generalization inherent in the definition of the fractional derivative that accurately represents interactions occurring over wide dynamic ranges of space or time. Thus, there is no need to segment or compartmentalize systems into subsystems or subunits all the way down to the molecular level – a process that often creates more complexity than can be experimentally evaluated. Also, viscoelastic experiments conducted on complex materials can shed light on the ‘‘meaning” of fractional order operations; when the results match the predictions of fractional order system dynamics the physics and the mathematics are in correspondence. Thus, we can begin to unravel the contributions to the physics that follow from the fundamental model dynamics, the geometry and the interaction between system components and physical barriers. In the following, we will present some initial steps taken to apply this methodology to the viscoelastic behavior of cartilage using MRI and MRE. 5. Biological tissue measurements MRI examination of tissue-engineered cartilage provides a non-invasive way to visualize changes in structure while MRE measures the inherent elastic properties (complex modulus) of the developing tissue. The objective of this ongoing investigation is to determine the key MR parameters that describe the essential structural characteristics of new cartilage tissue derived from engineered human mesenchymal stem cells (hMSC). The hMSCs used in this study were isolated from commercially available bone marrow, expanded in cell culture and seeded in porous poly(ethylene glycol) scaffold. The tissue engineered construct was exposed to chondrogenic induction media for 4 weeks and samples evaluated each week (Fig. 2). All samples were sealed in glass NMR tubes filled with media and loaded into a 10-mm diameter RF saddle coil for the MRI studies. T2 maps were calculated each week using a spin echo imaging sequence to acquire 32 echoes with a 7 ms echo spacing for the chosen axial slice. The data were analyzed to produce a pixel-by-pixel T2 map for the entire sample scope (Fig. 2). As the cartilage tissue developed contrast appeared throughout the slice and a monotonic increase in T2 was observed. To confirm that the T2 changes reflect progression of the hMSC along the developmental pathway to form cartilage, samples were removed for histochemical analysis: safranin O and alcian blue staining, collagen Type II immunohistochemistry, and quantitative glycosaminoglycan (GAG) analyses. The constructs stained positive for collagen Type II and GAGs as time increased from week 0 to week 4. Furthermore, the quantitative GAG assay established a larger increase of GAG in the experimental constructs relative to the controls [33]. Other recent studies of tissue-engineered cartilage examined other MR-based parameters: T1 relaxation in the rotating frame (T1q) and magnetization transfer (MT). T1q and MT are magnetic resonance derived parameters known to be sensitive to tissue composition and structure. In brief, both reflect the extent of water binding to proteins and other macromolecules. Our hypothesis for this study is that developing tissue will display differences in the relative amount of free and ‘‘so-called” bound water. The bound protein fraction or BPF can be derived from the MT data. Previous studies of bone and adipose tissues showed that the tissue water relaxation times (T1 and T2) do not always show early changes in tissue composition; changes that need to be measured in order for tissue engineers to optimize tissue growth. The bar graph shown in Fig. 3
Cartilage
MRI
Safranin-O
Superficial zone Superficial zone
Transitional zone
Chondrocytes
Radial zone
Proteoglycan
Transitional zone Radial zone Calcified zone
Calcified zone
Bone Marrow
140 120 100 80 60 40
Time (wk)
10 mm
0
1
2
3
4
Fig. 2. (Top) The three images of cartilage are: an anatomical drawing, a high resolution MRI, and a histological image (safranin-O/fast green-stained section (40)). In the MR image, three zones can be easily identified. MR acquisition parameters were: TR/TE = 1000/ 60 ms; field-of-view (FOV) = 0.6 cm; slice thickness = 0.5 mm; in-plane resolution = 22.7 lm 22.7 lm; and 64 averages [3]. (Bottom) Five T2 maps for engineered cartilage over a 4 week growth period. T2 contrast and a monotonic increase of T2 is observed with time.
662
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664 180 Control
160
Stimulated
140
T1p in milli-sec
120 100 80 60 40 20 0 Day 0
Day 3
Day 6
Day 9
Day 12
Fig. 3. T1q values derived from T1q weighted images of tissue-engineered cartilage. The data are means ± one standard deviation for four measurements.
shows the significant changes in the T1q parameter over the duration of this study. A similar change was also observed in BPF (see Fig. 4). In a tissue engineering study by Marion and coworkers [5] the growth stages of engineered fat were monitored using MR measurements of relaxation times and the water diffusion coefficient. In this study of MSC derived adipogenic cells the T2 values stayed unchanged for 4 weeks, but MR elastography measurements of the bulk modulus showed a significant decrease (Fig. 5). By week 4, for example, the shear modulus had decreased to less than one fourth of the initial value (5.8 to 1.4 kPa). The decrease in shear stiffness indicates that the engineered adipogenic tissue becomes ‘softer’ as adipogenic matrix production gradually increases. These results illustrate the complementary role of MRI and MR elastography in assessing the development of engineered tissue in vitro. 6. Combined hysteresis and fractional calculus models In the generalized elasticity model of tendons and ligaments [1] a number, say N, of springs are connected in parallel. An applied force induces a stress that as it grows engages a sequence of spring elements beginning with the weakest (lowest elastic constant En) extending first and the strongest (largest elastic constant EN) last. The resulting stress–strain behavior is not a straight line, but a curve with an initial flat section and an increasing slope. We propose to extend this idea by using a parallel combination of N individual Coulomb sliding friction element and one or more fractional order springpots (Fig. 5). The sliding friction Iwan parallel-series hysteresis model captures the rate-independent hysteresis behavior, while the springpot captures the power law dynamics observed in tissue mechanics.
Fig. 4. T2 relaxation time and shear stiffness data for engineered fat.
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
663
μh φ1
φ2
σ[t] μh
φm φN
μh εN[t,φN] ε[t]
Fig. 5. Combined hysteresis and fractional model.
Under quasistatic loading conditions, the springpot asymptotically carries no load and hysteretic relation between stress
r and strain e is captured via the elasto-slide elements [30]. Under dynamic loading conditions, elasto-slide element behavior is unchanged. However, now the springpot, acting in parallel, introduces a rate-dependent stress/strain component. In the time domain, the stress–strain relation for the hysteretic component can be expressed as follows. Consider a finite number N of parallel-series elasto-slide elements. Referring to Fig. 5 and separating the stress related to hysteresis rH[t] from that of the dashpot rF[t], such that r[t] = rH[t] + rF[t], we have: N X
lh ðe em ðt; /m ÞÞ ; lh /m sign½e_ m¼1 je em ðt; /m Þj 6 /m unchanged em ½t; /m ¼ : if je em ðt; /m Þj ¼ /m e /m sign½e_ rH ½t ¼
rm ½t; rm ½t ¼
ð1Þ
In the time domain, the stress–strain relation for the fractional component, or springpot, has dynamics that can be expressed as, a
rF ðtÞ ¼ Esa
d eðtÞ a ; dt
ð2Þ
where a can take values between 0 and 1. Here a = 0 represents a linear spring and a = 1 represents a damper (Es = g). Any specific value of a can be represented with hierarchical arrangements of springs and dashpots, such as ladders, trees and fractal networks [6]. 7. Conclusions These new multiple scale techniques extend the existing computational methods for predicting stiffness and strength, assessing short versus long term load response, and measuring static versus dynamic response to mechanical loads over a wide range of frequencies. While the specific application considered in this project involves tissue-engineered cartilage, the multi-scale computational methodology used to connect nano-scale structure with macro-scale function is relevant to a much wider range of biomechanical systems and biomedical problems. Our new methods could, for example, help identify early microscopic changes in connective tissue at the microscopic scale. Knowledge of such changes would enable early diagnostic monitoring for the onset of disease and better assessment of the effectiveness of new drugs or therapies. References [1] Ethier CR, Simmons CA. Introductory biomechanics: from cells to organisms. Cambridge, UK: Cambridge University Press; 2007. [2] Othman SF, Xu H, Royston TJ, Magin RL. Microscopic magnetic resonance elastography (microMRE). Magn Reson Med 2005;54:605–14. [3] Othman SF, Li J, Abdullah O, Moinness JJ, Magin RL, Muehleman C. High-resolution/high-contrast MRI of human articular cartilage lesions. Acta Orthop 2007;78:536–46. [4] Xu H, Othman SF, Liu H, Peptan IA, Magin RL. Magnetic resonance microscopy for monitoring osteogenesis of tissue-engineered constructs in vitro. Phys Med Biol 2006;51:719–32. [5] Marion N, Othman SF, Stosich MS, Magin RL, Mao JJ. Measurement of tissue engineered adipogenic constructs using magnetic resonance elastography. Abstr Annu Meet BMES 2006. [6] Magin RL. Fractional calculus in bioengineering. Redding, CT: Begell House; 2006. [7] Hilfer R. Fractional calculus in bioengineering. River Edge, NJ: World Scientific; 2000. [8] Grimnes S, Martinsen OG. Bioimpedance and bioelectricity basics. San Diego, CA: Academic Press; 2000. [9] Lakes RS. Viscoelastic solids. Boca Raton, FL: CRC Press; 1999. [10] Bard AJ, Faulkner LR. Electrochemical methods: fundamentals and applications. second ed. New York, NY: John Wiley & Sons; 2001. [11] Meachim G, Stockwell RA, The matrix, in: M.A.R. Freeman (Ed.), Adult articular cartilage, Pitman Medical, Kent; 1979. p. 1–50. [12] Hunziker EB. Articular cartilage repair: basic science and clinical progress, a review of the current status and prospects. Osteoarthr Cartil 2002;10:432–63. [13] Pietrzak WS. Musculoskeletal tissue regeneration: biological materials and methods. Totowa, NJ: Humana Press; 2008.
664
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
[14] Speer DP, Chvapil M, Volz RG, Holmes MD. Enhancement of healing in osteochondral defects by collagen sponge implants. Clin Orthop Relat Res 1979;144:326–35. [15] Elisseeff JH, Puleo C, Yang F, Sharma B. Advances in skeletal tissue engineering with hydrogels. Orthod Craniofac Res 2005;8:150–61. [16] Temenoff JS, Mikos AG. Review: tissue engineering for regeneration of articular cartilage. Biomaterials 2000;21:431–40. [17] Mao JJ. Stem-cell-driven regeneration of synovial joints. Biol Cell 2005;97:289–301. [18] Caplan AI. Tissue engineering designs for the future: new logics, old molecules. Tissue Eng 2000;6:1–8. [19] Sharma B, Elisseeff JH. Engineering structurally organized cartilage and bone tissues. Ann Biomed Eng 2004;32:148–59. [20] Muthupillai R, Lomas DJ, Rossman PJ, Greenleaf JF, Manduca A, Ehman RL. Magnetic resonance elastography by direct visualization of propagating acoustic strain waves. Science 1995;269:1854–7. [21] Oliphant TE, Manduca A, Ehman RL, Greenleaf JF. Complex-valued stiffness reconstruction for magnetic resonance elastography by algebraic inversion of the differential equation. Magn Reson Med 2001;45:299–310. [22] Van Houten EEW, Miga MI, Weaver JB, Kennedy FE, Paulsen KD. Three-dimensional subzone-based reconstruction algorithm for MR elastography. Magn Reson Med 2001;45:827–37. [23] Chandran PL, Barocas VH. Affine versus non-affine fibril kinematics in collagen networks: theoretical studies of network behaviour. ASME J Biomech Eng 2006;128:259–70. [24] Magin RL. Fractional calculus in bioengineering. Crit Rev Biomed Eng 2004;32(1–4):1–377. [25] Sinkus RK, Siegmann T, Xydeas M, Tanter C, Claussen C, Fink M. MR elastography of breast lesions: understanding the solid/liquid duality can improve the specificity of contrast-enhanced MR mammography. Magn Reson Med 2007;58:1135–44. [26] Coimbra CFM. Mechanics with variable-order differential operators. Ann Phys 2003;12:692–703. [27] Soon CM, Coimbra CFM, Kobayashi MH. The variable elasticity operator. Ann Phys 2005;14:378–89. [28] Ramirez LES, Coimbra CFM, A constitutive equation for linear viscoelastic thermoset materials undergoing compression, In: Proceedings of the of ASME IDETC/CIE conference. Paper No. DETC2007-34169; 2007. [29] Lorenzo CF, Hartley TT. Variable order and distributed order fractional operators. Nonlinear Dyn 2002;29:57–98. [30] Royston TJ. Leveraging the equivalence of hysteresis models from different fields for analysis and numerical simulation of jointed structures. ASME J Comput Nonlinear Dyn 2008;3:031006-1–6-8. [31] Mayergoyz ID. Mathematical models of hysteresis and their applications. NY: Elsevier, Academic Press; 2003. [32] Banks HT, Luke NS, Samuels Jr, JR, Viscoelasticity in polymers: phenomenological to molecular mathematical modelling. Numer Methods Partial Differ Eqn. doi:10.1002/num.20250; 2008. [33] Troken AJ, Abdullah OM, Othman SF, Magin RL, Mao JJ. Magnetic resonance imaging for monitoring chondrogenesis of stem-cell based constructs in vitro. Abstr Annu Meet BME 2006.
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Fractional-order elastic models of cartilage: A multi-scale approach Richard L. Magin a,*, Thomas J. Royston b a b
University of Illinois at Chicago, 851 S. Morgan St., Chicago, IL 60607, USA University of Illinois at Chicago, 842 West Taylor St., Chicago, IL 60607, USA
a r t i c l e
i n f o
Article history: Received 22 April 2009 Received in revised form 4 May 2009 Accepted 4 May 2009 Available online 12 May 2009 PACS: 43.80.Cs 46.35.+z 87.61.c 45.10.Hj Keywords: Mathematical models Mechanical stress Medical systems Stress Tissues Hysteresis
a b s t r a c t The objective of this research is to develop new quantitative methods to describe the elastic properties (e.g., shear modulus, viscosity) of biological tissues such as cartilage. Cartilage is a connective tissue that provides the lining for most of the joints in the body. Tissue histology of cartilage reveals a multi-scale architecture that spans a wide range from individual collagen and proteoglycan molecules to families of twisted macromolecular fibers and fibrils, and finally to a network of cells and extracellular matrix that form layers in the connective tissue. The principal cells in cartilage are chondrocytes that function at the microscopic scale by creating nano-scale networks of proteins whose biomechanical properties are ultimately expressed at the macroscopic scale in the tissue’s viscoelasticity. The challenge for the bioengineer is to develop multi-scale modeling tools that predict the three-dimensional macro-scale mechanical performance of cartilage from micro-scale models. Magnetic resonance imaging (MRI) and MR elastography (MRE) provide a basis for developing such models based on the nondestructive biomechanical assessment of cartilage in vitro and in vivo. This approach, for example, uses MRI to visualize developing proto-cartilage structure, MRE to characterize the shear modulus of such structures, and fractional calculus to describe the dynamic behavior. Such models can be extended using hysteresis modeling to account for the non-linear nature of the tissue. These techniques extend the existing computational methods to predict stiffness and strength, to assess short versus long term load response, and to measure static versus dynamic response to mechanical loads over a wide range of frequencies (50–1500 Hz). In the future, such methods can perhaps be used to help identify early changes in regenerative connective tissue at the microscopic scale and to enable more effective diagnostic monitoring of the onset of disease. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction For the repair and regeneration of load-bearing tissues, such as bone, spinal disks, or articular cartilage, knowledge of the mechanical strength and stiffness is critical. Such biomechanical properties depend, in the case of cartilage, on the distribution and micro-scale architecture of the component elastomer and protein adhesion molecules (collagen, elastin, proteoglycan and hyaluronan) that form the extracellular matrix (ECM). Chondrocytes, the cells of cartilage that mold the ECM, function at the microscopic scale by creating nano-scale networks of elastomers, while the tissue biomechanics of cartilage is expressed at the macroscopic scale through the free movement and shock absorbing properties of the articulating joints of the body [1]. Consequently, in designing knee or hip joint replacements from either synthetic or natural cell-based materials,
* Corresponding author. Tel.: +1 312 996 2331. E-mail addresses: [email protected] (R.L. Magin), [email protected] (T.J. Royston). 1007-5704/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2009.05.008
658
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
the challenge for the bioengineer is to develop visualization and assessment tools that predict the three-dimensional (3D) macro-scale mechanical performance from noninvasive micro- and nano-scale observations and measurements. In this paper, we describe novel multi-scale viscoelastic models for the design and in vivo analysis of engineered chondrogenic tissue using magnetic resonance imaging (MRI) and magnetic resonance elastography (MRE). MRI and MRE are 3D imaging methods that generate spatial maps of tissue relaxation times, diffusion coefficients, magnetization transfer ratios and estimates of the mechanical shear moduli. In the past we have applied MRI and MRE to monitor the growth of tissue engineered bone, fat and cartilage [2–4]. These previous studies demonstrated a reduction in T1 and T2 relaxation times during the stages of osteogenesis [4] and a decrease in tissue stiffness, but not T2 with adipogenesis in tissue-engineered fat [5]. This initial work highlighted the basic capabilities of microscopic MRI to measure the tissue T1 and T2 relaxation times and MRE to determine material stiffness. Our long term research goal – outlined in Fig. 1 – is to extend these methods to assess the structure and function of other complex multi-scale tissue engineered biomaterials at the nano-, micro- and macro-scales. However, in order to accurately interpret the microscopic basis for macroscopic measurements, an improved understanding of engineered tissue mechanics is needed. Conventional linear, fixed element macro-scale models (e.g., Voigt, Maxwell, etc.) are insufficient for describing many phenomena observed in biomaterials [1]. New models that account for behavior arising from micro-scale phenomena, including memory, hysteresis and fractional order dynamics (fractional calculus) need to be formulated to computationally predict macro-scale behavior. Fractional calculus models represent a relatively simple way to describe dynamics in complex, porous, or composite systems, such as cartilage [6–10]. Hence, they are finding increasing use in many areas of science and engineering from the nano- to the macro- scale. There is a multi-scale generalization inherent in the definition of the fractional derivative that accurately represents interactions occurring over wide dynamic ranges of space or time [6,7]. Thus, there is no need to segment or compartmentalize systems into thousands of subsystems or subunits – a system generalization that often creates more complexity than can be experimentally evaluated. Also, experiments and measurements conducted on complex materials can shed light on the meaning of fractional order operations, when the results match the intermediate order system dynamics predicted by fractional calculus [8,9]. Using fractional calculus, we can begin to unravel the contributions to the physics that follows from the fundamental model dynamics, the geometry, the interaction between system components, and the physical barriers [10]. Fractional order models often work well, particularly for dielectrics and viscoelastic materials over extended ranges of time and frequency [8,9]. In heat transfer and electrochemistry, for example, the half order fractional integral is the natural integral operator connecting the applied gradients (thermal or material) with the diffusion of ions of heat [6,10]. Thus, a focus of this current work is to investigate dynamic processes in cartilage using MRI and MRE and to apply fractional order derivative models where they are useful and appropriate. Articular cartilage is a thin heterogeneous tissue (with a typical thickness less than 2 mm) that coats the skeletal joint surfaces. Its physical and chemical integrity are essential for maintaining normal joint movements. Studies of the anatomy of articular cartilage show an organized structure, consisting of four distinct zones: (i) superficial zone with collagen fibrils (macromolecular arrangements of collagen fibers) oriented along the cartilage surface (5–10% of thickness in human), (ii) transitional zone with randomly oriented collagen (40–45% of thickness), (iii) radial zone with radially oriented collagen (40–45% of thickness) and (iv) calcified zone acting as a transition region from soft hyaline cartilage to subchondral bone (5–10% of thickness). Each zone has different biomechanical properties, which largely reflect the composition of the ECM [11]. All zones; however, are avascular, lack sensory neurons, and contain only a few mature chondrocytes, whose low met-
Fig. 1. A multi-scale research paradigm for cartilage tissue engineering with emphasis on the structure and composition of the extracellular matrix (ECM).
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
659
abolic activity precludes rapid proliferation [12]. In addition, the proteoglycans of cartilage inhibit cell adhesion. Thus, in the adult, articular cartilage has a very limited capability for self-repair. 2. Cartilage tissue engineering Cartilage tissue engineering begins by introducing cells to the injured site, usually embedded in a porous polymer matrix or scaffold, and proceeds with the administration of specific tissue growth and differentiation factors [13]. The ultimate goal is for the implant to take over the role of ECM production and thus to restore the function of damaged cartilage. Natural and synthetic scaffolds are available in various physical forms, but so far naturally derived collagen, the major structural component of cartilage, has received the most attention. This is largely due to the fact that cartilage can be recognized by cellular enzymes and be remodeled as needed [14]. Both chondrocytes and mesenchymal stem cells are used as cell sources in cartilage tissue engineering. The advantage of using chondrocytes is that the cells do not have to differentiate before they can secrete cartilage-specific ECM components (e,g., collagen, proteoglycans), and hence are able to produce cartilage-like tissue in a short time (<2 months) [15]. However, chondrocytes have the limitations of low expandability and complexity in their interactions with the dense connective tissue membrane (periosteum) that lines the surface of all bones [16]. Also, in clinical situations, a large number of donor chondrocytes are not always available. In contrast, stem cells are relatively easy to harvest, readily expandable from small quantities in vitro in cell culture, and reliably differentiate into chondrogenic cells [17]. Even though cartilage is a relatively uncomplicated tissue – given its low cell density, limited cell diversity and the absence of vascular structure or nerve supply – so far tissue engineering techniques have failed to produce a mature tissue with all the structural and mechanical properties needed for it to function in vivo at the joint surface [18,19]. Therefore, there is a need to provide the tissue engineer with additional information about the structure and composition of the ECM in engineered cartilage; information that can be used to optimize the selection of materials and administration of growth factors. 3. Magnetic resonance elastography The stiffness of biological tissue is a key regulator of cell differentiation that modulates embryonic development and tissue regeneration. Mechanical properties such as stiffness (shear modulus and viscosity) can be measured by using MRI to visualize acoustic shear wave motion in the tissue [20]. In this technique (commonly referred to as magnetic resonance elastography, MRE) the mechanical vibrations in the tissue are viewed in a phase contrast MR image that encodes the spatial and temporal pattern of strains associated with the propagation of the acoustic waves. These image data are combined with a viscoelastic model of the tissue to construct a map of the ‘‘shear stiffness” throughout the tissue. Mechanical shear waves, typically with amplitudes of less than 100 lm, and with frequencies of 100–1000 Hz, are generated in tissue by a piezoelectric or a speaker coil oscillator directly coupled to the region of interest. By using multiple phase offsets and motion encoding gradients, the shear wave images are reconstructed, from which the local values of the tissue viscoelastic properties can be calculated [20]. MRE studies using clinical 1.5 T MRI systems have established that changes in the mechanical properties of tissues are associated with developing disease: malignant tumors appear to be stiffer than benign tumors; fibrosis and cirrhosis tend to increase liver stiffness; and articular cartilage softens in developing osteoarthritis. Recent micro-MRE studies at 11.7 T have achieved higher resolution: generating stiffness maps of small biological specimens with a higher in-plane resolution (tens of microns, depending on the MR system resolution) and extending measurements of the shear modulus up to 100 kPa [2]. A necessary step in constructing the stiffness map of a tissue from the MRE data is the assumption of a linear viscoelastic model for the tissue. In the simplest case a single spring (Hookean) or a two element parallel spring-dashpot (Voigt) constitutive model are used. In addition to assuming linearity, many techniques have been proposed to solve the inverse problem of estimating the tissue elastic modulus, E, and viscosity, f, from MRE measurements of wavelength and attenuation coefficient. Some techniques use algorithms that make an assumption of tissue incompressibility, and then use a spatial filter to remove displacements caused by compression wave motion. Other MRE reconstruction techniques are based on an algebraic inversion of the differential equation (AIDE) for the medium (e.g., [21]). AIDE techniques directly estimate elastic modulus and attenuation with data smoothing to minimize noise problems. AIDE techniques may also assume incompressibility to reduce complexity and noise susceptibility of the inversion. Other techniques use finite element (FE) analysis iteratively within the material identification algorithm or apply other model-based reconstructions (e.g., [22]). In our work we have to date used the two element Voigt model, but one of the goals of the current work is to employ new multi-scale constitutive models of cartilage to better understand how microscopic MRE measurements can be related to macro-scale biomechanics. 4. Cartilage tissue multi-scale models In order to correctly interpret MRE measurements of cartilage in terms of its micro-scale tissue structure a better understanding of multi-scale tissue mechanics is needed. Recent work [23] has shown that the macro-scale mechanics of cartilage does not follow directly from the mechanics of the individual fibrils, but that the tissue’s strength strongly depends on the arrangement of the fibrils at the micro-scale. Cartilage tissue also shows hysteresis at the macro-scale, which is not evident in isolated fibril mechanics, but rather is indicative of a friction-like sliding motion of the tissue components [1]. Finally, the
660
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
mechanical properties of cartilage are time and history dependent, that is, the observed strain, e(t), will be different when the rate at which stress, r(t), is applied changes and when different stress paths are followed to arrive at the same final load stress. While the dynamics of individual collagen fibrils at the micro-scale may be reasonably modeled using conventional viscoelastic approaches incorporating Hookean (r(t) = E e(t)) and Newtonian fluid elements (r(t) = gde(t)/dt) for elasticity (E) and viscosity (g), respectively (e.g., Voigt, Maxwell, Kelvin models, etc.), extending this to the macro-scale with many fibrils arranged in complex orientations, where the orientation also affects the dynamics, requires some form of spatial averaging and approximation. Considering the problem to be linear (neglecting hysteresis), this is precisely the type of problem for which fractional calculus techniques have shown great promise [24,25]. If, in addition, non-linear hysteretic behavior can be incorporated into the fractional calculus model, then we can develop simplified mathematical expressions, of hysteresis models for cartilage that have a range of applications. Current work in applying fractional calculus to viscoelastic theory starts with the idea that it is the order of the derivative of the strain that – to a first approximation – characterizes the material’s behavior (assuming a one-dimensional stress– strain relation, r(t) = E*dae(t)/dta). In this model, the order of the derivative is zero for a Hookean solid (E* = E) and to one for a Newtonian fluid (E* = g). Viscoelastic materials occupy the intermediate range with a fractional order ‘‘a” between zero and one. Thus, one can build a multi-component fractional equivalent of the ‘‘standard” linear solid (or fluid) by replacing one or more springs and dashpots with ‘‘Springpots”. There are a number of papers (and one monograph) that outline this approach [24,6]. The rationale for the fractional calculus approach is that it encompasses a wider range of observed viscoelastic data (in the time or frequency domain) in a simpler (fewer elements), yet linear model. Many complex, composite and biomaterials exhibit fractional order behavior through power law decays in their creep or relaxation response and their dynamic elastic modulus. Alternatively, quasi-linear elastic and non-linear elastic models has been used to analyze these cases. In our work to ‘‘span” the nano-, micro- and meso-scales we choose to apply the fractional calculus models. Such fractional order constitutive models can also be applied to our elastography experiments. Note in particular, Coimbra and coworkers [26–28], however, go further; they generalize the fractional derivative parameter ‘‘a” so that it becomes a function of time: a(t). Such time varying derivative operators have only recently been introduced (e.g., [29]) and are little studied. There are many issues still to be resolved; but, this model of linear viscoelastic materials provides a bridge between the different fractional calculus models, and it can perhaps also include strain rate and strain path dependence in the formalism. Briefly, Coimbra and coworkers show that the so-called variable order definition of the fractional derivative is simply a properly weighted sum of a large number of fractional order derivative terms. This extends the idea that a single fractional order derivative is the sum of properly weighted integer order spring and dashpot elements [29]. In addition, the variable order fractional derivative model of stress/strain is postulated to be directly related to the rate of change of the long range order of the molecules within a viscoelastic material. In the case of simple compression of a carbon–epoxy composite, Ramirez and Coimbra [28] demonstrate that a(t*) is proportional to an ‘‘entropy” like function t*ln(t*) where t* is a dimensionless time for compression experiments spanning a wide range of strain rates. For example, in this paper, they present a normalized stress versus strain plot as a single curve that extends over eight orders of magnitude in strain rate using only one parameter. Thus, the variable order fractional derivative yields a constitutive model that on the one hand connects with the underlying molecular disorder (disorder that evolves with compression of the composite), and on the other hand with a macroscopic description of mechanical state through an extended viscoelastic model of the physical specimens under study. In our work, we apply these approaches to model complex biomaterials, which exhibit complex viscoelastic behavior. The hypothesis is that fractional and perhaps variable order calculus will provide a greater understanding of the molecular events that occur as cartilage grows, ages and responds to injury. Such understanding can lead to a simpler description of the underlying multi-scale processes that occur when composite materials, such as cartilage, are stressed or strained. Now, returning to our consideration of hysteresis modeling, recently Royston [30] established the relationship of an Iwan hysteresis model consisting of series Coulomb friction and elastic slide components arranged in parallel, to several other models of hysteretic behavior that have been used to simulate a wide range of physical phenomena, most significantly the classical Preisach model. These hysteresis models all possess the necessary and sufficient properties of ‘‘wiping-out” and ‘‘congruency”. Wiping-out means that the current state is only dependent on past extrema; congruency means that the response to a cyclic input is identical in form (forms the same shape hysteresis curve) regardless of previous extrema, which only serve to shift the entire output loop. It was then possible to extend analytical tools developed for the Preisach model to all of the models with which it is related, and vice versa. Such tools include experimental identification, inversion and analysis of dynamic energy flow and dissipation. This prior work provides a framework upon which there can be extensions in several directions relevant to biological tissues. For example, while Royston [30] focused on one-dimensional behavior, coupled multi-dimensional anisotropic behavior could be considered. Given the developed links between the elasto-slide and Preisach model formulations, a logical approach may be to extend techniques developed for multi-dimensional vector-based Preisach models to the anisotropic multi-scale problem [31]. In a more general sense, since hysteresis in and of itself, is actually a manifestation of the presence of multiple scales in a system [32], it is appropriate that it be included as a part of a multi-scale modeling approach to elastic biomaterials.
661
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
Fractional calculus models provide a relatively simple way to describe the biomechanics of complex, porous, or composite systems. There is a multi-scale generalization inherent in the definition of the fractional derivative that accurately represents interactions occurring over wide dynamic ranges of space or time. Thus, there is no need to segment or compartmentalize systems into subsystems or subunits all the way down to the molecular level – a process that often creates more complexity than can be experimentally evaluated. Also, viscoelastic experiments conducted on complex materials can shed light on the ‘‘meaning” of fractional order operations; when the results match the predictions of fractional order system dynamics the physics and the mathematics are in correspondence. Thus, we can begin to unravel the contributions to the physics that follow from the fundamental model dynamics, the geometry and the interaction between system components and physical barriers. In the following, we will present some initial steps taken to apply this methodology to the viscoelastic behavior of cartilage using MRI and MRE. 5. Biological tissue measurements MRI examination of tissue-engineered cartilage provides a non-invasive way to visualize changes in structure while MRE measures the inherent elastic properties (complex modulus) of the developing tissue. The objective of this ongoing investigation is to determine the key MR parameters that describe the essential structural characteristics of new cartilage tissue derived from engineered human mesenchymal stem cells (hMSC). The hMSCs used in this study were isolated from commercially available bone marrow, expanded in cell culture and seeded in porous poly(ethylene glycol) scaffold. The tissue engineered construct was exposed to chondrogenic induction media for 4 weeks and samples evaluated each week (Fig. 2). All samples were sealed in glass NMR tubes filled with media and loaded into a 10-mm diameter RF saddle coil for the MRI studies. T2 maps were calculated each week using a spin echo imaging sequence to acquire 32 echoes with a 7 ms echo spacing for the chosen axial slice. The data were analyzed to produce a pixel-by-pixel T2 map for the entire sample scope (Fig. 2). As the cartilage tissue developed contrast appeared throughout the slice and a monotonic increase in T2 was observed. To confirm that the T2 changes reflect progression of the hMSC along the developmental pathway to form cartilage, samples were removed for histochemical analysis: safranin O and alcian blue staining, collagen Type II immunohistochemistry, and quantitative glycosaminoglycan (GAG) analyses. The constructs stained positive for collagen Type II and GAGs as time increased from week 0 to week 4. Furthermore, the quantitative GAG assay established a larger increase of GAG in the experimental constructs relative to the controls [33]. Other recent studies of tissue-engineered cartilage examined other MR-based parameters: T1 relaxation in the rotating frame (T1q) and magnetization transfer (MT). T1q and MT are magnetic resonance derived parameters known to be sensitive to tissue composition and structure. In brief, both reflect the extent of water binding to proteins and other macromolecules. Our hypothesis for this study is that developing tissue will display differences in the relative amount of free and ‘‘so-called” bound water. The bound protein fraction or BPF can be derived from the MT data. Previous studies of bone and adipose tissues showed that the tissue water relaxation times (T1 and T2) do not always show early changes in tissue composition; changes that need to be measured in order for tissue engineers to optimize tissue growth. The bar graph shown in Fig. 3
Cartilage
MRI
Safranin-O
Superficial zone Superficial zone
Transitional zone
Chondrocytes
Radial zone
Proteoglycan
Transitional zone Radial zone Calcified zone
Calcified zone
Bone Marrow
140 120 100 80 60 40
Time (wk)
10 mm
0
1
2
3
4
Fig. 2. (Top) The three images of cartilage are: an anatomical drawing, a high resolution MRI, and a histological image (safranin-O/fast green-stained section (40)). In the MR image, three zones can be easily identified. MR acquisition parameters were: TR/TE = 1000/ 60 ms; field-of-view (FOV) = 0.6 cm; slice thickness = 0.5 mm; in-plane resolution = 22.7 lm 22.7 lm; and 64 averages [3]. (Bottom) Five T2 maps for engineered cartilage over a 4 week growth period. T2 contrast and a monotonic increase of T2 is observed with time.
662
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664 180 Control
160
Stimulated
140
T1p in milli-sec
120 100 80 60 40 20 0 Day 0
Day 3
Day 6
Day 9
Day 12
Fig. 3. T1q values derived from T1q weighted images of tissue-engineered cartilage. The data are means ± one standard deviation for four measurements.
shows the significant changes in the T1q parameter over the duration of this study. A similar change was also observed in BPF (see Fig. 4). In a tissue engineering study by Marion and coworkers [5] the growth stages of engineered fat were monitored using MR measurements of relaxation times and the water diffusion coefficient. In this study of MSC derived adipogenic cells the T2 values stayed unchanged for 4 weeks, but MR elastography measurements of the bulk modulus showed a significant decrease (Fig. 5). By week 4, for example, the shear modulus had decreased to less than one fourth of the initial value (5.8 to 1.4 kPa). The decrease in shear stiffness indicates that the engineered adipogenic tissue becomes ‘softer’ as adipogenic matrix production gradually increases. These results illustrate the complementary role of MRI and MR elastography in assessing the development of engineered tissue in vitro. 6. Combined hysteresis and fractional calculus models In the generalized elasticity model of tendons and ligaments [1] a number, say N, of springs are connected in parallel. An applied force induces a stress that as it grows engages a sequence of spring elements beginning with the weakest (lowest elastic constant En) extending first and the strongest (largest elastic constant EN) last. The resulting stress–strain behavior is not a straight line, but a curve with an initial flat section and an increasing slope. We propose to extend this idea by using a parallel combination of N individual Coulomb sliding friction element and one or more fractional order springpots (Fig. 5). The sliding friction Iwan parallel-series hysteresis model captures the rate-independent hysteresis behavior, while the springpot captures the power law dynamics observed in tissue mechanics.
Fig. 4. T2 relaxation time and shear stiffness data for engineered fat.
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
663
μh φ1
φ2
σ[t] μh
φm φN
μh εN[t,φN] ε[t]
Fig. 5. Combined hysteresis and fractional model.
Under quasistatic loading conditions, the springpot asymptotically carries no load and hysteretic relation between stress
r and strain e is captured via the elasto-slide elements [30]. Under dynamic loading conditions, elasto-slide element behavior is unchanged. However, now the springpot, acting in parallel, introduces a rate-dependent stress/strain component. In the time domain, the stress–strain relation for the hysteretic component can be expressed as follows. Consider a finite number N of parallel-series elasto-slide elements. Referring to Fig. 5 and separating the stress related to hysteresis rH[t] from that of the dashpot rF[t], such that r[t] = rH[t] + rF[t], we have: N X
lh ðe em ðt; /m ÞÞ ; lh /m sign½e_ m¼1 je em ðt; /m Þj 6 /m unchanged em ½t; /m ¼ : if je em ðt; /m Þj ¼ /m e /m sign½e_ rH ½t ¼
rm ½t; rm ½t ¼
ð1Þ
In the time domain, the stress–strain relation for the fractional component, or springpot, has dynamics that can be expressed as, a
rF ðtÞ ¼ Esa
d eðtÞ a ; dt
ð2Þ
where a can take values between 0 and 1. Here a = 0 represents a linear spring and a = 1 represents a damper (Es = g). Any specific value of a can be represented with hierarchical arrangements of springs and dashpots, such as ladders, trees and fractal networks [6]. 7. Conclusions These new multiple scale techniques extend the existing computational methods for predicting stiffness and strength, assessing short versus long term load response, and measuring static versus dynamic response to mechanical loads over a wide range of frequencies. While the specific application considered in this project involves tissue-engineered cartilage, the multi-scale computational methodology used to connect nano-scale structure with macro-scale function is relevant to a much wider range of biomechanical systems and biomedical problems. Our new methods could, for example, help identify early microscopic changes in connective tissue at the microscopic scale. Knowledge of such changes would enable early diagnostic monitoring for the onset of disease and better assessment of the effectiveness of new drugs or therapies. References [1] Ethier CR, Simmons CA. Introductory biomechanics: from cells to organisms. Cambridge, UK: Cambridge University Press; 2007. [2] Othman SF, Xu H, Royston TJ, Magin RL. Microscopic magnetic resonance elastography (microMRE). Magn Reson Med 2005;54:605–14. [3] Othman SF, Li J, Abdullah O, Moinness JJ, Magin RL, Muehleman C. High-resolution/high-contrast MRI of human articular cartilage lesions. Acta Orthop 2007;78:536–46. [4] Xu H, Othman SF, Liu H, Peptan IA, Magin RL. Magnetic resonance microscopy for monitoring osteogenesis of tissue-engineered constructs in vitro. Phys Med Biol 2006;51:719–32. [5] Marion N, Othman SF, Stosich MS, Magin RL, Mao JJ. Measurement of tissue engineered adipogenic constructs using magnetic resonance elastography. Abstr Annu Meet BMES 2006. [6] Magin RL. Fractional calculus in bioengineering. Redding, CT: Begell House; 2006. [7] Hilfer R. Fractional calculus in bioengineering. River Edge, NJ: World Scientific; 2000. [8] Grimnes S, Martinsen OG. Bioimpedance and bioelectricity basics. San Diego, CA: Academic Press; 2000. [9] Lakes RS. Viscoelastic solids. Boca Raton, FL: CRC Press; 1999. [10] Bard AJ, Faulkner LR. Electrochemical methods: fundamentals and applications. second ed. New York, NY: John Wiley & Sons; 2001. [11] Meachim G, Stockwell RA, The matrix, in: M.A.R. Freeman (Ed.), Adult articular cartilage, Pitman Medical, Kent; 1979. p. 1–50. [12] Hunziker EB. Articular cartilage repair: basic science and clinical progress, a review of the current status and prospects. Osteoarthr Cartil 2002;10:432–63. [13] Pietrzak WS. Musculoskeletal tissue regeneration: biological materials and methods. Totowa, NJ: Humana Press; 2008.
664
R.L. Magin, T.J. Royston / Commun Nonlinear Sci Numer Simulat 15 (2010) 657–664
[14] Speer DP, Chvapil M, Volz RG, Holmes MD. Enhancement of healing in osteochondral defects by collagen sponge implants. Clin Orthop Relat Res 1979;144:326–35. [15] Elisseeff JH, Puleo C, Yang F, Sharma B. Advances in skeletal tissue engineering with hydrogels. Orthod Craniofac Res 2005;8:150–61. [16] Temenoff JS, Mikos AG. Review: tissue engineering for regeneration of articular cartilage. Biomaterials 2000;21:431–40. [17] Mao JJ. Stem-cell-driven regeneration of synovial joints. Biol Cell 2005;97:289–301. [18] Caplan AI. Tissue engineering designs for the future: new logics, old molecules. Tissue Eng 2000;6:1–8. [19] Sharma B, Elisseeff JH. Engineering structurally organized cartilage and bone tissues. Ann Biomed Eng 2004;32:148–59. [20] Muthupillai R, Lomas DJ, Rossman PJ, Greenleaf JF, Manduca A, Ehman RL. Magnetic resonance elastography by direct visualization of propagating acoustic strain waves. Science 1995;269:1854–7. [21] Oliphant TE, Manduca A, Ehman RL, Greenleaf JF. Complex-valued stiffness reconstruction for magnetic resonance elastography by algebraic inversion of the differential equation. Magn Reson Med 2001;45:299–310. [22] Van Houten EEW, Miga MI, Weaver JB, Kennedy FE, Paulsen KD. Three-dimensional subzone-based reconstruction algorithm for MR elastography. Magn Reson Med 2001;45:827–37. [23] Chandran PL, Barocas VH. Affine versus non-affine fibril kinematics in collagen networks: theoretical studies of network behaviour. ASME J Biomech Eng 2006;128:259–70. [24] Magin RL. Fractional calculus in bioengineering. Crit Rev Biomed Eng 2004;32(1–4):1–377. [25] Sinkus RK, Siegmann T, Xydeas M, Tanter C, Claussen C, Fink M. MR elastography of breast lesions: understanding the solid/liquid duality can improve the specificity of contrast-enhanced MR mammography. Magn Reson Med 2007;58:1135–44. [26] Coimbra CFM. Mechanics with variable-order differential operators. Ann Phys 2003;12:692–703. [27] Soon CM, Coimbra CFM, Kobayashi MH. The variable elasticity operator. Ann Phys 2005;14:378–89. [28] Ramirez LES, Coimbra CFM, A constitutive equation for linear viscoelastic thermoset materials undergoing compression, In: Proceedings of the of ASME IDETC/CIE conference. Paper No. DETC2007-34169; 2007. [29] Lorenzo CF, Hartley TT. Variable order and distributed order fractional operators. Nonlinear Dyn 2002;29:57–98. [30] Royston TJ. Leveraging the equivalence of hysteresis models from different fields for analysis and numerical simulation of jointed structures. ASME J Comput Nonlinear Dyn 2008;3:031006-1–6-8. [31] Mayergoyz ID. Mathematical models of hysteresis and their applications. NY: Elsevier, Academic Press; 2003. [32] Banks HT, Luke NS, Samuels Jr, JR, Viscoelasticity in polymers: phenomenological to molecular mathematical modelling. Numer Methods Partial Differ Eqn. doi:10.1002/num.20250; 2008. [33] Troken AJ, Abdullah OM, Othman SF, Magin RL, Mao JJ. Magnetic resonance imaging for monitoring chondrogenesis of stem-cell based constructs in vitro. Abstr Annu Meet BME 2006.