A Validated Model of the Pro- and Anti-Inflammatory Cytokine Balancing Act in Articular Cartilage Lesion Formation

Traumatic injuries of articular cartilage result in the formation of a cartilage lesion and contribute to cartilage degeneration and the risk of osteoarthritis (OA). A better understanding of the framework for the formation of a cartilage lesion formation would be helpful in therapy development. Toward this end, we present an age and space-structured model of articular cartilage lesion formation after a single blunt impact. This model modifies the reaction-diffusion-delay models in Graham et al. (2012) (single impact) and Wang et al. (2014) (cyclic loading), focusing on the “balancing act” between pro- and anti-inflammatory cytokines. Age structure is introduced to replace the delay terms for cell transitions used in these earlier models; we find age structured models to be more flexible in representing the underlying biological system and more tractable computationally. Numerical results show a successful capture of chondrocyte behavior and chemical activities associated with the cartilage lesion after the initial injury; experimental validation of our computational results is presented. We anticipate that our in silico model of cartilage damage from a single blunt impact can be used to provide information that may not be easily obtained through in in vivo or in vitro studies.


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 3 ). 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.
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.

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.

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 τ 1 and τ 2 . 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 ).

Frontiers in Bioengineering and Biotechnology | Biomechanics
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 4 (Farndale et al., 1982).

Model equations
The equations for the chemical components of our system are The initial and boundary conditions are 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 4 . The average concentration of SO 4 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, The equation for the ECM dynamics is 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 .

www.frontiersin.org
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 where σ is the spread parameter and γ 0 is the height parameter. For a fixed γ 0 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 = 1 2 for S T → S A . The equation for the population density of healthy cells not yet signaled by ROS is The equation for the population density of healthy cells signaled by ROS and in the process of becoming catabolic is with age boundary condition 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 with initial condition C E (r, 0) = 0. The equations for the population density of sick cells in the catabolic state is with age boundary condition and initial condition S T (r, a, 0) = 0. The equation for the population density of EPOR-active sick cells is with initial condition S A (r, 0) = 0. The equation for the necrotic cell population is with initial condition Frontiers in Bioengineering and Biotechnology | Biomechanics

Numerical methods
The computational methods used to solve the age-and spacestructured differential equations were developed and analyzed in Ayati (2000Ayati ( , 2007a and Dupont (2002, 2009) and have been used effectively and efficiently in the modeling and simulation of biofilms Klapper, 2007, 2012;Klapper et al., 2007;Ayati, 2011), avascular tumor invasion , and Proteus mirabilis swarm colony development (Ayati, 2006(Ayati, , 2007b(Ayati, , 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%.

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 2 , 5% CO 2 ). 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 2 . 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.

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. www.frontiersin.org 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 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 2 . 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 2 (approximately 20-60 pixels), while the lower threshold of the positive cell size is 25-30 µm 2 (approximately 10-15 pixels).

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 2 .  τ 2 1 days Graham et al. (2012) 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 α 1 , β 11 , β 12 , β 13 , κ 1 , and κ 2 , 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 α 2 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 EPOproducing healthy cells (C E ) is neglected; this transition is assumed to be dominated by ROS (R) in these simulations. 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 EPORactive 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 antiinflammatory 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.

Frontiers in Bioengineering and Biotechnology | Biomechanics
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 www.frontiersin.org 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 EPORactive 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.

Parameter sensitivity analysis
We examined the approximated parameters in Table 2 for sensitivity, with the exception of α 2 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 α 1 , all variables, except for EPOR-active cells (S A ) and ECM (U ); • for β 11 , all healthy cells (C U , C T , C E ), ECM (U ), and EPO (P); • for β 12 , all healthy cells (C U , C T , C E ) and EPO (P); • for β 13 , normal healthy cells (C U ) and ECM (U ); • for κ 1 , healthy pre-catabolic cells (C T ), catabolic cells (S T ), EPOR-active cells (S A ), IL-6 (F ), ROS (R), and DAMPs (M ); • for κ 2 , healthy pre-catabolic cells (C T ), EPOR-active cells (S A ), IL-6 (F ), ROS (R), and DAMPs (M ); • for µ S T , none (we see no appreciable change for the interval around the base parameter); • for µ S A , none (we see no appreciable change for the interval around the base parameter); • for µ D N , normal healthy cells (C U ), pre-catabolic healthy cells (C T ), and all chemicals (R, M, F, P); 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.

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) (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.
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.
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.

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 Frontiers in Bioengineering and Biotechnology | Biomechanics

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.
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.