A validated model of the pro- and anti-inflammatory cytokine balancing act in articular cartilage lesion formation

1. Introduction

The degenerative joint disease known as osteoarthritis (OA) is among the most common causes of disability worldwide. While OA involves multiple joint tissues including bone, tendons, ligaments and synovium, articular cartilage degeneration, and erosion is the proximal cause of loss of joint function. Articular cartilage is a thin layer of connective tissue that covers the ends of long bones in synovial joints such as the shoulder, hip, knee, and ankle, where it distribute mechanical loads and allows for smooth joint motion. These functions are attributable to the unique composition and structure of cartilage extracellular matrix (ECM), which consists of water (>70%), proteoglycan (15%), and collagen (15%) (Martin et al., 2011). Aggrecan, the major cartilage proteoglycan, is heavily decorated with negatively charged sulfate groups, which retain water in the tissue. Large complexes of aggrecan are trapped within the matrix by a collagen fibril network. The resistance of cartilage to compression and its ability to distribute loads is largely due to this macromolecular arrangement, and aggrecan depletion or collagen degradation radically reduces mechanical strength (Farndale et al., 1982; Lu et al., 2011).

Even though articular cartilage is only 1–3 mm thick, it has four distinct zones (superficial, transitional, radial, calcified). Different zones have different cell morphologies, matrix composition, and collagen fibril properties. Compared to other zones, the superficial zone is specialized to resist tensile stresses and minimize surface friction. In this paper, we focus on the properties of the superficial zone. We assume the solid cartilage matrix to be a homogeneous system, so that the chemical and cell properties remain the same inside the space. Cartilage cells, known as chondrocytes, are distributed sparsely within the tissue (104 cells/mm³). They are largely trapped inside the ECM, so there is no appreciable cell motility. Chondrocytes are solely responsible for the maintenance of the cartilage matrix, and engage in complex biochemical signaling/regulation via the synthesis and release or recognition of signaling molecules such as cytokines and oxidants, among others. We classify chondrocytes to be in different states in our model with respect to the chemical signaling processes. In addition to synthesizing ECM components, chondrocytes also release matrix proteases that cause matrix degradation.

There are various ways to cause traumatic articular cartilage injury, but they all share high loading rates and a high peak stress amplitude, which initiates the damage (lesion). The damage done by injuries seldom heals spontaneously and often leads to the progressive cartilage degeneration characteristic of post-traumatic osteoarthritis (PTOA). The strain in the superficial zone under a single blunt impact can easily exceed 40%, and when combined with the high loading rate and excessive stress, is lethal to chondrocytes and detrimental to the ECM. Chondrocyte death inside the superficial zone can be assumed to be directly caused by this impact. The main focus of most therapies for PTOA has been to prevent chondrocyte death and dysregulation. The chondrocyte depletion and ECM degradation process is illustrated in Figure 1. Cartilage damage initiates the production of alarmins, such as damage associated molecular patterns (DAMPs) (Bianchi, 2007), which can induce the release of pro-inflammatory cytokines such as interleukin-6 (IL-6), TNF-α, IL-1 α, and IL-1 β. A previous microarray study suggests that TNF-α might be the early mediator; however, IL-1 β is the cytokine that sustains the degradation (Martin et al., 2011). In this model, we use IL-6 to represent the entire cytokine family for model simplicity. Other cytokines such as TNF-α or IL-1 β can be easily added to the equations, but doing so would provide limited benefit with the need for additional parameter estimation and model complexity.

FIGURE 1

Figure 1. Schematic of the articular cartilage lesion formation process due to a single blunt impact.

Pro-inflammatory cytokines are the main reason for cell apoptosis. Moreover, the pro-inflammatory cytokines can cause severe aggrecan depletion, which leads to loss of strength and elevated strain in affected cartilage. Even though it is quite limited, chondrocytes still have some self-repair ability. Anti-inflammatory cytokines such as erythropoietin (EPO) can antagonize the effect of pro-inflammatory cytokines and result in reduced cell apoptosis and ECM degradation. The “balancing act” (Graham et al., 2012) between pro-inflammatory and anti-inflammatory cytokines is the essential mechanism determining whether a cartilage lesion will heal or expand in size.

This paper is organized as follows. We describe the mathematical models and numerical methods used to solve the model equations. We then describe the materials and methods used for the experimental validation of our computational results. We then present the computational results and the experimental validation. We finish with a discussion of these results.

2. Materials and Methods

2.1. Mathematical Model and Numerical Methods

In this section, we describe the age- and space-structured model developed for the inflammatory response after a single blunt impact injury. The cartilage lesion caused by a single severe traumatic event was described in a reaction-diffusion-delay model by Graham et al. (2012) and Graham (2013). In our model, we use age structure instead of delay terms to model the delayed transitions between cell states. This change in modeling approach was necessitated by the increased complexity of the interactions we are representing mathematically. Age structure is more efficient computationally than the use of delay terms, and is more flexible in representing transitions between cell states.

2.1.1. Components of the system

We assume circular symmetry so that the system can be reduced to a one-dimensional model with respect to space. The components of the system depend on radius (r), age (a), and time (t). The time scales for our cell transitions are denoted by τ₁ and τ₂. We also assume that the initial blunt impact occurred on a small region near the origin with radius smaller than 0.25 cm, in a spatial domain with a total radius of 2.5 cm.

There are two main categories of components in our mathematical model, cells and chemicals. A schematic of the system is presented in Figure 1.

Let r denote radius, a time since signaled (i.e., age in current state), and t time (i.e., wall clock). The cell components are

• C_U (r, t) = population density (cells per unit area) of healthy cells not yet signaled by ROS.

• C_T (r, a, t) = population density of healthy cells signaled by ROS and in the process of becoming catabolic.

• C_E (r, t) = population density of healthy cells signaled by ROS and starting to produce EPO (roughly 20–24 h after being signaled).

• S_T (r, a, t) = population density of sick cells in the catabolic state. Healthy cells signaled by alarmins (DAMPs) and inflammatory cytokines (IL-6, TNF-α) enter into the catabolic state. Catabolic cells start to synthesize inflammatory cytokines and ROS.

• S_A (r, t) = population density of EPOR-active cells. Catabolic cells signaled by inflammatory cytokines enter the EPOR-active state. EPOR-active cells express a receptor (EPOR) for EPO; however, there is a 8–12 h time gap before a cell expresses the EPO receptor after it was signaled to become EPOR-active (Brines and Cerami, 2008). After EPOR-active cells express a receptor for EPO, they may switch back to the healthy state C_U when signaled by EPO.

• D_A (r, t) = population density of apoptotic cells. Catabolic cells signaled by alarmins and inflammatory cytokines enter the apoptotic state. EPOR-active cells will also become apoptotic after signaled by inflammatory cytokines. Apoptotic cells do not play an explicit role in the mathematical model.

• D_N (r, t) = population density of necrotic cells. In our model, necrotic cells are only created by the initial blunt impact and release alarmins (DAMPs) into the system. Before cells become necrotic, cells release a small amount of ROS. Fully necrotic cells cannot produce ROS and are eventually cleared from the system.

We assume negligible motility on the part of chondrocytes, so there are no diffusion terms for the cells equations. The different cell states correspond to the different chemical signals. The injury kills cells inside the impact area by rendering them necrotic (D_N). The other cells adjacent to the impact area transit from “healthy” (C_U, C_T, C_E) to “sick” (S_T, S_A) and then “dead” (apoptotic D_A) states. Since apoptotic cells (D_A) do not play an explicit role in the system, they are not expressed explicitly in our mathematical model, but are instead represented by the sink terms inside the S_T and S_A equations.

The chemical components are

• R (r, t) = concentration of reactive oxygen species (ROS). ROS triggers healthy cells to change states.

• M (r, t) = concentration of alarmins (DAMPs) released by necrotic cells and ECM degradation. DAMPs signal healthy cells C_T to enter the catabolic state S_T, and together with inflammatory cytokines such as IL-6, cause catabolic cells S_T to become apoptotic.

• F (r, t) = concentration of inflammatory cytokines (IL-6) produced by catabolic cells (S_T). These inflammatory cytokines

– signals healthy cells (C_T) to enter the catabolic state (S_T),

– signal catabolic cells (S_T) to enter the EPOR-active state (Brines and Cerami, 2008),

– cause both catabolic and EPOR-active cells to become apoptotic,

– degrade the extracellular matrix, which in turn increases the level of DAMPs, resulting in further damage of the system. The degradation of ECM is a slow and complex process. However, we assume for mathematical convenience that inflammatory cytokines directly damage ECM,

– limit the production of EPO.

• P (r, t) = concentration of erythropoietin (EPO). EPO is produced exclusively by C_E in our model. Inflammatory cytokines suppress this process. EPO helps EPOR-active cells (S_A) switch back to the healthy state C_U. The effects of EPO depend on its concentration. When the concentration of EPO passes the threshold P_c (Brines and Cerami, 2008), the spread of inflammation can be slowed by terminating the effect of inflammatory cytokines and alarmins on the system. We also assume that C_E cells revert to C_U when the EPO level exceeds the P_c threshold.

We assume that the chemicals diffuse throughout the whole region. The diffusion coefficients were estimated by Graham et al. (2012) and Leddy and Guilak (2003). Chemicals decay after a certain time, and the decay rate can be approximated by their half lives (Eckardt et al., 1989; Wedlock et al., 1996; Ito et al., 2008). However, the decay of ROS is almost instantaneous under the superoxide dismute SOD, so the decay rate needs to be adjusted to fit the mathematical model. The inflammatory cytokines such as IL-6 and TNF-α are the main cause of cartilage lesion formation. EPO plays an opposing role by helping cell recovery and limiting the inflammation (Brines and Cerami, 2008). Our model captures the balance between these pro-inflammatory and anti-inflammatory cytokines.

In addition to the chemicals, we track the extracellular matrix density:

• U(r, t) = density of extracellular matrix (ECM). ECM is degraded by inflammatory cytokines and releases alarmins. The degradation of ECM is measured by the decrease in the concentration of SO₄ (Farndale et al., 1982).

2.1.2. Model equations

The equations for the chemical components of our system are

$\begin{array}{ll}\hfill {\partial }_t\underset{\mathrm{ROS}}{\underbrace{R(r,t)}}& =\underset{\mathrm{diffusion}}{\underbrace{\frac{1}{r}{\partial }_r\left(r{D}_R{R}_r\right)}}-\underset{\mathrm{natural}\mathrm{decay}}{\underbrace{{\mathrm{\delta }}_{R}R}}+\underset{\mathrm{production}\mathrm{by}{S}_T}{\underbrace{{\mathrm{\sigma }}_R{S}_T}},\end{array}$ (1a)

$\begin{array}{ll}\hfill {\partial }_t\underset{\mathrm{DAMPs}}{\underbrace{M\left(r,t\right)}}& =\underset{\mathrm{diffusion}}{\underbrace{\frac{1}{r}{\partial }_r\left(r{D}_M{M}_r\right)}}-\underset{\mathrm{natural}\mathrm{decay}}{\underbrace{{\mathrm{\delta }}_{M}M}}+\underset{\mathrm{production}\mathrm{by}{D}_N}{\underbrace{{\mathrm{\sigma }}_M{D}_N}}+\underset{\mathrm{production}\mathrm{by}\mathrm{ECM}}{\underbrace{{\mathrm{\delta }}_{U}U\frac{F}{{\mathrm{\lambda }}_F+F}}},\end{array}$ (1b)

$\begin{array}{ll}\hfill {\partial }_t\underset{\mathrm{IL}-6}{\underbrace{F\left(r,t\right)}}& =\underset{\mathrm{diffusion}}{\underbrace{\frac{1}{r}{\partial }_r\left(r{D}_F{F}_r\right)}}-\underset{\mathrm{natural}\mathrm{decay}}{\underbrace{{\mathrm{\delta }}_{F}F}}+\underset{\mathrm{production}\mathrm{by}{S}_T}{\underbrace{{\mathrm{\sigma }}_F{S}_T}},\end{array}$ (1c)

$\begin{array}{l}\hfill \hfill {\partial }_t\underset{\mathrm{EPO}}{\underbrace{P\left(r,t\right)}}& =\underset{\mathrm{diffusion}}{\underbrace{\frac{1}{r}{\partial }_r\left(r{D}_P{P}_r\right)}}-\underset{\mathrm{natural}\mathrm{decay}}{\underbrace{{\mathrm{\delta }}_{P}P}}+\underset{\mathrm{production}\mathrm{by}{C}_E,\mathrm{controlled}\mathrm{by}F}{\underbrace{{\mathrm{\sigma }}_P{C}_E\frac{R}{{\mathrm{\lambda }}_R+R}\frac{\Lambda }{\Lambda +F}}}.\end{array}$ (1d)

The initial and boundary conditions are

${\partial }_{r}R\left(0,t\right)={\partial }_{r}M\left(0,t\right)={\partial }_{r}F\left(0,t\right)={\partial }_{r}P\left(0,t\right)=0,$

R (r, 0) = σ_R · D_N (r, 0), and M (r, 0) = F (r, 0) = P (r, 0) = 0.

Our chemical equations are similar to those in Graham et al. (2012). A very significant change is that the system is triggered by ROS released by cells as they become necrotic. The initial condition of ROS is not zero. We use the diffusion and natural decay rates of chemicals as estimated in Wang et al. (2014).

The ECM is assumed to be degraded by inflammatory cytokines such as IL-6, measured in terms of decreased proteoglycan concentration in the matrix. When ECM is intact, the sulfate groups are kept inside the ECM. The release of sulfate groups is an indication of ECM degradation, which can be estimated by the decrease in concentration of SO₄. The average concentration of SO₄ in normal undamaged cartilage is 30 g/L (Farndale et al., 1982), which is the initial density of ECM in this model. EPO concentration also affects ECM degradation.

We define the Heaviside function,

$H\left(\mathrm{\theta }\right)=\left\{\begin{array}{cc}1,\hfill & \mathrm{\theta }>0,\hfill \\ 0,\hfill & \mathrm{\theta }<1.\hfill \end{array}\right.$ (2)

The equation for the ECM dynamics is

${\partial }_t\underset{\mathit{ECM}}{\underbrace{U\left(r,t\right)}}=\underset{\mathrm{degradated}\mathrm{by}\mathrm{alarmins},\mathrm{checked}\mathrm{by}\mathrm{EPO}}{\underbrace{-{\mathrm{\delta }}_{U}U\frac{F}{{\mathrm{\lambda }}_F+F}H\left(P_c-P\right)}},$ (3)

with initial condition U(r, 0) = 30 mg. The Heaviside function H(P_c − P) represents the property that ECM degradation can be terminated when the level of P exceed P_c.

A substantial difference between our model and Graham et al. (2012) is that we use age structure to represent the time delays for cells to become EPO producing and EPOR-active, instead of delay terms. This gives us a system that is more tractable computationally and more flexible in representing the transitions. Another major difference is that we capture more fully the dynamics among healthy cells. We separate healthy cells into three states, C_U, C_T, and C_E. Since there are only two transitions, C_T → C_E and S_T → S_A, that involve a time delay, we include age structure a in only the equations for C_T and S_T. Following Ayati (2007a), we capture the sharp age of transition using the function

$\mathrm{\gamma }\left(a-a_{\max }\right)=\frac{{\mathrm{\gamma }}_0}{\mathrm{\sigma }}\left(\tanh \left(\frac{a-a_{\max }}{\mathrm{\sigma }}\right)+1\right),$ (4)

where σ is the spread parameter and γ₀ is the height parameter. For a fixed γ₀ and σ → 0, cells C_T and S_T switch their states to C_E and S_A all at the same age a_max. When σ is small but not zero, this behavior is smoothed out, better representing the underlying biology. For definiteness, we choose a_max = 1 for C_T → C_E and $a_{\max }=\frac{1}{2}$ for S_T → S_A.

The equation for the population density of healthy cells not yet signaled by ROS is

${\partial }_t{C}_U\left(r,t\right)=\underset{S_A\underset{}{\overset{\mathrm{EPO}}{\to }}{C}_U}{\underbrace{{\mathrm{\alpha }}_1{S}_A\frac{P}{{\mathrm{\lambda }}_P+P}}}+\underset{C_E\underset{}{\overset{\mathrm{EPO}}{\to }}{C}_U}{\underbrace{{\mathrm{\alpha }}_{2}H\left(P-P_c\right)C_E}}-\underset{C_U\underset{}{\overset{\mathrm{ROS}}{\to }}{C}_T}{\underbrace{{\mathrm{\beta }}_{13}{C}_U\frac{R}{{\mathrm{\lambda }}_R+R}}},$ (5a)

with initial condition

$C_U\left(r,0\right)=\left\{\begin{array}{cc}0,\hfill & 0\le r\le 0.25\mathrm{cm}\hfill \\ 100000\mathrm{cells}∕\mathrm{c}{\mathrm{m}}^2,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$ (5b)

The equation for the population density of healthy cells signaled by ROS and in the process of becoming catabolic is

$\begin{array}{llll}\hfill & {\partial }_t{C}_T\left(r,a,t\right)+{\partial }_a{C}_T\left(r,a,t\right)=& \hfill & \\ \hfill & -\underset{C_T\underset{}{\overset{\mathrm{DAMPs}}{\to }}{S}_T}{\underbrace{{\mathrm{\beta }}_{11}\frac{M}{{\mathrm{\lambda }}_M+M}H\left(P_c-P\right)C_T\left(r,a,t\right)}}-\underset{C_T\underset{}{\overset{\mathrm{TNF}-\mathrm{\alpha }}{\to }}{S}_T}{\underbrace{{\mathrm{\beta }}_{12}\frac{F}{{\mathrm{\lambda }}_F+F}H\left(P_c-P\right)C_T\left(r,a,t\right)}}& \hfill & \\ \hfill & -\underset{C_T\underset{}{\overset{\mathrm{ROS}}{\to }}{C}_E}{\underbrace{{\mathrm{\kappa }}_2\mathrm{\gamma }\left(a-{\mathrm{\tau }}_2\right)\frac{R}{{\mathrm{\lambda }}_R+R}{C}_T\left(r,a,t\right)}},\end{array}$ (6a)

with age boundary condition

$C_T\left(r,0,t\right)=\underset{C_U\underset{}{\overset{\mathrm{ROS}}{\to }}{C}_T}{\underbrace{{\mathrm{\beta }}_{13}{C}_U\frac{R}{{\mathrm{\lambda }}_R+R}}},$ (6b)

and initial condition C_T (r, a, 0) = 0.

The equation for the population density of healthy cells signaled by ROS and producing EPO is

$\begin{array}{ll}\hfill {\partial }_t{C}_E\left(r,t\right)& =\underset{C_T\underset{}{\overset{{\mathrm{\tau }}_2\mathrm{delay}}{\to }}{C}_E}{\underbrace{{\int }_0^{\infty }{\mathrm{\kappa }}_2\mathrm{\gamma }\left(a-{\mathrm{\tau }}_2\right)\frac{R\left(r,t\right)}{{\mathrm{\lambda }}_R+R\left(r,t\right)}{C}_T\left(r,a,t\right)\mathit{da}}}-\underset{C_E\underset{}{\overset{\mathrm{EPO}}{\to }}{C}_U}{\underbrace{{\mathrm{\alpha }}_{2}H\left(P-P_c\right)C_E}},\end{array}$ (7)

with initial condition C_E (r, 0) = 0.

The equations for the population density of sick cells in the catabolic state is

$\begin{array}{ll}\hfill {\partial }_t{S}_T\left(r,a,t\right)+{\partial }_a{S}_T\left(r,a,t\right)& =-\underset{S_T\underset{}{\overset{\mathrm{TNF}-\mathrm{\alpha },\mathrm{DAMPs}}{\to }}{D}_A}{\underbrace{{\mathrm{\mu }}_{\mathit{ST}}\frac{F}{{\mathrm{\lambda }}_F+F}\frac{M}{{\mathrm{\lambda }}_M+M}{S}_T\left(r,a,t\right)}}-\underset{S_T\underset{}{\overset{\mathrm{TNF}-\mathrm{\alpha }}{\to }}{S}_A}{\underbrace{\mathrm{\eta }\cdot \mathrm{\gamma }\left(a-{\mathrm{\tau }}_1\right)\frac{F}{{\mathrm{\lambda }}_F+F}{S}_T\left(r,a,t\right)}},\end{array}$ (8a)

with age boundary condition

$S_T\left(r,0,t\right)={\displaystyle {\int }_0^{\infty }\underset{C_T\stackrel{\text{DAMPs}}{\to }{S}_T}{\underbrace{({\beta }_{11}{C}_U\frac{M}{{\lambda }_M+M}H(P_c-P)}}+\underset{C_T\stackrel{\text{TNF}-\alpha }{\to }{S}_T}{\underbrace{{\beta }_{12}{C}_U\frac{F}{{\lambda }_F+F}H(P_c-P)}})C_T(r,a,t)}da,$ (8b)

and initial condition S_T (r, a, 0) = 0.

The equation for the population density of EPOR-active sick cells is

$\begin{array}{l}\hfill {\partial }_t{S}_A\left(r,t\right)=\underset{S_T\underset{}{\overset{\mathrm{TNF}-\mathrm{\alpha }}{\to }}{S}_A}{\underbrace{{\int }_0^{\infty }\mathrm{\eta }\cdot \mathrm{\gamma }\left(a-{\mathrm{\tau }}_1\right)\frac{F\left(r,t\right)}{{\mathrm{\lambda }}_F+F\left(r,t\right)}{S}_T\left(r,a,t\right)\mathit{da}}}-\underset{S_A\underset{}{\overset{\mathrm{EPO}}{\to }}{C}_U}{\underbrace{{\mathrm{\alpha }}_1{S}_A\frac{P}{{\mathrm{\lambda }}_P+P}}}-\underset{S_A\underset{}{\overset{\mathrm{TNF}-\mathrm{\alpha }}{\to }}{D}_A}{\underbrace{{\mathrm{\mu }}_{S_A}\frac{F}{{\mathrm{\lambda }}_F+F}{S}_A}},\end{array}$ (9)

with initial condition S_A (r, 0) = 0.

The equation for the necrotic cell population is

${\partial }_t{D}_N\left(r,t\right)=-\underset{\mathrm{natural}\mathrm{decay}}{\underbrace{{\mathrm{\mu }}_{D_N}{D}_N}},$ (10a)

with initial condition

$D_N\left(r,0\right)=\left\{\begin{array}{cc}100000\mathrm{cells}∕\mathrm{c}{\mathrm{m}}^2,\hfill & 0\le r\le 0.25\mathrm{cm},\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$ (10b)

2.1.3. Numerical methods

The computational methods used to solve the age- and space-structured differential equations were developed and analyzed in Ayati (2000, 2007a) and Ayati and Dupont (2002, 2009) and have been used effectively and efficiently in the modeling and simulation of biofilms (Ayati and Klapper, 2007, 2012; Klapper et al., 2007; Ayati, 2011), avascular tumor invasion (Ayati et al., 2006), and Proteus mirabilis swarm colony development (Ayati, 2006, 2007b, 2009). Numerical efficiency is the main reason that we use age structure rather than delay terms in the model in this paper. To solve the model equations in Graham et al. (2012) and Wang et al. (2014), the authors first did a semi-discretization in space and then solved the resulting system of ordinary delay differential equations using the Matlab function dde23 (Shampine and Thompson, 2001). The function dde23 uses an explicit time integration method. As a result, the stability constraints of the semi-discretized system result in very small time steps. Although another numerical method for delay differential equations could be constructed without such constraints, the need to store the past history of a large spatial system renders the delay approach much less efficient than using age structure; the extra dimension of age is accurately discretized with far less data storage and computation (we use 129 age intervals in the simulations in this paper). The time discretization is adaptive (Ayati and Dupont, 2005) and controlled by a tolerance parameters. We use center finite differences to discretize the diffusion terms.

The relative errors obtained in our computational convergence studies of the numerical results are given in Table 1. For simplicity, we vary the number of spatial intervals. We see that even for a fairly coarse discretization using 100 intervals, we have relative errors of <0.3%.

TABLE 1

Table 1. Table of relative errors.

2.2. Experimental Material and Methods

2.2.1. Generation of cartilage specimens

Osteochondral explants were surgically excised from bovine lateral tibial plateaus (25 mm × 25 mm × 10 mm). All specimens were allowed to equilibrate for two days in culture media in a low oxygen environment (5% O₂, 5% CO₂). After equilibration (day 0), explants either underwent a sham impact procedure or were impacted with a 5.5 mm rounded edge indenter from a drop tower, imparting an energy of 2.18 J/cm². Unimpacted explants were removed on day 0, preserved and embedded in paraffin. Impacted explants were placed back in fresh culture media, and removed at day 1, 7, or 14 for preservation and paraffin embedding. Culture media was changed every two days for explants cultured out to day 7 or 14.

2.2.2. Immunohistochemistry

Paraffin embedded cartilage explants were cut into 5 μm thick sections on glass slides to be processed for immunohistochemical detection of interleukin-6 (IL-6), erythropoietin (EPO), or the erythropoietin receptor (EPOR). All cartilage sections for each stain were processed in batch at the same time under identical conditions. The sections were deparaffinized, quenched of endogenous peroxidase activity, and then underwent antigen retrieval using a 0.01 M citric acid solution for 4 h at 65°C. Non-specific antigen binding was then blocked for 1 h with PBS containing 10% goat serum, 1% bovine serum albumin (w/v), and 0.1% tween 20. This blocking solution was then used to dilute the primary antibodies against IL-6 (Abcam ab193799), EPO (Abcam ab65394), and the EPOR (Abcam ab83696) at a ratio of 1:100 from their starting concentration. The solution containing the primary antibody was incubated on the slides overnight at 4°C. The slides were then rinsed and blocking solution applied for 30 min at room temperature. A biotinylated secondary antibody against the primary antibody was then incubated on the slides at a 1:250 dilution (Vector Labs). The slides were then rinsed and incubated with the VECTASTAIN^® ABC Reagent (Vector Labs) for 30 min at room temperature. After another rinsing, DAB (3, 3-diaminobenzidine) HRP substrate solution (Vector Labs DAB substrate kit) was incubated on the sections for 2–10 min (depending on primary antibody). The sections were then rinsed for a final time, counter stained with eosin, dehydrated, and a glass coverslip was attached and sealed with permount (Fisher).

Slides were then digitally scanned on an Olympus VS110. The Olympus software controlled the slide stage in all axes and automatically focused/captured high resolution images (322 nm/pixel) which were then tiled together, resulting in full explant immunohistochemically labeled images. The resolution of the resulting images was then reduced to 20% of their original size and exported as tiff images for analysis.

To obtain the validation data, we analyzed 12 cartilage slices. The slices were separated into six categories with two slices per category: EPO at days 1, 7, 14; and IL-6 at days 1, 7, 14. The length of each slice is around 1″ (2.5 cm approximately). The sham impact was placed in the middle of each cartilage slice and formed a dent (see Figure 2). Cells inside the cartilage slices were stained by the antibodies, and the number of cells per area corresponds to the relative concentrations of EPO and IL-6.

FIGURE 2

Figure 2. Slices of articular cartilage created at 1, 7, and 14 days after a single blunt impact. The impact site appears as a dent at the top of each image. The squares in each region are the zoomed areas shown in Figure 3. These slices have been stained to measure EPO concentrations, which are shown to increase as time progresses.

For image processing, we split each slice into two half-slices, divided at the point of impact. From these we could obtain four measurements of each of EPO and IL-6 at each radial value. We further subdivided each half-slice into pieces approximately 0.15 cm long; the pieces correspond roughly to the radius intervals [0, 0.15] cm, [0.15, 0.3] cm, [0.3, 0.45] cm, [0.45, 0.6] cm, [0.6, 0.75] cm, and [0.75, 0.9] cm.

We cropped an 800 × 800 pixel sample image from each piece and estimated the cell numbers in each sample. Figure 3 shows example samples from slices shown in Figure 2. The 800 × 800 pixel sample image has an area of 0.0166 cm². We have four such sample images per radius interval. Average cell count is used as the estimator of cell numbers in each radius interval. The EM algorithm and K-mean methods were applied for the cell segmentation. The upper threshold of the positive cell size is 50–150 μm² (approximately 20–60 pixels), while the lower threshold of the positive cell size is 25–30 μm² (approximately 10–15 pixels).

FIGURE 3

Figure 3. Zoomed areas from Figure 2. These samples have been stained to measure EPO concentrations, which are shown to increase as time progresses.

3. Results

3.1. Computational Results

We simulate the change of the chondrocyte population densities over a 14-day period after an initial cartilage injury in the center of a disk with a radius of 2.5 cm. The radius of the area of impact is 0.25 cm. We assume the initial cells density is 100,000 cells/cm².

The parameters used for the simulations are shown in Table 2. Most of the parameters are the same as in Wang et al. (2014), which contains detailed descriptions of how many of the parameters were determined from the literature. The parameters α₁, β₁₁, β₁₂, β₁₃, κ₁, and κ₂, which determine the cells transition rates, have been changed or added to account for the transition from delay terms to age structure in our model. Because our trigger P > P_c is never reached, the value of α₂ is irrelevant. We include it for completeness of the general model. Similarly, Λ is set sufficiently high that the effect of inflammatory cytokines on the transition to EPO-producing healthy cells (C_E) is neglected; this transition is assumed to be dominated by ROS (R) in these simulations.

TABLE 2

Table 2. Table of parameters.

When we initiate the simulation at t = 0, a population of necrotic cells occupies the impact area. We show our results as six series of snapshots in time, taken at t = 0, 1, 5, 7, 10, 14 days. We plot days 1, 7, and 14 in particular since those correspond to the attendant experiments. The time length of 14 days is used in our simulation as it is typical of such experimental set-ups.

Figure 4 shows the population densities of the three states of healthy cells. Figure 5 shows the healthy, catabolic, and EPOR-active cell population densities. Figure 6 shows the live cell population densities, i.e., the populations of all healthy and sick cell classes. Figure 7 compares the concentrations of IL-6 (our representative inflammatory cytokine) and EPO (our primary anti-inflammatory cytokine). Our focus is on the balance of pro- and anti-inflammatory cytokines. We choose IL-6 as representative since it a known inflammatory cytokine. We remark that TNF-α is also a viable choice. Figure 8 shows the concentration of DAMPs. Figure 9 shows the concentration of ROS.

FIGURE 4

Figure 4. The densities of healthy cells (C_U, C_T, C_E) at t = 0, 1, 5, 7, 10, 14 days. We show the integral of C_T over all ages to obtain a total population density. C_U is graphed with a solid line, C_T with a dash-dot line, and C_E with a dashed line.

FIGURE 5

Figure 5. The densities of healthy, catabolic and EPOR-active cells (C_U + C_T + C_E, S_T, S_A) at t = 0, 1, 5, 7, 10, 14 days. The contribution of pre-catabolic cells is the integral of C_T over all ages to obtain a total population density. We show the integral of S_T over all ages as well to obtain a total population density. Healthy cells are graphed with a solid line, catabolic cells S_T with a dash-dot line, and EPOR-active cells S_A with a dashed line.

FIGURE 6

Figure 6. The density of live cells, (C_U + C_T + C_E + S_T + S_A) at t = 0, 1, 5, 7, 10, 14 days. We show the integrals of C_T and S_T over all ages to obtain total population densities.

FIGURE 7

Figure 7. The concentrations of IL-6 (F) and EPO (P) at t = 0, 1, 5, 7, 10, 14 days. IL-6 is graphed using a solid line and EPO is graphed using a dashed line.

FIGURE 8

Figure 8. The concentration of DAMPs (M) at t = 0, 1, 5, 7, 10, 14 days.

FIGURE 9

Figure 9. The concentration of ROS (R) at t = 0, 1, 5, 7, 10, 14 days.

Simulation results show that after the initial injury, the cells adjacent to the impact area quickly sense the release of ROS and begin to switch states. At the same time, DAMPs released by necrotic cells start to trigger the pre-catabolic cell population (C_T) to become catabolic (S_T). These catabolic cells then start to produce more DAMPs, inflammatory cytokines, and ROS that aggravate the inflammation. Because of the non-zero transition times in the C_T → C_E and S_T → S_A transitions, cells cannot produce the anti-inflammatory cytokine EPO nor express a receptor for EPO immediately. The earliest an anti-inflammatory response can occur is day 1 (shown in the simulation figures). In Figures 4–6, we see that the population density of healthy cells C_U starts to decline as those cells switch to the pre-catabolic and catabolic states (C_T and S_T). By day 5, we observe a significant density of pre-catabolic (C_T), EPO producing (C_E) and catabolic (S_T) cells. From day 5 onward, there is some cell death adjacent to the impact area followed by the continuing spread of the cartilage lesion. Due to the low release rate of EPO by the C_E population, we do neither have significant conversion of C_E cells to C_U cells nor do EPOR-active cells (S_A) recover to a healthy state. By day 10, there is an appreciable density of EPO-producing cells (C_E), and EPOR-active cells (S_A) awaiting EPO to finish their conversion to the healthy state (C_U). However, the proportion of sick cells to healthy cells continues to grow during the simulation period, which was chosen to match common experimental set-ups.

In the simulation period of 14 days, the decrease of ECM density is quite small. This is because ECM degradation is a much slower process than the apoptosis of cells.

Figures 7–9 show the concentration of the chemical concentrations: IL-6, EPO, DAMPs, and ROS. The EPO concentration is too low to trigger the anti-inflammatory process efficiently, which is an expected result (Brines and Cerami, 2008). Since cells will release a small amount of ROS before they become necrotic, the concentration of ROS is not zero inside the impact area at t = 0. Although the initial burst of ROS decreases rapidly due to the fast decay rate of ROS, the production of ROS is continued by catabolic cells S_T. This results in a lower but broader concentration of ROS, as expected (Brouillette et al., 2013). DAMPs concentrations peak are early at the site of the impact.

The main competition between pro- and anti-inflammatory cytokines results in a steady increase of IL-6 across the cartilage, whereas EPO is more heavily concentrated just outside the penumbra of the inflamed region. The net result is a slowing, but not full cessation, of the spread of the inflammation.

We remark that in the simulations presented, we never reach the threshold P > P_c. This is consistent with (Brines and Cerami, 2008) that chondrocytes alone cannot produce enough EPO to stop the cartilage inflammation, so that H(P_c − P) = 1 and H(P − P_c) = 0 for all r and t. Brines and Cerami (2008) described the possibility of EPO therapy for cartilage injury, in which case P may exceed P_c. We note that the model we have formulated is set up to be modified for such therapies.

3.1.1. Parameter sensitivity analysis

We examined the approximated parameters in Table 2 for sensitivity, with the exception of α₂ and Λ, which do not affect the system in these simulations, but are included in the model for generality. Holding all other parameters to their base value in Table 2, we perturbed a given parameter. We found a region around each parameter base value where the qualitative behavior of the system did not change. The perturbed values for each parameter and an interval in which the system did not change qualitatively or appreciably quantitatively are shown in Table 3. Specifically, as each parameter was increased within its perturbation range, we found most noticeably that the following densities and concentrations decrease:

• for λ_R, EPOR-active cells (S_A) and DAMPs (M);

• for λ_M, sick cells (S_T, S_A), IL-6 (F), ROS (R), and DAMPs (M);

• for λ_F, EPOR-active cells (S_A) and DAMPs (M);

• for λ_P, normal healthy cells (C_U), pre-catabolic healthy cells (C_T), catabolic cells (S_T), IL-6 (F), EPO (P), and DAMPs (M);

• for α₁, all variables, except for EPOR-active cells (S_A) and ECM (U);

• for β₁₁, all healthy cells (C_U, C_T, C_E), ECM (U), and EPO (P);

• for β₁₂, all healthy cells (C_U, C_T, C_E) and EPO (P);

• for β₁₃, normal healthy cells (C_U) and ECM (U);

• for κ₁, healthy pre-catabolic cells (C_T), catabolic cells (S_T), EPOR-active cells (S_A), IL-6 (F), ROS (R), and DAMPs (M);

• for κ₂, healthy pre-catabolic cells (C_T), EPOR-active cells (S_A), IL-6 (F), ROS (R), and DAMPs (M);

• for ${\mathrm{\mu }}_{S_T}$, none (we see no appreciable change for the interval around the base parameter);

• for ${\mathrm{\mu }}_{S_A}$, none (we see no appreciable change for the interval around the base parameter);

• for ${\mathrm{\mu }}_{D_N}$, normal healthy cells (C_U), pre-catabolic healthy cells (C_T), and all chemicals (R, M, F, P);

TABLE 3

Table 3. Table of perturbed parameter values.

The most sensitive parameters in the system are λ_R, λ_M, and λ_F. In particular, when these parameters are relatively small compared to their base value, solutions that are otherwise monotone become non-monotone or somewhat oscillatory.

3.2. Experimental Validation

The sequence and distributions of the cell states match what is expected from past and current understanding of the “penumbra” around a lesion, in which an initial burst of pro-inflammatory cytokine expression is eventually attenuated by an increase in anti-inflammatory cytokine expression. Our major simulation result predicting this progression (Figure 7) was mainly in agreement with experimental results showing an acute post-injury increase in IL-6 expression and decrease in EPO expression, with both factors gradually returning to baseline at 14 days (Figures 10 and 11, respectively), signaling the end of the injury response cycle. It is noteworthy that our simulation does not predict a runaway pro-inflammatory response. This is also consistent with experimental findings, which showed that the initial damage done to the cartilage matrix by impact does not progress over time. Importantly, these results imply that additional factors beyond the single impact insult are required to drive further cartilage degeneration. Such factors may include synovitis or abnormal joint loading, which accompany most joint injuries.

FIGURE 10

Figure 10. The concentration of IL-6 measured experimentally. Cell counts per 0.15 cm² is the density of chondrocytes expressing IL-6 at a level detectable by immunohistochemistry.

FIGURE 11

Figure 11. The concentration of EPO measured experimentally. Cell counts per 0.15 cm² is the density of chondrocytes expressing EPO at a level detectable by immunohistochemistry.

We bring attention to some specific features of the validation result. Figure 11 shows higher concentrations of EPO in day 14 than earlier. In Figure 10, the concentration of IL-6 in day 7 is higher than the concentration in day 1 or day 14. These correspond to our simulation results in Figure 7.

We clarify that our current experimental techniques cannot extract precise concentration values for cytokines throughout the cartilage matrix. Instead, they measure the numbers of chondrocytes expressing each cytokine at a level detectable by immunohistochemistry (“positive cells”). Figure 12 shows these proteins being expressed and then released into extracellular space. Although there are limitations to counting only positive cells in this fashion, for validation purposes of this model this type of analysis is more than sufficient.

FIGURE 12

Figure 12. The concentration of EPO as a function of “positive cells.” Each scale bar is 100 μ. Immunohistochemistry images were analyzed for “positive cells,” indicating that they were above the background thresholding technique and within our cell size criteria. (A) shows concentration at day 1 at the impact site, (B) at day 7 at the impact site, (C) at day 14 at the impact site, (D) at day 1 a distance 2 mm from the impact site, (E) at day 7 at 2 mm away, and (F) at day 14 at 2 mm away. We see clearly EPO being less constitutively expressed at day 1, when compared to day 7 or 14.

The immunohistochemistry analysis indicates relative rises and declines in the amounts of cytokine present in and around individual cells where they may act in autocrine or paracrine fashion. Thus, the immunohistochemical approach gives a rough estimation of the distribution of bioavailable cytokines. In this sense, we see agreement with our predictions: the broad, nearly uniform, distributions of IL-6, peaking near day 7; and the more unimodal distributions of EPO, increasing steadily.

4. Discussion

We present an age- and space-structured model to simulate the development of an articular cartilage lesion after a single blunt impact. We simplified the model to a radially symmetric geometry. Based on previous experimental findings, we hypothesized that the consequences of mechanical trauma to cartilage depend on the balance between competing chondrolytic and chondroprotective responses of local chondrocytes. Whether damaged cartilage is stabilized, or begins down the path to progressive degeneration, has been regarded as a matter of great biologic complexity. However, our mathematical simulation confirms that the complex behavior of this system can be modeled using just two cytokines with opposing pro- and anti-inflammatory activities. The post-injury pattern of expression for IL-6 and EPO indicated by immunohistology suggests they may play roles in this binary system.

Our model did not predict runaway degeneration after a simulated impact injury. This result is appropriate and matches bench experiments showing that explanted cartilage recovers from impacts of similar magnitude, which cause minimal structural damage. Structurally damaging impacts deserve attention, but will require consideration of spatial heterogeneity in the impact site, which is beyond the scope of the current work. In addition, extrinsic factors that are sure to influence cartilage recovery in vivo, including mechanical stress and synovitis, must eventually be incorporated in the simulation to predict in vivo outcomes. Adding synovitis to our model as a second source of pro-inflammatory/pro-chondrolytic cytokines is computationally straightforward, but awaits parameterization data from in vivo models. Adding mechanical stress effects involves a much more extensive effort to combine our biomathematical model with a finite element model, which is the subject of ongoing research.

In summary, although our simulation results were obtained using a relatively small number of approximated parameters, our calculations provide useful predictions of the formation of cartilage lesions after a single blunt impact. Although we limited simulations to only 2 weeks in this work, it is possible to predict outcomes over the much longer time frames. Thus, our in silico approach may ultimately enable us to examine the long-term effects of various joint injury scenarios in people, which is difficult or impossible to address in experimental models.

Author Contributions

XW and BA developed the mathematical model and software used to solve the model equations. XW conducted the numerical simulations. MB conducted the immunohistochemistry, and slide preparation and scanning for the validation experiments. MB and JM designed the validation experiments. XW, MB, BA, and JM contributed to the writing and editing of the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Anneliese D. Heiner conducted the tissue culture experiments that generated the sample set used for model validation. Funding: BA was partially supported by the NSF under award DMS-0914514. XW, BA, MB, and JM were partially supported by NIAMS grant #1 P50 AR055533. MB was supported by a Merit Review Award from the Department of Veterans Affairs.

References

Ayati, B. P. (2000). A variable time step method for an age-dependent population model with nonlinear diffusion. SIAM J. Numer. Anal. 37, 1571–1589. doi: 10.1137/S003614299733010X

CrossRef Full Text | Google Scholar

Ayati, B. P. (2006). A structured-population model of Proteus mirabilis swarm-colony development. J. Math. Biol. 52, 93–114. doi:10.1007/s00285-005-0345-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ayati, B. P. (2007a). “Modeling and simulation of age- and space-structured biological systems,” in Current Developments in Mathematical Biology, eds K. Mahdavi, R. Culshaw, and J. Boucher (Singapore: World Scientific Publishing), 107–130.

Google Scholar

Ayati, B. P. (2007b). Modeling the role of the cell cycle in regulating Proteus mirabilis swarm-colony development. Appl. Math. Lett. 20, 913–918. doi:10.1016/j.aml.2006.11.003

CrossRef Full Text | Google Scholar

Ayati, B. P. (2009). A comparison of the dynamics of the structured cell population in virtual and experimental Proteus mirabilis swarm colonies. Appl. Numer. Math. 59, 487–494. doi:10.1016/j.apnum.2008.03.023

CrossRef Full Text | Google Scholar

Ayati, B. P. (2011). Microbial dormancy in batch cultures as a function of substrate-dependent mortality. J. Theor. Biol. 293, 34–40. doi:10.1016/j.jtbi.2011.10.003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ayati, B. P., and Dupont, T. F. (2002). Galerkin methods in age and space for a population model with nonlinear diffusion. SIAM J. Numer. Anal. 40, 1064–1076. doi:10.1137/S0036142900379679

CrossRef Full Text | Google Scholar

Ayati, B. P., and Dupont, T. F. (2005). Convergence of a step-doubling Galerkin method for parabolic problems. Math. Comput. 74, 1053–1066. doi:10.1090/S0025-5718-04-01696-5

CrossRef Full Text | Google Scholar

Ayati, B. P., and Dupont, T. F. (2009). Mollified birth in natural-age-grid Galerkin methods for age-structured biological systems. Nonlinearity 22, 1983–1995. doi:10.1088/0951-7715/22/8/012

CrossRef Full Text | Google Scholar

Ayati, B. P., and Klapper, I. (2007). A multiscale model of biofilm as a senescence-structured fluid. Multiscale Model. Simul. 6, 347–365. doi:10.1137/060669796

CrossRef Full Text | Google Scholar

Ayati, B. P., and Klapper, I. (2012). Models of microbial dormancy in biofilms and planktonic cultures. Commun. Math. Sci. 10, 493–511. doi:10.4310/CMS.2012.v10.n2.a4

CrossRef Full Text | Google Scholar

Ayati, B. P., Webb, G. F., and Anderson, A. R. A. (2006). Computational methods and results for structured multiscale models of tumor invasion. Multiscale Model. Simul. 5, 1–20. doi:10.1137/050629215

CrossRef Full Text | Google Scholar

Bianchi, M. E. (2007). DAMPs, PAMPs and alarmins: all we need to know about danger. J. Leukoc. Biol. 81, 1–5. doi:10.1189/jlb.0306164

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Brines, M., and Cerami, A. (2008). Erythropoietin-mediated tissue protection: reducing collateral damage from the primary injury response. J. Intern. Med. 264, 405–432. doi:10.1111/j.1365-2796.2008.02024.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Brouillette, M. J., Ramakrishnan, P. S., Wagner, V. M., Sauter, E. E., Journot, B. J., McKinley, T. O., et al. (2013). Strain-dependent oxidant release in articular cartilage originates from mitochondria. Biomech. Model. Mechanobiol. 3, 565–572. doi:10.1007/s10237-013-0518-8

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Eckardt, K.-U., Boutellier, U., Kurtz, A., Schopen, M., Koller, E. A., and Bauer, C. (1989). Rate of erythropoietin formation in humans in response to acute hypobaric hypoxia. J. Appl. Physiol. 66, 1785–1788.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Farndale, R. W., Sayers, C. A., and Barrett, A. J. (1982). A direct spectrophotometric microassay for sulfated glycosaminoglycans in cartilage cultures. Connect. Tissue Res. 9, 247–248. doi:10.3109/03008208209160269

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Graham, J. M. (2013). A measure of control for secondary cytokine-induced injury of articular cartilage: a computational study. Appl. Math. Comput. 223, 225–236. doi:10.1016/j.amc.2013.07.072

CrossRef Full Text | Google Scholar

Graham, J. M., Ayati, B. P., Ding, L., Ramakrishnan, P. S., and Martin, J. A. (2012). Reaction-diffusion-delay model for EPO/TNF-α interaction in articular cartilage lesion abatement. Biol. Direct 7, 9. doi:10.1186/1745-6150-7-9

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ito, T., Kawahara, K., Okamoto, K., Yamada, S., Yasuda, M., Imaizumi, H., et al. (2008). Proteolytic cleavage of high mobility group box 1 protein by thrombin-thrombomodulin complexes. Arterioscler. Thromb. Vasc. Biol. 28, 1825–1830. doi:10.1161/ATVBAHA.107.150631

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Klapper, I., Gilbert, P., Ayati, B. P., Dockery, J., and Stewart, P. S. (2007). Senescence can explain microbial persistence. Microbiology 153, 3623–3630. doi:10.1099/mic.0.2007/006734-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Leddy, H. A., and Guilak, F. (2003). Site-specific molecular diffusion in articular cartilage measured using fluorescence recovery after photobleaching. Ann. Biomed. Eng. 31, 753–760. doi:10.1114/1.1581879

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Lu, Y. C. S., Evans, C. H., and Grodzinsky, A. J. (2011). Effects of short-term glucocorticoid treatment on changes in cartilage matrix degradation and chondrocyte gene expression induced by mechanical injury and inflammatory cytokines. Arthritis Res. Ther. 13, R142. doi:10.1186/ar3456

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Martin, J. A., Ramakrishnan, P. S., Lim, T.-H., Thedens, D., and Buckwalter, J. A. (2011). “Articular cartilage and intervertebral disk,” in Orthopaedic Knowledge Update 10, ed. J. M. Flynn Vol. 3 (Rosemont, IL: American Academy of Orthopaedic Surgeons), 23–36.

Google Scholar

Shampine, L. F., and Thompson, S. (2001). Solving DDEs in MATLAB. Appl. Numer. Math. 37, 441–458. doi:10.1016/j.cmpb.2013.05.026

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Terada, C., Yoshida, A., Nasu, Y., Mori, S., Tomono, Y., Tanaka, M., et al. (2011). Gene expression and localization of high-mobility group box chromosomal protein-1 (HMGB-1) in human osteoarthritic cartilage. Acta Med. Okayama 65, 369–377.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Wang, X., Ayati, B. P., Brouillete, M. J., Graham, J. M., Ramakrishnan, P. S., and Martin, J. A. (2014). Modeling and simulation of the effects of cyclic loading on articular cartilage lesion formation. Int. J. Numer. Meth. Biomed. Eng. 30, 927–941. doi:10.1002/cnm.2636

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wedlock, D. N., Aldwell, F. E., and Buddle, B. M. (1996). Molecular cloning and characterization of tumor necrosis factor alpha (TNF-alpha) from the Australian common brushtail possum, Trichosurus vulpecula. Immunol. Cell Biol. 74, 151–158. doi:10.1038/icb.1996.20

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Zhou, S., Cui, Z., and Urban, J. P. G. (2004). Factors influencing the oxygen concentration gradient from the synovial surface of articular cartilage to the cartilage-bone interface: a modeling study. Arthritis Rheum. 50, 3915–3924. doi:10.1002/art.20675

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar