Reaction-Diffusion-Delay Model for EPO/TNF-α Interaction in articular cartilage lesion abatement

Background Injuries to articular cartilage result in the development of lesions that form on the surface of the cartilage. Such lesions are associated with articular cartilage degeneration and osteoarthritis. The typical injury response often causes collateral damage, primarily an effect of inflammation, which results in the spread of lesions beyond the region where the initial injury occurs. Results and discussion We present a minimal mathematical model based on known mechanisms to investigate the spread and abatement of such lesions. The first case corresponds to the parameter values listed in Table 1, while the second case has parameter values as in Table 2. In particular we represent the "balancing act" between pro-inflammatory and anti-inflammatory cytokines that is hypothesized to be a principal mechanism in the expansion properties of cartilage damage during the typical injury response. We present preliminary results of in vitro studies that confirm the anti-inflammatory activities of the cytokine erythropoietin (EPO). We assume that the diffusion of cytokines determine the spatial behavior of injury response and lesion expansion so that a reaction diffusion system involving chemical species and chondrocyte cell state population densities is a natural way to represent cartilage injury response. We present computational results using the mathematical model showing that our representation is successful in capturing much of the interesting spatial behavior of injury associated lesion development and abatement in articular cartilage. Further, we discuss the use of this model to study the possibility of using EPO as a therapy for reducing the amount of inflammation induced collateral damage to cartilage during the typical injury response.Table 1 Model Parameter Values for Results in Figure 5 Table of Parameter Values Corresponding to Simulations in Figure 5 Parameter Value Units Reason DR 0.1 cm2day Determined from [13] DM 0.05 cm2day Determined from [13] DF 0.05 cm2day Determined from [13] DP 0.005 cm2day Determined from [13] δR 0.01 1day Approximated δM 0.6 1day Approximated δF 0.6 1day Approximated δP 0.0087 1day Approximated δU 0.0001 1day Approximated σR 0.0001 micromolar⋅cm2day⋅cells Approximated σM 0.00001 micromolar⋅cm2day⋅cells Approximated σF 0.0001 micromolar⋅cm2day⋅cells Approximated σP 0 micromolar⋅cm2day⋅cells Case with no anti-inflammatory response Λ 10 micromolar Approximated λR 10 micromolar Approximated λM 10 micromolar Approximated λF 10 micromolar Approximated λP 10 micromolar Approximated α 0 1day Case with no anti-inflammatory response β1 100 1day Approximated Β2 50 1day Approximated γ 10 1day Approximated ν 0.5 1day Approximated μSA 1 1day Approximated μDN 0.5 1day Approximated τ1 0.5 days Taken from [5] τ2 1 days Taken from [5]Table 2 Model Parameter Values for Results in Figure 6 Table of Parameter Values Corresponding to Simulations in Figure 6 Parameter Value Units Reason DR 0.1 cm2day Determined from [13] DM 0.05 cm2day Determined from [13] DF 0.05 cm2day Determined from [13] DP 0.005 cm2day Determined from [13] δR 0.01 1day Approximated δM 0.6 1day Approximated δF 0.6 1day Approximated δP 0.0087 1day Approximated δU 0.0001 1day Approximated σR 0.0001 micromolar⋅cm2day⋅cells Approximated σM 0.00001 micromolar⋅cm2day⋅cells Approximated σF 0.0001 micromolar⋅cm2day⋅cells Approximated σP 0.001 micromolar⋅cm2day⋅cells Approximated Λ 10 micromolar Approximated λR 10 micromolar Approximated λM 10 micromolar Approximated λF 10 micromolar Approximated λP 10 micromolar Approximated α 10 1day Approximated β1 100 1day Approximated β2 50 1day Approximated γ 10 1day Approximated ν 0.5 1day Approximated μSA 1 1day Approximated μDN 0.5 1day Approximated τ1 0.5 days Taken from [5] τ2 1 days Taken from [5] Conclusions The mathematical model presented herein suggests that not only are anti-inflammatory cy-tokines, such as EPO necessary to prevent chondrocytes signaled by pro-inflammatory cytokines from entering apoptosis, they may also influence how chondrocytes respond to signaling by pro-inflammatory cytokines. Reviewers This paper has been reviewed by Yang Kuang, James Faeder and Anna Marciniak-Czochra.


Background
Articular cartilage is composed of cells known as chondrocytes. Mechanical stress and injury is known to kill (via necrosis) chondrocytes and results in the formation of lesions on the cartilage surface [1][2][3][4]. The limited capacity of chondrocytes to self-repair, together with certain aspects of the typical injury response such as inflammation, can cause the spread of lesions and development of osteoarthritis. Recent research [5] suggests that inflammatory cytokines such as tumor necrosis factor a (TNF-a) play a significant role in causing the spread of cartilage lesions. Anti-inflammatory cytokines such as erythropoietin (EPO) play an antagonistic role to TNF-a, limiting the area over which a lesion can spread by counteracting some of the effects of inflammation [5]. Moreover, the authors in [5] describe the potential use of EPO as a therapy for cartilage injury and lesion abatement The full clinical implementation of an EPO derived therapy for osteoarthritis may be a difficult task. One initial step is to establish experimentally that EPO suppresses inflammation in vitro. The potential for EPO to suppress catabolic and inflammatory responses of chondrocytes to various alarmins, i.e. molecules that trigger the innate immune response, was explored in a cell culture model (Figure 1). Bovine articular chondrocytes in Table 1 Model Parameter Values for Results in Figure 5 Table of Parameter Values Corresponding to Simulations in Figure 5 Parameter Value Units Reason

cm 2 day
Determined from [13] D M 0.05

cm 2 day
Determined from [13] D F 0.05

cm 2 day
Determined from [13] D P 0.005

cm 2 day
Determined from [13] δ R 0.01  1 day serum-free culture were pre-treated with EPO and then exposed to alarmins in the form of cell lysates or the 29 kDa fibronectin fragment (Fn-f). Interleukin-1b (IL-1b) was used as a control. Cell lysates were made by repeated freeze-thawing of excess cells from the same culture. Culture media were collected for immunoblot analysis of matrix metalloproteinase-3 (MMP-3) and A disintegrin and metalloproteinase with thrombospondin motifs 5 (ADAMTS-5). Cell layers were extracted for analysis of nuclear factor B (NF-B) and inducible nitric oxide synthase (iNOS). The results revealed that lysates, Fn-f, and IL-1b induced NF-B phosphorylation ( Figure 1A). Fn-f, and IL-1b but not lysates also induced iNOS, MMP-3 and ADAMTS-5 expression ( Figure 1B, C, D). In the case of NF-B, EPO treatment blocked activation by lysates, but not by Fn-f or IL-1b. EPO also dramatically reduced iNOS, MMP-3, and ADAMTS-5 expression in response to Fn-f and IL-1b. These results underscore the ability of EPO to counter injury effects on multiple pathways. Figure 2 illustrates the balance between pro-inflammatory (e.g. TNF-a) and anti-inflammatory (e.g. EPO) cytokines during the typical injury response. Damage due to chemical or mechanical stress initiates production of alarmins, such as damage associated molecular pattern molecules (DAMPs) [5][6][7]. These alarmins trigger production of pro-inflammatory cytokines such as TNF-a by cells near the initial injury site (see (b) in figure 2) [5]. Thus there is a penumbra ((c) in figure 2) of "sick" cells formed at the boundary of the initial injury as a result of inflammation. While still viable, these cells are at risk of dying themselves due to the presence of the pro-inflammatory cytokines, resulting in the spread of injury as shown in (d)-(e) in figure 2. However, cells far enough from the injury are able to produce antiinflammatory cytokines such as EPO to check the spread of injury and promote healing. This is shown (d*)-(e*) in figure 2. It is also suggested in [5] that TNFa acts in such a way to limit EPO signaling. Thus there is a critical balance between pro-inflammatory and antiinflammatory cytokines which determines the spreading properties of an articular cartilage lesion associated with injury. Moreover, the authors in [5] suggest that introduction of exogeneous EPO at an early time after an initial injury could be an effective therapy for minimizing the amount of collateral damage to cartilage that occurs as a result of inflammation.
The goal of this work is to study via mathematical modeling the EPO/TNF-a interaction described above, and in particular to provide a conceptual framework for a truly mechanistic understanding of some aspects of cartilage injury closely related to lesion expansion. Herein we develop a model capable of simulating the scenario illustrated in figure 2. Moreover, the model is developed as a potential for studying the use EPO or other anti-inflammatory cytokines as a therapy for lesion abatement. The type of model, and in particular the one described below, in this article is of value as it provides a link between observed phenomena and the mechanisms involved in driving that phenomena. Models such as those presented below can also help to generate hypotheses and suggest experiments or fundamental quantities to be measured. More practically, mathematical modeling such as we use in this paper could provide quick and inexpensive initial screening of potential therapies.

Methods
In this section we describe the mathematical model developed to represent the biological interactions of chondrocytes and cytokines described in [5]. We aim for a minimal model based on known mechanisms considered to be the dominant factors in articular cartilage lesion abatement. By a minimal model we mean one in which the removal of any component results in behavior that is inconsistent with the typical injury responses in cartilage as discussed in [5]. We note that the model described below can be expanded or modified to include further interactions as new results from experiment and observation dictate. In particular, further chemical pathways such as interleukin-1 b (IL1-b) can easily be incorporated if necessary. The chemical species included     below are chosen for their functionality that is their action on chondrocytes during injury response. Hence they can be replaced with any other chemical species whose effects on chondrocytes are functionally the same. The chondrocyte cell states described below represent the biological actions of the cells in response to signaling by a particular cytokine during the typical injury response. These actions are assumed to be analogous to those of other cell types, for which aspects of innate immune response, such as local inflammation considered in this paper, are well established.
We consider a population of chondrocytes fixed in matrix, which is typical of articular cartilage. Sub-populations of chondrocytes are assumed to exist in different states corresponding to the chemical signals being received by the cells during injury response. Figure 3 shows the states in which subpopulations of chondrocytes may exist. We refer to the normal state of a subpopulation of chondrocytes as the healthy state. We denote by C the population density (cells per unit area) of healthy chondrocytes at a given time and location. As a result of inflammation and injury healthy chondrocytes can enter into a "sick" state. Cells in this state are at risk death (via apoptosis) unless their signaling by TNF-a is somehow limited. We consider two subpopulations of cells in the sick state. We denote by S T the population density of cells in the "catabolic" state. Catabolic cells are chondrocytes that have been signaled by alarmins and are capable of synthesizing TNF-a and other cytokines associated with inflammation. Healthy cells signaled by DAMPs or TNF-a enter into the catabolic state and begin to synthesize TNF-a and produce reactive oxygen species (ROS). Catabolic cells that are signaled by TNF-a express a receptor (EPOR) for EPO and make up the subpopulation of sick cells we refer to as EPOR active. It should be noted that there is a time delay of 8-12 hours before a cell expresses the EPO receptor after being signaled to become EPOR active [5]. We denote the population density of EPOR active cells by S A . Since EPOR active cells express a receptor for EPO, they may switch back to the healthy state if signaled by EPO. However, as discussed in [5] TNF-a limits production of EPO. Thus there is a balance between EPO and TNF-a that determines the spreading behavior of cartilage lesions. The catabolic and EPOR active cells together make up the population of cells forming the penumbra as illustrated in figure 2. We also consider a "dead" state for subpopulations of chondrocytes. This includes necrotic cells D N and apoptotic cells D A . We note that for the purposes considered herein that apoptotic cells do not feed back into the system and thus are not explicitly represented in the mathematical model. Due to the abrupt nature of the injury, we assume that the initial injury results in necrosis of cells at the injury site. Furthermore, we assume that cell death due to secondary cytokineinduced injury is strictly through apoptosis. The reasoning here is that necrosis is a nonspecific event that occurs in cases of severe pathological cell and tissue damage, whereas secondary cytokine-induced injury corresponds with a physiologic form of cell death used to remove cells in a more orderly and regulated fashion and there is evidence that often, this is via apoptosis [8]. While these may be simplifying assumption, it is not entirely clear what type of cell death dominates in osteoarthritis, and there is some evidence to support our assumptions [8]. Necrotic cells release alarmins such as DAMPs that initiate the injury response. Either catabolic or EPOR active cells become apoptotic if signaled by inflammatory cytokines. Figure 4 details the chemical signaling and the switching of chondrocyte cell states represented in the mathematical model presented in the following. An initial injury creates a population of necrotic cells which release alarmins (such as DAMPs) [6,7]. We denote by M the concentration of DAMPs at a given time and location. The DAMPs signal healthy cells near the injury to enter the catabolic state resulting in the production of TNF-a. We denote the concentration of TNF-a at a given time and location by F. The inflammatory cytokine TNF-a has several effects on the system: It 1. feeds back to continue to switch healthy cells into the catabolic state, 2. causes catabolic cells to enter the EPOR active state [5], 3. influences apoptosis of catabolic and EPOR active cells, 4. degrades extracellular matrix (denoted by U) which results in increased concentrations of DAMPs, 5. has a limiting effect on production of EPO [5].
Catabolic cells also produce reactive oxygen species (ROS) which influences the production of EPO by healthy cells. We denote the concentration of ROS at a given time and location by R. There is a time delay of 20-24 hours before a healthy cell signaled by ROS will begin to produce EPO.
In developing a mathematical model to represent the scenario described above we assume that the chemicals diffuse throughout a domain. This diffusion is essential in determining the spatial behavior, i.e. the expansion or abatement of the lesion. Since the chondrocytes are fixed in the matrix we do not consider cell motility. First the mathematical model contains four equations (one for each chemical species) describing the dynamics (1) Next the model consists of four equations for the population densities of cells in the healthy, catabolic, EPOR active, and necrotic states. Each of these equations are of the form: change in population density of particular cell state = population density of cells switching into that state due to appropriate signaling -population density of cells switching out of that state due to appropriate signaling.

(2)
Finally there is an equation corresponding to the degradation of extracellular matrix by TNF-a which feeds back into the production of DAMPs.
The model equations for the concentrations of the chemical species are: The equation describing matrix degradation by TNF-a is: We note that the right hand side of this equation appears in (4) as part of the production of DAMPs as we have assumed that degraded matrix releases alarmins. The equations for the switching of cell states are: In (8) and (9) the function H(·) is given by The constant P c represents a critical level of EPO above which the effects of DAMPs and TNF-a on healthy cells is limited.
Based on the diagram in figure 4 and the assumptions underlying that diagram, one could assume that the function H(·) is identically one, i.e. H ≡ 1. However, if we take H ≡ 1 then figure 4 implies, and computational results confirm, that there is infinite feedback into the system by alarmins, TNF-a, and catabolic cells. As noted in [5] "the pro-inflammatory arm of the injury response is inherently self-amplifying". This does not allow for lesion abatement or limitation of secondary pro-inflammatory cytokine induced injury described in [5]. This suggests that signaling by anti-inflammatory cytokines such as EPO not only promote the switch from the EPOR active state to the healthy state but also influences the response of healthy cells to alarmins and pro-inflammatory cytokines. Thus the mathematical model suggests that in some way the anti-inflammatory cytokines limit the switch from the healthy state to the catabolic state so that there is not an infinite feedback into the system by alarmins, TNF-a, and catabolic cells. This is consistent with the observation of Brines and Cerami on the antagonistic relationship between TNF-a and EPO, that "each is capable of suppressing the biological activity of the other" [5]. In the model equations (3-11) we take the function H(·) as in (12) for convenience. However, there may be a more appropriate form for the function H(·) that must be determined through experiment to discover the nature of the effects that sufficient concentrations of EPO or other antiinflammatory cytokines have on healthy chondrocytes.
We note some features of the mathematical model (equations (3)(4)(5)(6)(7)(8)(9)(10)(11)): 1. we have incorporated the time delays for activation of the EPO receptor and synthesis of EPO (6), (9),(10), 2. it requires a concentration of TNF-a and DAMPs together for apoptosis of catabolic cells (see (9)), 3. after some time necrotic cells decay to an inert state (see (11)) so that there is not a continuous production of DAMPs for all time from the initial injury. This also corresponds to the loss of cartilage such as is sometimes associated with osteoarthritis.

Results and discussion
Using the mathematical model (equations (3-11)) we can simulate the two scenarios depicted in figure 2. That is, we can manipulate parameter values that correspond first to the case where inflammation is not limited by the presence of anti-inflammatory cytokines (d) and (e) in figure 2 and second to the case where anti-inflammatory cytokines limit the development of secondary inflammation induced injury (d*) and (e*) in figure 2. For the computations and simulations we assume a two dimensional dimensional domain. Moreover, as depicted in figure 2 we assume circular symmetry so that changes only occur in the radial direction. In order to simulate the injury response we obtain numerical approximations to the system of partial differential equations (equations (3)(4)(5)(6)(7)(8)(9)(10)(11)). This is done by first employing the scheme for discretizing the symmetric diffusion operator from appendix C of [9]. This gives a semi-discrete system of delay-differential equations which are then solved in MATLAB using the dde23 code for solving delay-differential equations. For details on the methods and software for solving delay-differential equations see [10][11][12]. One of the goals of this model is to provide a conceptual framework for a truly mechanistic understanding of certain aspects of cartilage injury closely related to lesion expansion. This requires biologists to determine parameters experimentally, but we need to know the nature of those parameters before we can go out and measure them. Since many of the parameter values have yet to be determined experimentally, we choose those values to be such that the simulations give quantitatively reasonable cell populations and chemical concentrations. We also choose the parameters so that the expansion or abatement of the lesion occurs over a time period consistent with the experimental results described in [5]. We note that diffusion of the chemical species is the mechanism which drives the expansion of the lesion in our model. The diffusion coefficients have been determined through the methods presented in [13].
First, using the mathematical model (equations (3-11)) we simulate the evolution of chondrocyte population density (per area, and assuming circular symmetry, as a function of radius) over a period of ten days after an initial cartilage injury. Figure   The results shown here correspond to a situation in which there is no production of or response to EPO. We see here that the healthy cell population is decimated. In accord with figure 4 and the corresponding mathematical model (equations (3)(4)(5)(6)(7)(8)(9)(10)(11)) this is due to the uninhibited influence of pro-inflammatory cytokines.
Following the initial injury there is a quick and steady increase in the population density of catabolic and EPOR active cells and decline in the population density of healthy cells (figures 5(b)-(i)). After ten days there has formed a thick penumbra of catabolic cells which are beginning to die off (via apoptosis, see figure 5(i)) showing that without the presence of EPO there is unlimited secondary pro-inflammatory cytokine-induced (i) t=10 days Figure 6 Healthy and Penumbra cells with EPO. This figure shows the evolution of the healthy and penumbra cell population densities over a period of 10 days. Initially there is the development of a population of penumbra cells (c) due to the time lags for EPO production and EPOR activation. However, once healthy cells at the edge of the penumbra begin to produce EPO and penumbra cells become EPOR active there is an increase in the population density of healthy cells near the injury and there ceases to be expansion of the penumbra ((d) -(i)).
injury and unabated spread of the cartilage lesion. This corresponds qualitatively well with the expected results illustrated in figure 2(d), (e). Figure 6 shows the time evolution of chondrocyte population density (per area, and assuming circular symmetry, as a function of radius) over a period of ten days after an initial cartilage injury with the antagonistic effects of EPO on TNF-a. The initial conditions are the same as previously described for the simulations without EPO. The dynamics of the chondrocyte cell states population densities are the same over the first day or so as in the case with no EPO. This is due to the time delays for cells to become EPO producing and EPOR active. However, we see by day 3 or so (see figure 6(d) and compare with 6(d)) there is a significant difference in the population dynamics compared with that in the no EPO case. We begin to see a reversion of cells in the EPOR active state back to the healthy state. While there is still some amount of secondary injury due to the inflammatory process and presence of pro-inflammatory cytokines, there is a fixed radius (see figures 6(f)-(i)) beyond which the penumbra, and hence the lesion and secondary proinflammatory cytokine induced injury, cannot expand. Thus there is lesion abatement by the action of the anit-inflammatory cytokine EPO. This corresponds qualitatively well with the expected results illustrated in figure 2(d*), (e*). Now we describe the behavior of the model under small changes to the parameter values. In order to do so we vary the parameters one at a time for a small range and observe how this effects the behavior of the healthy and penumbral cells. This exploration of the parameter space gives some indication of which parameters dominate the dynamical behavior for the model equations. For a given parameter z we vary the value over a range of [0, 2z], with the exception of l R , l M , l F , l P , a which must be strictly positive. Here we focus on those parameters associated with nonlinear terms in the model equations. Table 3 summarizes the results described in what follows. Note that below and in table 3 when we say "robust" we mean that the qualitative behavior of the simulations in unchanged as the parameter value is varied. We first consider varying the values of s R , s M , s F , s P , i.e. the production rates for each of the four chemical species. For both the production of ROS, s R , and the production of EPO, s P , decreasing these values leads to simulations which steadily approach the case with no EPO. That is, lesion expansion becomes more and more pronounced until inflammation and secondary cytokine induced injury become unchecked. This behavior is as to be expected. In contrast, increasing the values for these two parameters results in the abatement of lesion expansion at earlier times. Increasing or decreasing the production rate, s F of TNF-a has the exact opposite behavior to increasing or decreasing s R or s P . For positive values of the production rate of DAMPs, s M , the long term behavior of the simulations are stable. This is because the principal role of alarmins is to initiate inflammation and thus are most influential for early time. Setting s M = 0 results in no inflammation which

Robust Robust
This tables summarizes the results of the parameter study described in the main body of the text. When we say "robust" we mean that the qualitative behavior of the simulations in unchanged as the parameter value is varied.
is as to be expected. Next we note two interesting features of the simulations which occur when the values of a, and g are decreased.
As the values of these two parameters are decreased, we observe the formation of a penumbra of smaller and smaller radius. However, the cells forming the penumbra are much slower to respond to EPO. This makes sense for the parameter a since this is the rate at which EPOR active cells convert back to the healthy state. On the other hand, γ is the rate at which catabolic cells switch to the EPOR active state. In this case, small values of g correspond to a situation in which the penumbra is made up of predominantly catabolic cells which do not respond to EPO. We note that increasing l R , l F has the same effect as decreasing s R , s F and similar for decreasing l R , l F . Also varying l M has an expected effect since it is the saturation constant for DAMPs. Finally, the other parameters are robust as noted in table 3.

Conclusions
Based on the biological assumptions summarized in the diagram shown in figure 2 we have developed a mathematical representation of some of the principal features of the injury response in articular cartilage. In particular we have captured, qualitatively, secondary pro-inflammatory cytokine induced injury and lesion based on the balance between TNF-a and EPO as described in Brines and Cermani [5]. This modeling effort suggests that anti-inflammatory cytokines may act in such a way as to influence the response of chondrocytes to alarmins and pro-inflammatory cytokines. In particular sufficient concentrations of antiinflammatory cytokines may limit healthy chondrocytes conversion to the catabolic state. Provided dosages, dosage responses, delivery method, etc. the model (equations (3)(4)(5)(6)(7)(8)(9)(10)(11)) can be applied to study the effectiveness of the introduction of exogenous EPO as a potential therapy for limiting secondary pro-inflammatory cytokine induced injury and promoting healing in articular cartilage. We note however, that due to the time delays for cells to become EPOR active, there is little or no response to EPO at early times after an injury has occurred. Thus there is a window of time during which introduction of tissue-protective EPO derivatives such as discussed in [5] are ineffective as therapies. This is evidenced by comparing figures 5 (a)-(b) with figures 6(a)-(b) in which there is little or no difference between the case with EPO and the case without. The model presented herein also has potential applications for studying other issues related to post-traumatic stress in articular cartilage. By coupling this model with mechanical models, or including geometric features of cartilage one may be able develop a more complete theoretical description of aspects of cartilage injury of relevance in the biomedical sciences.

Reviewer Comments and Responses
Reviewer number 1: Yang Kuang The authors claim to present a minimal mathematical model based on known mechanisms to investigate the spread and abatement of such lesions. This is debatable as the model involves both diffusion and time delay yet the simulation does not address the spatial structure. The simulated dynamics is rather simple which suggests that some simple nonlinear system of ODEs may generate comparable outcomes. While modeling anything biological of this complexity will inevitably involve many simplifying assumptions and even the plausible models are ad hoc, I think the authors can present much more details on their model formulation and validation. Why can not some simple ODE models accomplish what they did here? The authors mentioned that specific functional forms appearing in their model system (equations (3)(4)(5)(6)(7)(8)(9)(10)(11)) have been chosen to capture the critical thresholds and represent the function shapes qualitatively. Later the authors refer that information as figure 2 which is hardly illustrative for thresholds or function shape. It is also important to provide meaningful parameter estimation and the source of information or the lack of it.
Response to review 1: Reviewer 1 mentions that the simulated dynamics do not address the spatial structure, but the simulations (Figures 5 and 6) show time slices of radial profiles cell populations which corresponds to the radial expansion of a lesion. Thus the authors believe that the simulations do address spatial structure. Moreover, diffusion of cytokines is the principal mechanism driving the expansion of lesions in articular cartilage injuries so that omitting diffusion ignores a principal mechanism behind one of the facets of cartilage injury response that is of primary interest. We would also like to clarify the inclusion of time delays in our model. Delays are included since, as described in the biological literature, there are time lags between when cells are signaled to become either EPOR active and/or EPO producing. By using delay terms rather than ODE's that capture the lesion abatement process phenomenologically, we maintain a clear connection between model parameters and the underlying physiology.
Reviewer 1 brings up a concern regarding the parametrization of the model. Reviewer 2 also has similar concerns so we address both simultaneously. First, in the revised manuscript we have provided a table listing the parameter values for the mathematical model used in the simulations. A small number of these parameters have been determined from literature which is cited in the revised manuscript. The remaining parameters have been approximated to give reasonable results. By this we mean cell population densities and chemical concentration levels that are consistent with those observed in other, but not necessarily injury related, cartilage experiments. Also we choose approximated parameter values that give lesion expansion and abatement occurring over a time frame as is observed experimentally. Second, one of the goals of our model is to provide a conceptual framework for the orthopaedics community to understand this phenomena from the perspective of cytokine interactions and determine what information is missing, so that we can eventually have a truly mechanistic understanding of lesion abatement. Moreover, we have attempted to develop a model that will motivate additional experiments to determine parameters, since to our knowledge many of the parameter values for this type of mechanistic model have yet to be carried out. We feel that it is of value to have these parameter values determined experimentally rather than employing statistical parameter estimation techniques. The level of additional parametrization called for in this modeling effort is consistent with other similar models in the literature, for an example see the recent article [14].

Summary
This paper presents a model of inflammation and tissue damage in cartilage. The model considers secretion and diffusion of soluble inflammatory and anti-inflammatory mediators and their effects on tissue, which is modeled as a set of spatial continuum variables. Time courses of tissue populations (healthy vs. dead or damaged cells) are determined for two different scenarios starting from the same initial damage -first with negative feedback from anti-inflammatory mediators and second with antiinflammatory mediators. The simulations show that negative feedback from anti-inflammatory mediators is necessary for recovery of damaged tissue to occur. Some experimental results from cell culture experiments are presented in the Background section and show that some effects of inflammation can be mitigated by EPO, which provides support for a key assumption of the model. Major comments 1) Missing parameters. Although the structure of the model is reasonably well described, unless I somehow missed them model parameters are not given in the paper. This should be required for publication. Parameters need to be provided along with some discussion of how they were chosen and references supporting those choices.
2) The results seem very preliminary. No attempt is made to compare the results of the simulation with experiments. Results are only presented for two scenarios and no characterization of the effect of varying model parameters is given. Have wound recovery experiments been performed, for example, in which the anti-inflammatory response is somehow manipulated -either enhanced or suppressed? Can the behavior of the model be modulated by changing various parameters that could either vary or be manipulated in actual wounds? The main point seems to be that without anti-inflammatory responses acting as a negative feedback, the inflammatory response, which has built-in positive feedbacks, will proceed unchecked. Although it is useful to a degree to demonstrate that the expected behavior does occur for a realistic model (something that has not completely been shown here because of point (1)), I think that is just the jumping off point for a reasonable model and does not by itself justify publication. As the paper stands currently, I do not think there is enough here to warrant a journal publication. Minor comments p. 2. The background starts immediately with presentation of results in Figure 1 without sufficient motivation. p. 4. "We assume chondrocytes become necrotic only at the time of an initial injury." What is the basis for this assumption? p. 5 and Figure 4. The diagram is somewhat confusing. TNF-alpha is shown as stimulating ECM, whereas in the model it stimulates degradation of ECM. Why are the DAMPs released by ECM degradation shown only as increasing apoptosis of ST cells and not also of healthy cells (or EPOR active cells). p. 6. Eq. 6. Is the mathematical form of the delay equation correct for the physical scenario? The production of EPO at the current time is assumed to be proportional to the population of healthy cells at a time τ 2 earlier (taken to be 24 hours), but what if the healthy cell population is rapidly decaying during that time? It would seem that the healthy cells that become catabolic cells should not be producing EPO. p. 8. Solution of partial differential equations for system with circular symmetry. References given are for numerical solution of ordinary diff. eqs., but the scheme for solving the PDEs is not described. A one or two sentence description and a reference should be given for the scheme used. Response to review 2: As reviewer 2 points out, the experimental results discussed in this manuscript are presented in the background without motivation. In the revised manuscript we have tried to briefly motivate these results. Reviewer 2 has asked for a justification for our assumption that chondrocytes become necrotic only at the time of initial injury. Our statement in the manuscript quoted by reviewer 2 is a little ambiguous. What we meant is that the initial injury results in cell death through necrosis, while secondary cytokine-induced injury results in cell death through apoptosis. While these may be simplifying assumptions, it is not entirely clear which type of cell death dominates in osteoarthritis. Moreover, we have provided an additional reference to lend some support to our assumptions and have clarified the quoted statement in the manuscript with a more detailed discussion.
Reviewer 2 has pointed out that Figure 4 is somewhat confusing, in that it is difficult to tell from the diagram that TNF-a degrades the extracellular matrix (ECM). Although, it is clearly stated in the caption to that figure that this is the case, and the mathematical model reflects the figure is confusing. We have modified the diagram shown as Figure 4 in the manuscript to try to make it more clear that extracellular matrix is not stimulated by TNF-a but is a "catalyst" of sorts of DAMPs due to degradation of ECM by TNF-a. DAMPs influence cells in different states in different ways. DAMPs do not cause healthy cells to enter directly into apoptosis. Healthy cells signaled by alarmins must first become "sick", i.e. catabolic, before they can enter into apoptosis.
Reviewer 2 has pointed out an error in the form of the production term for EPO. The production of EPO should be proportional to the healthy cell population at the current time, not at a τ 2 earlier time. We have modified the equations in the revised manuscript and carried out and included new simulations. The results are qualitatively the same, however in the second scenario it takes a slightly longer time for lesion abatement to occur since there are fewer (at least near the wound) healthy cells to produce EPO. We note that in the first scenario we take EPO production to be zero so that there is no different from the results in the previous version of the manuscript.
Reviewer 2 has pointed out that we neglected to describe the scheme for discretizing the partial differential equations that make up our mathematical model in carrying out numerical simulations. We have addressed this by stating the scheme used and providing a reference for this scheme.
Reviewer 3: Anna Marciniak-Czochra Provided no comments for publication.