Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Inhomogeneous Response of Articular Cartilage: A Three-Dimensional Multiphasic Heterogeneous Study

  • Sara Manzano,

    Affiliations Mechanical Engineering Department, School of Engineering and Architecture (EINA), University of Zaragoza, Spain, Aragón Institute of Engineering Research (I3A), University of Zaragoza, Spain, Biomedical Research Networking Center in Bioengineering, Biomaterials and Nanomedicine (CIBER-BBN), Spain

  • Monica Armengol,

    Affiliations Department of Engineering Science, University of Oxford, Oxford, United Kingdom, Nuffield Department of Orthopaedics, Rheumatology and Musculoskeletal Sciences, University of Oxford, Nuffield Orthopaedic Centre, Oxford, United Kingdom

  • Andrew J. Price,

    Affiliation Nuffield Department of Orthopaedics, Rheumatology and Musculoskeletal Sciences, University of Oxford, Nuffield Orthopaedic Centre, Oxford, United Kingdom

  • Philippa A. Hulley,

    Affiliation Nuffield Department of Orthopaedics, Rheumatology and Musculoskeletal Sciences, University of Oxford, Nuffield Orthopaedic Centre, Oxford, United Kingdom

  • Harinderjit S. Gill,

    Affiliations Nuffield Department of Orthopaedics, Rheumatology and Musculoskeletal Sciences, University of Oxford, Nuffield Orthopaedic Centre, Oxford, United Kingdom, Department of Mechanical Engineering, University of Bath, Bath, United Kingdom

  • Manuel Doblaré,

    Affiliations Mechanical Engineering Department, School of Engineering and Architecture (EINA), University of Zaragoza, Spain, Aragón Institute of Engineering Research (I3A), University of Zaragoza, Spain, Biomedical Research Networking Center in Bioengineering, Biomaterials and Nanomedicine (CIBER-BBN), Spain

  • Mohamed Hamdy Doweidar

    mohamed@unizar.es

    Affiliations Mechanical Engineering Department, School of Engineering and Architecture (EINA), University of Zaragoza, Spain, Aragón Institute of Engineering Research (I3A), University of Zaragoza, Spain, Biomedical Research Networking Center in Bioengineering, Biomaterials and Nanomedicine (CIBER-BBN), Spain

Abstract

Articular cartilage exhibits complex mechano-electrochemical behaviour due to its anisotropy, inhomogeneity and material non-linearity. In this work, the thickness and radial dependence of cartilage properties are incorporated into a 3D mechano-electrochemical model to explore the relevance of heterogeneity in the behaviour of the tissue. The model considers four essential phenomena: (i) osmotic pressure, (ii) convective and diffusive processes, (iii) chemical expansion and (iv) three-dimensional through-the-thickness heterogeneity of the tissue. The need to consider heterogeneity in computational simulations of cartilage behaviour and in manufacturing biomaterials mimicking this tissue is discussed. To this end, healthy tibial plateaus from pigs were mechanically and biochemically tested in-vitro. Heterogeneous properties were included in the mechano-electrochemical computational model to simulate tissue swelling. The simulation results demonstrated that swelling of the heterogeneous samples was significantly lower than swelling under homogeneous and isotropic conditions. Furthermore, there was a significant reduction in the flux of water and ions in the former samples. In conclusion, the computational model presented here can be considered as a valuable tool for predicting how the variation of cartilage properties affects its behaviour, opening up possibilities for exploring the requirements of cartilage-mimicking biomaterials for tissue engineering. Besides, the model also allows the establishment of behavioural patterns of swelling and of water and ion fluxes in articular cartilage.

Introduction

Articular cartilage is a highly heterogeneous, anisotropic and multiphasic material mainly composed of a solid phase (collagen fibrils and glycosaminoglycans (GAGs)), a fluid phase (interstitial liquid) and an ionic phase (Na+, Cl-, Ca+ and K+ among others) [13]. This composition and the associated structure determine the heterogeneous mechano-electrochemical properties and the inhomogeneous behaviour of the tissue [4, 5]. The main components involved in the swelling process of this tissue are a dense collagen network with attached GAGs, ions and the interstitial liquid that transports them, and oxygen and nutrients for tissue maintenance [6, 7]. The difference in ion concentration between the extra-cellular matrix and the interstitial liquid induces an osmotic pressure variation, also known as Donnan osmotic pressure. From the biological point of view, swelling is an essential process for the exchange of essential nutrients, oxygen, carbon dioxide and waste products in articular cartilage [8]. It has been reported that swelling is a good indicator of tissue quality [911]. When cartilage is degraded (for instance in advanced grades of osteoarthritis, long-term immobilized patients and severe cartilage degenerative diseases), collagen fibers and proteoglycans (PGs) lose their ability to retain water molecules. As a result, the physiological hydration of the cartilage increases [1217]. In contrast, in the case of short-term immobilized patients and early osteoarthritis, the matrix becomes stiffer, exerting higher resistance to fluid passage [12, 13, 18, 19]. Despite intensive research carried out on cartilage-like biomaterials, the success of new biomaterials in mimicking cartilage behaviour remains elusive. No effective constructs for chondrocyte growth have been developed and advances in cartilage-mimicking biomaterials are limited [20]. The main reason seems to be directly related to the lack of consensus on simple answers to basic questions such as: what is the role of the inhomogeneity of the tissue in swelling? How do variations in tissue properties affect ion fluxes? Does neglecting these variations, when designing cartilage-mimicking materials, partly account for their lack of success? Three-dimensional (3D) computational models can provide additional insight into how the variation in the biomaterial composition affect its behaviour when implanted in the human body. Numerical modelling allows this to be done in an accurate and time efficient manner [8, 21, 22]. In the last decade, the majority of works related to articular cartilage simulations where developed in 2D. However, in the last years, advances in this field toward more accurate 3D simulations have achieved, which leads to improvement in the studying of specific cartilage phenomena as swelling [23, 24], neglecting important phenomena which occur under more realistic 3D environmental conditions. Moreover, most currently reported models consider cartilage as a homogeneous material, so the spatial variation in its properties is neglected [18]. More models that are recently incorporate transverse anisotropy, considering different cartilage layers, while neglecting the existing radial variations [6, 7]. In this paper, we present a novel 3D mechano-electrochemical model (based on our previous model [12, 13]) which aims to address these important questions and to assist in the analysis of the effects of articular cartilage anisotropy and inhomogeneity on tissue behaviour. In this respect, we present a novel 3D model that includes all the essential features of the articular cartilage, namely: (i) Donnan osmotic pressure, (ii) convective and diffusive processes, (iii) variations of chemical expansion due to the heterogeneous distribution of PGs and (iv) three-dimension through-the-thickness heterogeneity of the tissue. To achieve these goals, several experiments were conducted to measure variations in the properties along the cartilage surface and through the thickness. Specifically, indentation tests were carried out to obtain Young´s modulus (Es) and histological analysis of tissue samples was undertaken to qualitatively establish the GAG distribution in six tibial plateaus from porcine specimens. These properties were included in the model to analyse their influence on free swelling, on water and cation fluxes and ion distribution. A comparison between the obtained results and the results of the homogeneous behaviour were performed. They showed significant lower swelling in stiffer areas. Indeed, swelling was altered by the lower water and ion fluxes that appear in such specific areas of samples. Furthermore, the cation distribution in the heterogeneous configuration has more significant alterations compared to that obtained with homogeneous simulations. In the former case, the cation distribution is significantly reduced in magnitude.

Material and Methods

Detailed experiments were conducted to measure Es and the qualitative GAG distribution in the three-dimensional environment of the tissue for six porcine tibial plateaus obtained commercially from an abattoir (W.H.Alder in Cowley Road, Oxford). A Whole Articular Surface Indentation Machine (WASIM), specially designed to scan and indent articular surfaces, was used to measure the Es. Qualitative GAG distributions were obtained from histological assays of samples extracted from these tibial plateaus after the mechanical testing had been completed. Note that specific values for GAG concentration within the sample and the subsequent initial fixed charge density () are included in the computational model were estimated following the method of Lai et al. [25] as well as our own histological results. This means that qualitative GAG distribution and the subsequent distribution in cartilage zones are included in our model.

The WASIM was specifically designed to measure mechanical properties of intact specimens. Due to its five degrees of freedom, in-situ mechanical properties can be obtained over the tibial surface using indentation normal to the surface at any given point on the articular surface.

The WASIM consists of a stage that can translate in the X and Y directions and rotate in pitch. The sample holder located in the middle of the stage can rotate about its central longitudinal axis and thus provide yaw rotations. The indenter is mounted so that it can translate in the Z direction, but is fixed in the X and Y directions. Indentation tests were carried out using a hemispherical indenter (1.35 mm radius). The WASIM also incorporates a high resolution laser (LK-G32, Keyence, Higashi-Nakajima, Japan) for measuring the height of the specimen surface (i.e. measurement in the Z direction).

The high resolution laser combined with raster scanning along the X and Y axes provides a high-accuracy topographical map of the tibial surface. Using this topographical map, test points over the tibial plateau were selected and normal vectors of test points were calculated. Using these normal vectors, the WASIM rotation angles and translations were calculated to achieve normal indentation at the selected test points. Es was obtained using the Field and Swann method and was adjusted by cartilage thickness measured post-hoc:

Where ν is the Poisson’s modulus (set to 0.4), P is the force, h is the penetration and dP/dh corresponds to the unloading curve of the force-displacement curve. Ac is the area of compression given by πa2 and was calculated from Filed and Swain method [26, 27].

The solid phase heterogeneity is an important aspect to be included in accurate cartilage simulation. This was indirectly considered in this model through the heterogeneous Es and distribution. Both parameters vary though the thickness and in radial direction due to the different collagen fibers orientation and distribution. In this sense, experimental measures reveal this collagen distribution providing different rigidity in both directions.

Sample description

Experimental tests were performed on six porcine tibial plateaux. For this reason, pig legs were purchased from a commercial abbatoir (W.H.Alder in Cowley Road, Oxford) following routine humane slaughter. No approval process or legislation governs use of commercial butcher’s meat for experimental use post slaughter. The porcine specimens were males aged between 6 month and 1 year. Porcine specimens were selected due to their greater similarities with human tibial plateau in terms of size and shape compared with ovine or equine species [28]. The porcine knee joints were dislocated and dissected to obtain intact tibial plateau with subchondral bone and tibial shank.

Cartilage sample preparation

Soft tissues, ligaments and menisci were removed from tibial plateau samples. They were then fixed onto the sample holder of the WASIM using bone cement (VersoCit-2 kit, Struers, UK). The cartilage was kept hydrated throughout the indentation process, by regularly applying saline solution (0.15 M) to the sample. After mechanical tests were performed, four millimetre diameter cartilage/subchrondral bone plugs, centred on previous indentation points, were extracted from the tibial plateau using an osteotome.

Indentation test

Displacement controlled indentation was carried out at 16 points of measurement in the upper surface of the cartilage samples up to 10% maximum strain (Figs 1 and 2). Thickness of the sample was measured after mechanical tests were performed, therefore an adjustment due to articular cartilage thickness was carried out. Displacement controlled indentation was applied at a rate of 10 percent per second (pps). This rate allows the assumption cartilage behaviour as an elastic material at each test point. Note that in this experiment we consider that the majority of water content placed in the superficial zone of the samples is exuded at the very beginning of the process. Thus, the resulting measurement of mechanical properties refer to a superficial zone conformed by collapsed collagen fibrils, assuming almost null water content.

thumbnail
Fig 1. Intact porcine tibial plateau sample used for indentation and histological tests.

The red circle indicates the area of interest where indentation tests were performed.

https://doi.org/10.1371/journal.pone.0157967.g001

thumbnail
Fig 2. Axial view of the topographical map of the joint surface and the specific points selected for indentation.

https://doi.org/10.1371/journal.pone.0157967.g002

To estimate the mechanical properties in sample thickness (Edepth and νdepth), the following equations were used [28], (1) (2) where Esup and νsup are the Young´s modulus and Poisson ratio at the articular cartilage surface, respectively; h is the thickness of the sample; z denotes the depth of the layer considered and αE and αν represent specific material constants (αE ≥ 0; αν ≥ 0). For this simulation, αE = 3.2 and αν = 2.8 [29]. Indentation tests were performed to measure top layer properties (Esup and νsup) Then, these superficial values were applied in Eqs (1) and (2) to estimate mechanical properties of the sample through the thickness. For the homogeneous simulation, E was obtained from a long time experimental indentation, which provides this parameter considering the sample as homogeneous material. In contrast, for heterogeneous simulations, this value was extracted for each discrete region of cartilage top layer and then mathematically extrapolate to get values in depth.

Histological analysis

Samples were prepared for histological analysis to observe cross-sectional GAG distribution. The cartilage/bone plugs were immersed in 10% neutral buffered formalin and 4% EDTA for 72 hours to decalcify. Samples were then cryofrozen and sectioned into 5 μm thick sections for staining with Safranin O. It has been shown that these techniques of fixation and decalcification has a small effect in GAG content [30, 31]. Finally, the samples were observed under 40x and 20x magnification using a light Olympus BX40F4 microscope. Histological images were recorded using an Olympus U-READDB 60X camera. Colour gradients were obtained indicating the degree of GAG concentration within the samples. Based on this obtained qualitative distribution, the specific values for the concentrations of the fixed charges attached to the GAGs were obtained from Lai et al. [25].

Computational formulation of the Mechano-electrochemical model

The main model formulation is explained below, building further on our previous work [12, 13]. The model was based on the triphasic theory for soft, charged and hydrated tissues [21]. Four specific components were considered: negatively charged porous-elastic solid (s), fluid (f), cations (+) and anions (-). The interaction between these phases results in the mechano-electrochemical phenomena required for cartilage maintenance (further information about the model formulation may be found in [7, 12, 13, 32]). A summary of the main model formulation is presented below.

Governing equations of the computational model of cartilage behaviour.

It is important to note that the model has been validated and established in previous works [12, 13, 14]. Considering the four basic unknowns (us the displacement of the solid matrix, εw the chemical potential of water, ε+ and ε the electro-chemical potential for cations and anions, respectively), the four constitutive equations were established as:

  1. Momentum balance equation of the mixture (3)
  2. Mass balance equation of the mixture (4)
  3. Charge balance equations of ions (5) (6)

In these equations, σ refers to the total mixture stress tensor whereas σf, σc and σs are the stress tensors related to the fluid, chemical and solid phases, respectively. The velocity of the fluid matrix is defined as . Regarding the ion concentrations, c+ and c denote cation and anion concentrations, respectively. Finally, Φw represents the porosity of the tissue [12, 13, 14]. The water flux, Jw, cation flux, J+, and anion flux, J, were expressed as a combination of the electrochemical potentials: (7) (8) (9) Here, α refers to the drag coefficient between the solid and the water phases. The cation and anion diffusivities are represented by D+ and D, respectively. The experimental atmospheric conditions were included in the model by R, the universal gas constant, and the absolute temperature, T [7]. Note that JW refers to advection of the species with the moving pore fluid. However, J+/− are the diffusion within and advection with the fluid flowing relative to the solid matrix.

The state variables appearing in Eqs 36 in terms of the basic variables of the problem may be also written as follows: (10) (11) (12) (13) Where Fc is the Faraday constant, Ψ the electrical potential, Bw the fluid-solid coupling coefficient, Φ the osmotic coefficient, γ+ and γ the activity coefficient of the anions and cations, respectively and I the identity tensor. P is the fluid pressure, θ = div[us] the solid matrix dilatation relating to the infinitesimal strain tensor of the solid matrix, λS and μS the Lame constants, and ϵ the infinitesimal strain tensor of the solid matrix [12, 13]. To consider the chemical expansion, TC, a new term was included in the mathematical formulation. This expansion was due to the presence of charge-to-charge electrostatic repulsive forces exerted on the GAG-collagen network (solid phase). Note that the equilibrium of cartilage swelling depends on the combined action of the chemical expansion stress, σC, and the osmotic pressure. Both are related to the internal concentration and distribution of ions and their interaction with each other. Thus, TC is represented as follows: (14) Where a0 and k are the PG repulsion coefficients, and cF the fixed charge density attached to PGs. γ± and γ±* are the mean activity coefficients of ions during the process and in the reference state, respectively. c refers to the neutral external salt concentration [7].

Discretization.

The tissue response in confined condition was solved using the finite element method. The primary unknowns of the model [u,εw,ε+,ε] were interpolated from nodal values. The time derivatives were approximated with the Crank-Nicolson method which yields an implicit approximation to the solution of the initial value problem y’ = f(x,y) with y(x0) = y0 at x for a given time step h [29]. To obtain the fully coupled non-linear system of equations describing the discretized model, we first established the weak formulation of the governing equations: (15) (16) (17) (18) where n is the unit vector normal to the boundary and δu, δε and are the so-called test functions. The superscript * stands for the quantities in the bathing solution and ck = c+ + c Eqs 1417 were expressed in terms of primary variables in our previous works related to this model formulation and validation [12, 13, 14]

Numerical implementation.

Tri-linear 8-noded hexahedral elements with 2×2×2 Gaussian integration points were used. The selected mesh resulted in 1680 elements. The finite element formulation described above was implemented in a user-defined element subroutine of the commercial software package Abaqus 6.11 (Dassault Systemes, Paris). To study the effect of through-the-thickness heterogeneity in cartilage behaviour, the experimental design described by Lai et al. [25] was computationally reproduced (Fig 3). Thus, a cartilage specimen of 30 mm diameter and 3 mm depth placed inside a circular confining ring was considered. The concentration of the external bath solution was decreased from 0.15 M to 0.125 M to mimic cartilage swelling observed in physiological conditions. No loads were applied on the sample in the z-direction. The properties of healthy homogeneous isotropic cartilage are listed in Table 1 while the properties of the through-the-thickness heterogeneous cases are listed in Tables 2 and 3. A comparison between the simulated results of the homogenous and the heterogeneous cases were performed. Quantification of cation fluxes and distributions within the samples and monitoring the tissue changes during the swelling processes were also achieved. Note that the homogeneous model was extracted from the literature [28] since this value was directly measured in the same kind of samples (porcine tibial plateaus) through a specific experimental procedure that provides the measure of the samples considering them as homogeneous as a whole. In this manner, side errors related to perform a mathematical approximation (medium of the measures of in layers) were avoided.

thumbnail
Fig 3. Schematic representation of the transient free swelling experiment (bottom right) and the boundary conditions of the cartilage sample applied to the computational simulation.

A cartilage sample of 3 mm thickness and 30 mm diameter is immersed in a NaCl solution with an initial concentration of 0.15 M. The sample is confined in an impermeable chamber.

https://doi.org/10.1371/journal.pone.0157967.g003

thumbnail
Table 1. 3D Model parameters for healthy cartilage considered as homogeneous material.

https://doi.org/10.1371/journal.pone.0157967.t001

thumbnail
Table 2. Initial permeability (k) introduced in the 3D computational model for healthy cartilage considered as through-the-thickness heterogeneous material.

https://doi.org/10.1371/journal.pone.0157967.t002

thumbnail
Table 3. Initial fixed charged density () introduced in the 3D computational model for healthy cartilage considered as heterogeneous material based on the work of Lai et al. [25] as well as following our GAG-histological results.

Property distribution similar to Young´s modulus represented in the zone of the sample shown in Fig 5. To clarify Table 3–Fig 5 correspondence, note that, for instance, first value in the first line and first column of table 3, 0.150, corresponds to the superficial layer and the point 13 of Fig 5. Values collected in table 3 are obtained in discrete material regions.

https://doi.org/10.1371/journal.pone.0157967.t003

Initial conditions.

Initially, the cartilage sample was equilibrated with a single salt (NaCl) solution with concentration c*. The initial conditions for the computational model were: In this study, the free-swelling state of tissue equilibrated with the bathing solution was chosen as the reference configuration for strain (time 0 seconds, undeformed configuration).

Boundary conditions.

The boundary conditions applied in the confined conditions are summarized below (Fig 3):

  1. Free surface:
  2. Lateral surface:
  3. Lower surface:

At t = 0 seconds the concentration of the external solution, c*, was increased from 0.15 M to 0.125 M. The transient response of the solid displacement was solved by using the new 3D model and compared to simulation results where heterogeneity and anisotropy of the cartilage are neglected.

Results and Discussion

The 3D computational model described above was used to simulate healthy articular cartilage behaviour in an external bath solution in confined conditions considering tissue through-the-thickness properties variation. First, an indentation test and histological analysis were performed on whole undamaged tibial plateau samples extracted from pigs to measure Es along the cartilage surface. The GAG distribution was qualitatively and quantitatively assessed in the same locations as Es Second, both parameters together with those extracted from the literature were introduced into the 3D mechano-electrochemical model to simulate tissue swelling. The results of the cartilage considered as heterogeneous were compared with those obtained in homogeneous and isotropic conditions (Fig 4).

thumbnail
Fig 4. Schematic representation of the experimental tests carried out to obtain the Young´s modulus and GAG distribution (phase 1) for the 3D computational simulation (phase 2) of articular cartilage behaviour.

Note that four random points were chosen to measure (qualitatively) the GAGs distribution.

https://doi.org/10.1371/journal.pone.0157967.g004

From these results, we conclude that the use of computational models allows the swelling behaviour in cartilage to be studied in a non-invasive and inexpensive manner. In addition, the need to incorporate tissue inhomogeneities in the manufacturing process of cartilage-mimicking biomaterial can be guessed in line with the obtained results for cartilage where heterogeneity highly affects tissue simulations [30, 31, 33]. Similarly, this model provides healthy cartilage behavioural patterns that allow specific ranges of swelling to be established as a measure of tissue health.

Due to the lack of data for cartilage pig knee, in homogeneous case was extracted from literature where the same kind of sample were used and they are considered as homogeneous material in its whole. Thus, mathematical errors, derived from medium of values measured in layers, are avoided. On the other hand, in homogeneous case is slightly higher than maximum value of heterogeneous case, this suppose higher swelling of the sample. Importantly, differences between homogeneous and heterogeneous swelling, ion and water cases are observed providing essential information about how important is considering heterogeneity in cartilage simulation. Importantly, a model limitation could be found in the measurement of the mechanical properties of the samples superficial zone. This zone is considered to be solid matrix due to the preliminary exudation of water at the very beginning of the experiment. Other authors [34] offer different manner to perform this mechanical properties quantification. Similarly, it is important to remark that in the case of the homogeneous study, average fixed charged density has been extracted from literature, this can also affect in accuracy when comparing with heterogeneous case.

Indentation tests

The sixteen Es measured on the whole cartilage surface are shown in Fig 5, where the blue columns correspond to the superficial zone. The resulting variation of E through the thickness was calculated by Eq 1 (Figs 5 and 6, green columns corresponding to the middle cartilage zone and red columns to the deep zone). Coincident with observations reported in the literature, zones located in the latero-posterior parts exhibited slightly lower values of E while this parameter was increased in deeper layers of the cartilage [29, 35, 36].

thumbnail
Fig 5. Young´s modulus experimentally obtained from the indentation test for the superficial layer (blue columns), and mathematically obtained for the middle zone (green columns) and deep zone (red columns) [19].

https://doi.org/10.1371/journal.pone.0157967.g005

thumbnail
Fig 6. Panoramic view of a longitudinal section of the cartilage sample extracted from porcine tibial plateaus.

Light microscopy of cartilage and subchondral bone, Safranin O staining, 100X.

https://doi.org/10.1371/journal.pone.0157967.g006

Histological analysis

Histological sections stained with Safranin O were observed by light microscopy. Five zones were clearly identified according to the red colour grading. In the superficial zone, the grey tone indicates low GAG content whereas in deeper areas a high GAG concentration was detected. These findings corroborate experimental observations reported in the literature [25]. The resulting approximate pattern of the GAG arrangement and the subsequent fixed charge distribution in the superficial, middle and deep layers of the articular cartilage are shown in Fig 6. Specific fixed charge values are listed in Table 3, obtained experimentally from the histological GAG distribution. Values collected in Table 3 correspond to discrete material regions.

Heterogeneous cartilage free swelling

To mimic physiological cartilage free swelling, the external bath concentration was decreased from 0.15 M to 0.125 M. Under this condition, the 3D computational model exhibited three distinct phases which were also observed in our previous works: I) an initial sample shrinking; II) a massive entrance of water into the sample during the first 1650 seconds of the free swelling simulation and III) a saturation of the sample, reaching the equilibrium state.

Z-displacement.

The physiological properties of healthy knee cartilage gave rise to a maximal displacement of 0.426 mm and 0.187 mm for the homogeneous and heterogeneous cases, respectively, within 1650 seconds of simulated time; subsequently in both simulations these values slightly decreased due to the minimum outgoing flux of cations, reaching a new equilibrium state at 2500 seconds at the end of simulation. This final state corresponds to 0.284 mm of z-displacement for the homogeneous case and 0.151 mm and 0.076 mm for maximal and minimal z-displacement, respectively, for the heterogeneous case (Fig 7). It is of importance to note in the latter case the substantially lower range of the z-displacement. Consistent with the behaviour of composite materials, when their properties have a heterogeneous distribution, the swelling of these materials is limited. Hence, irregular E and GAG distributions within the sample generate lower z-displacements. This finding has special relevance for biomaterials designed for the replacement of cartilage-damaged areas.

thumbnail
Fig 7. Comparison of the upper surface displacement in the free swelling test from the model considering articular cartilage as homogeneous material (green line) and through-the-thickness heterogeneous material (range marked by blue line, maximum value, and red line, minimum value).

The dark area represents the saturation phase.

https://doi.org/10.1371/journal.pone.0157967.g007

Cation distribution.

Under these conditions, the cation gradient was also monitored after 1650 seconds of free swelling, coincident with the maximum displacement time (Fig 8). The results are again compared with those obtained for the cartilage considered as homogeneous material. The latter case leads to a reduction in the cation concentration towards the surface of the sample, from 255 mol/m3 at z = 0 mm to 160 mol/m3 in the upper surface of the sample (z = 3 mm). In contrast, with varying through-the-thickness properties, this resulted in a wide range of cation distribution within the specimen, ranging from a maximum value of 410 mol/m3 in the upper surface to a minimum value of 245 mol/m3 at z = 0 mm. It is of importance to note that this alteration in the distribution of cations is due to the different resistance that cartilage offers to flow passage, which depends on the Es, GAG distribution and k values. This is an essential factor to consider when designing biomaterials for tissue repair, since it directly influences the distribution of nutrients and oxygen in the articular cartilage.

thumbnail
Fig 8. Range of cation concentration (c+) distribution in sample depth for the free swelling test considering cartilage through-the-thickness heterogeneity.

https://doi.org/10.1371/journal.pone.0157967.g008

Fluxes and swelling.

The influence of anisotropy on cation fluxes and their distribution within the tissue was also considered. The simulation results are again shown again at 1650 seconds of the swelling for both numerical situations (through-the-thickness varying properties and homogeneous), corresponding to the main entrance of water (Fig 7). These results demonstrate that, similar to swelling, important variations in water and cation fluxes appear when considering anisotropy. In Fig 9, the heterogeneous results show how areas with lower Es exhibit the highest entrance of water into the sample, 2.5·10−9 m3/s and subsequent higher z-displacements, whereas stiffer areas offer higher resistance to water entrance. Within the sample, the flux was reduced, the minimum being at the bottom of the sample (z = 0 mm). Biologically, this Es variation within the cartilage is correlated with the collagen fibre distribution which enables resistance to higher shear stress. The GAG content influences the compression capacity of articular cartilage in the joint regions where direct cartilage-cartilage contact is found. Regarding the cation flux, the simulation shows how the cations flow out of the sample. Note that negative values indicate outflow. In the cartilage surface samples, the maximum and minimum outflows of cations were 9.57·10−4 mol/s and 2.37·10−4 mol/s, respectively. In Fig 9, the simulation results corresponding to the homogeneous case show homogeneous layers of water and cation fluxes. As expected, a regular water flux occurred in the cartilage, reaching its maximum value of 2.5·10−9 m3/s at the upper surface and decreasing towards the bottom of the sample. A maximum outgoing flux of 1.2·10−3 mol/s was observed in the upper surface. In both cases, anisotropy and isotropy water and cation fluxes were consistent with the external bath variation (from 0.15 m to 0.125 m). To equilibrate the imbalance between the internal sample medium and the external bath solution, an incoming flux of water and an outflow of cation fluxes was generated. These results show significant differences in water and cation fluxes when modelling cartilage as a heterogeneous material or as a simplified homogeneous material.

thumbnail
Fig 9. z-displacement of the upper surface of the cartilage, and the flux of water and cations obtained with the model after 1650 seconds (coincident with the maximum swelling of the sample) of the free swelling simulation for through-the-thickness heterogeneous materials and homogeneous materials.

Note that the positive flux refers to the entrance of water or cations into the sample and the negative flux indicates the outflow of the substance.

https://doi.org/10.1371/journal.pone.0157967.g009

Conclusions

This work presents an extension of a previously developed 3D mechano-electrochemical model used to analyse and quantify the effects of variations in through-the-thickness as well as radial properties on articular cartilage behaviour. The model is capable of predicting swelling and flow of water and cations as well as their distributions within samples when external bath solution varied. Furthermore, it enables us to quantify the influence of variations in the mechano-electrochemical properties on articular cartilage behaviour. This, in combination with its capability of displaying the results in easily interpretable 3D images, makes this model an interesting novel tool for designing and assessing the efficacy of biomaterials for cartilage repair. Additionally, it allows the establishment of baseline patterns from which tissue degradation due to disease can be assessed. The main novelty of the model is the incorporation of experimentally measured heterogeneous properties (Es and GAGs distribution) which have been shown not only to vary in depth but also in the radial direction. In the absence of data in the literature regarding an accurate Es distribution along the surface, we performed a systematic indentation test with the WASIM, an instrument specially designed to perform normal indentation on intact articular surfaces. A qualitative GAG distribution was obtained from histological assays of porcine specimens. The results show how the Es as well as the GAG distribution increased in latero-posterior areas, coincident with zones that support higher compression and shear stress. These properties were introduced into the 3D computational model to determine their effects on the mechano-electrochemical events occurring in articular cartilage. The results obtained demonstrate that when considering cartilage as a heterogeneous material with variations in its through-the-thickness properties, the model predicts marked differences in swelling, fluxes and cation distribution compared to those obtained in homogeneous material simulations. There are significant reductions in the swelling and in the water and ion flux values compared to those previously obtained with samples having homogeneous properties. The computational results also showed that the upper surface displacement was twice as high in the isotropic case compared to the heterogeneous case where the stiffer areas offer higher resistance to material deformation. Similarly, the water and ion fluxes in the homogeneous case lead to higher flows that provoke the higher swelling. In the case of the cation distribution within the sample, the simulation evidenced how the heterogeneous case offers a wide range of cation distributions whereas in the homogeneous case, as expected, a homogeneous distribution within each layer was obtained. This indicates that fluxes are greatly affected by material heterogeneity and anisotropy providing areas with higher concentration of cations. In the light of the results obtained, mediated diffusive processes can be said to be highly influenced by the variation in properties generating zones with greater nutrient and oxygen up-take. It is therefore important to take into account through-the-thickness as well as radial variations in the material properties when manufacturing cartilage-mimicking biomaterials in order to achieve success in cartilage repair process.

Acknowledgments

The authors gratefully acknowledge the financial support from the Spanish Ministry of Economy and Competitiveness (MINECO MAT2013-46467-C4-3-R and FPU graduate research program AP2010/2557), the Government of Aragon (DGA) and the Biomedical Research Networking Center in Bioengineering, Biomaterials and Nanomedicine (CIBER-BBN). CIBER-BBN is financed by the Instituto de Salud Carlos III with the assistance from the European Regional Development Fund. The authors also acknowledge Nuffield Orthopaedic Centre NIHR Biomedical Research Unit funding and Orthopaedic Research of UK (ORUK).

Author Contributions

Conceived and designed the experiments: SM MA AJP PAH HSG MD MHD. Performed the experiments: SM MA AJP PAH HSG MD MHD. Analyzed the data: SM MA AJP PAH HSG MD MHD. Contributed reagents/materials/analysis tools: SM MA AJP PAH HSG MD MHD. Wrote the paper: SM MA AJP PAH HSG MD MHD.

References

  1. 1. Tomkoria S, Patel RV, Mao JJ. Heterogeneous nanomechanical properties of superficial and zonal regions of articular cartilage of the rabbit proximal radius condyle by atomic force microscopy. Med Eng Phys, 2004. 26(10): 815–822. pmid:15567698
  2. 2. Jurvelin JS, Buschmann MD, Hunziker EB. Mechanical anisotropy of the human knee articular cartilage in compression. Proc Inst Mech Eng H, 2003. 217(3): 215–219. pmid:12807162
  3. 3. Zheng SK, Xia Y, Badar F. Further Studies on the Anisotropic Distribution of Collagen in Articular Cartilage by mu MRI. Magnetic Resonance in Medicine, 2011. 65(3): 656–663. pmid:20939069
  4. 4. Mow V, Guo XE. Mechano-electrochemical properties of articular cartilage: Their inhomogeneities and anisotropies. Annual Review of Biomedical Engineering, 2002. 4: 175–209. pmid:12117756
  5. 5. Wang CC, Hung CT, Mow VC. An analysis of the effects of depth-dependent aggregate modulus on articular cartilage stress-relaxation behavior in compression. J Biomech, 2001. 34(1): 75–84. pmid:11425083
  6. 6. Federico S, Grillo A, La Rosa G, Giaquinta G, Herzog W. A transversely isotropic, transversely homogeneous microstructural-statistical model of articular cartilage. J Biomech, 2005. 38(10): 2008–2018. pmid:16084201
  7. 7. Garcia JJ, Cortes DH. A nonlinear biphasic viscohyperelastic model for articular cartilage. J Biomech, 2006. 39(16): 2991–2998. pmid:16316659
  8. 8. Wilson W, Van Donkelaar CC, Van Rietbergen R, Huiskes R. The role of computational models in the search for the mechanical behavior and damage mechanisms of articular cartilage. Med Eng Phys, 2005. 27(10): 810–826. pmid:16287601
  9. 9. Broom ND, Oloyede A. The importance of physicochemical swelling in cartilage illustrated with a model hydrogel system. Biomat, 1998. 19(13): 1179–1188.
  10. 10. Martin JA, Buckwalter JA. Aging, articular cartilage chondrocyte senescence and osteoarthritis. Bioger, 2002. 3(5): 257–264.
  11. 11. Watson PJ, Carpenter TA, Hall LD, Tyler JA. Cartilage swelling and loss in a spontaneous model of osteoarthritis visualized by magnetic resonance imaging. Osteoarthritis and Cartilage, 1996. 4(3): 197–207. pmid:8895221
  12. 12. Manzano S, Gaffney EA, Doblare M, Hamdy Doweidar M. Cartilage Dysfunction in ALS Patients as Side Effect of Motion Loss: 3D Mechano-Electrochemical Computational Model. BioMed Research International, 2014.
  13. 13. Manzano S, Manzano R, Doblaré M, Doweidar MH. Altered swelling and ion fluxes in articular cartilage as a biomarker in osteoarthritis and joint immobilization: a computational analysis. Journal of The Royal Society Interface, 2015. 12(102): 20141090.
  14. 14. Manzano S, Doblaré M, Doweidar MH. Parameter-dependent behavior of articular cartilage: 3D mechano-electrochemical computational model. Computer Methods and Programs in Biomedicine, 2015. 112(3): 491–502.
  15. 15. Maroudas A, Venn M. (1977). Chemical composition and swelling of normal and osteoarthrotic femoral head cartilage. II. Swelling. Annals of the rheumatic diseases, 36(5), 399–406. pmid:200188
  16. 16. Summers GC, Merrill A, Sharif M, Adams MA. (2008). Swelling of articular cartilage depends on the integrity of adjacent cartilage and bone.Biorheology, 45(3), 365.
  17. 17. Wang Q, Yang YY, Niu HJ, Zhang WJ, Feng QJ, Chen WF. (2013). An ultrasound study of altered hydration behaviour of proteoglycan-degraded articular cartilage. BMC musculoskeletal disorders, 14(1), 289.
  18. 18. Richard F, Villars M, Thibaud S. Viscoelastic modeling and quantitative experimental characterization of normal and osteoarthritic human articular cartilage using indentation. J Mech Behav Biomed Mater, 2013. 24: 41–52. pmid:23684353
  19. 19. Calvo E, Palacios I, Delgado E, Sánchez-Pernaute O, Largo R, Egido J, et al. Histopathological correlation of cartilage swelling detected by magnetic resonance imaging in early experimental osteoarthritis. Osteoarthritis Cartilage, 2004. 12(11): 878–886. pmid:15501403
  20. 20. Huey DJ, Hu JC, Athanasiou KA. Unlike bone, cartilage regeneration remains elusive. Science, 2012. 338(6109): 917–921. pmid:23161992
  21. 21. Guo H, Maher SA, Torzilli PA, A biphasic multiscale study of the mechanical microenvironment of chondrocytes within articular cartilage under unconfined compression. J Biomech, 2014.
  22. 22. Zhang J, Wang X. The three-dimensional simulation of arytenoid cartilage movement. Lin Chung Er Bi Yan Hou Tou Jing Wai Ke Za Zhi, 2011. 25(15): 687–689. pmid:22010338
  23. 23. Babalola OM, Bonassar LJ. Parametric finite element analysis of physical stimuli resulting from mechanical stimulation of tissue engineered cartilage. J Biomech Eng, 2009. 131(6): 061014. pmid:19449968
  24. 24. Rasanen LP, Mononen ME, Nieminen MT, Lammentausta E, Jurvelin JS, Korhonen RK, et al. Implementation of subject-specific collagen architecture of cartilage into a 2D computational model of a knee joint—data from the Osteoarthritis Initiative (OAI). J Orthop Res, 2013. 31(1): 10–22. pmid:22767415
  25. 25. Lai WM, Hou JS, Mow VC. A triphasic theory for the swelling and deformational behaviors of articular cartilage. J Biomech Eng, 1991. 113: 245–258. pmid:1921350
  26. 26. Pharr GM, Oliver WC, Brotzen FR. On the generatlity of the relationship among contact stiffness, contact area and elastic modulus during indentation. Journal of Material Research. 2011, 7(1): 613–617.
  27. 27. Field JS, Swain MV. A simple predictive model for spherical indentation. Journal of Materials Research, 1993, 8(2): 297–306.
  28. 28. Wang T, Wen CY, Yan CH, Lu WW, Chiu KY. Spatial and temporal changes of subchondral bone proceed to microscopic articular cartilage degeneration in guinea pigs with spontaneous osteoarthritis. Osteoarthritis Cartilage, 2013. 21(4): 574–81. pmid:23313833
  29. 29. Li LP, Buschmann MD, Shirazi-Adl A. A fibril reinforced nonhomogeneous poroelastic model for articular cartilage: inhomogeneous response in unconfined compression. J Biomech, 2000. 33(12): 1533–1541. pmid:11006376
  30. 30. Albro MB, Nims RJ, Durney KM, Cigan AD, Shim JJ, Vunjak-Novakovic G, et al. Heterogeneous engineered cartilage growth results from gradients of media-supplemented active TGF-β and is ameliorated by the alternative supplementation of latent TGF-β. Biomaterials, 2016. 77: 173–185. pmid:26599624
  31. 31. Duarte Campos DF, Drescher W, Rath B, Tingart M, Fischer H. Supporting Biomaterials for Articular Cartilage Repair. Cartilage, 2012: 205–221. pmid:26069634
  32. 32. Sun DN, Gu WY, Guo XE, Lai WM, Mow VC. A mixed finite element formulation of triphasic mechano-electrochemical theory for charged, hydrated biological soft tissues. International Journal for Numerical Methods in Engineering, 1999. 45(10): 1375–1402.
  33. 33. Zadpoor A. Nanomechanical characterization of heterogeneous and hierarchical biomaterials and tissues using nanoindentation: the role of finite mixture model. Mater Sci Eng C Mater Biol Appl., 2015 48:150–157. pmid:25579908
  34. 34. Armstrong CG, Lai WM, Mow VC. An analysis of the unconfined compression of articular cartilage. Journal of biomechanical engineering, 1984, 106(2), 165–173. pmid:6738022
  35. 35. Korhonen RK, Wong M, Arokoski J, Lindgren R, Helminen HJ, Hunziker EB, et al. Importance of the superficial tissue layer for the indentation stiffness of articular cartilage. Med Eng Phys, 2002. 24(2): 99–108. pmid:11886828
  36. 36. Schinagl RM, Gurskis D, Chen AC, Sah RL. Depth-dependent confined compression modulus of full-thickness bovine articular cartilage. J Orthop Res, 1997. 15(4): 499–506. pmid:9379258