Modeling the Effect of Blunt Impact on Mitochondrial Dysfunction in Cartilage

Mounting evidence for the role of oxidative stress in the degeneration of articular cartilage after an injurious impact requires our modeling&simulation efforts to temporarily shift from just describing the effect of mechanical stress and inflammation on osteoarthritis (OA). The hypothesis that the injurious impact causes irreversible damage to chondrocyte mitochondria, which in turn increase their production of free radicals, affecting their energy production and their ability to rebuild the extracellular matrix, has to be modeled and the processes quantified in order to further the understanding of OA, its causes, and viable treatment options. The current article presents a calibrated model that captures the damage oxidative stress incurs on the cell viability, ATP production, and cartilage stability in a cartilage explant after a drop-tower impact. The model validates the biological hypothesis and will be used in further efforts to identify possibilities for treatment and be a part of a bigger modeling&simulation framework for the development of OA.


Introduction
Osteoarthritis (OA) is a degenerative disease characterized by thinning of the joint cartilage and is associated with disability and pain. A large portion of the cases are age-related -progressive overloading of the cartilage, as well as the reduced abilities of the body to regenerate with age, result in cartilage lesions and joint inflammation [5]. Osteoarthritis can also occur after a single forceful impact injury, in which case it is referred to as post-traumatic osteoarthritis (PTOA). We hypothesize that the biochemical processes associated with PTOA are the same as those that lead to OA, only occurring on a different time scale while age-related OA can take decades to occur, PTOA can develop in a matter of months [1,4].
There are generally two hypotheses as to the biochemical source of degeneration in OA. One is the role of joint inflammation and particularly the disruption of the balance between pro-and anti-inflammatory cytokines in the joint [14]. Our group has been leading the modeling efforts in this area with several publications [10,15,16,3,11]. However, although non-invasive anti-inflammatory therapies reduce joint pain, they have not been found to slow the development of OA [12].
Recent research efforts by the University of Iowa Department of Orthopedics and Rehabilitation have identified a different chondrocyte-centered hypothesis for the development of OA: oxidative stress and particularly the disruption of mitochondrial function as a result of joint overload [6,8,7]. The present work is the first attempt known to us to model this aspect of the underlying biochemistry after blunt impact. The article is organized as follows: Section 2 presents the laboratory experimental set-up from already published work, the hypotheses, and the resulting mathematical model; Section 3 presents our results; and Section 4 is a discussion.

Model
This section details the experimental set-up we are modeling, the biological hypotheses involved, the resulting equations, and the computational work involved in solving and fitting the model to the experimental data.

Laboratory Experiments
Our modeling efforts revolve around laboratory experiments outlined in [13] and [6]. Briefly, osteochondral explants (pieces of cartilage harvested from cattle), were subjected to high-energy blunt impacts of different magnitudes, comparable to those estimated to occur in serious joint injuries (7 J/cm 2 and 14 J/cm 2 ). The effects of impacts on cell viability within 72 h post-impact were recorded [13]. Effects of 7 J/cm 2 impacts on glucosaminoglycan (GAG) content, an indicator of cartilage stiffness and stability [6], were measured at 7 and 14 days post-impact. Metabolic activity as revealed by adenosine triphosphate (ATP) content was assessed at three different sites (impact, near impact, remote) at 24h and 48h post 7 J/cm 2 impact.

Biological Hypothesis
High energy impacts to cartilage cause local oxidative chondrocyte death [13], as well as a decline in ATP production (Figure 17.1). These effects were found to be related to excessive production of reactive oxygen species (ROS) by the mitochondrial electron transport chain, which causes damage to chondrocyte mitochondria and oxidative stress that inhibits glycolytic activity. The loss of ATP affects many cellular activities, but most notably it diminishes the production of GAGs, which undermines the stability of the cartilage matrix. Describing these processes mathematically and predicting their effects on the progression of PTOA is the subject of the present work.

Mathematical Equations
To understand the processes that occur at the impact area of the cartilage explant, we formulated a mathematical model to reflect the hypothesis outlined above. Since most of the data in [13] is relative (% of viable cells, for example), we chose to use a unitless representation for these quantities. Generally, we consider 1 to be a level of the relevant variable that is considered optimal for cartilage function. Some deviations in that assumption are outlined later in this section for the choices of parameter values. We include the following variables: The following dynamical system describes the overall cellular/biochemical processes that occur after impact. Due to impact, we assume a sudden release of ROS (captured in different initial amount of ROS depending on the impact energy). This abnormal release of ROS in the explant causes oxidative stress.
The effect of oxidative stress on the chondrocytes is transitioning of cells with normal, functional mitochondria into cells that have dysfunctional mitochondria; resulting in cell death or metabolic impairment. The difference between cells with functional mitochondria versus dysfunctional mitochondria is that the latter release twice the amount of ROS released by the former. The amount of ROS available within the cell determines the amount of ATP produced (explained in more detail later), and the amount of ATP available determines the production of GAG. We assume no cell proliferation -none was observed during the experiments, and the experiment's time frame would suggest no significant number of new cells has been added.
apoptosis due to ox. stress The function S(R) represents the effect of oxidative stress on the system. It only triggers when an excessive amount of ROS is present.
The constant s C represents the direct effect of oxidative stress, represented by ROS being above the optimal level of 1, on the mitochondrial function and viability. We choose the constant α to be greater than 1. This ensures that S(R) is smooth at R =1 and simplifies the equilibrium analysis. The function f E (x) describes the energy (ATP) production. In our model x is the ratio between available ROS and viable cells R/(M + D + ). The parameter > 0 is there to avoid division by 0.
What the function f E describes is that if the relative amount of ROS is below some optimal level R 0 or above it, then the energy production is lower than optimal. Furthermore, no available ROS (R = 0) or too much ROS (R > 2R 0 ) shuts down ATP production. This idea is presented in [6]. A plot of the function can be seen in Figure 2.

Results
This section includes mathematical analysis of the system equilibria and the computational results of the model.

Mathematical Analysis
In deterministic mathematical models, like the present one, the equivalent of statistical analysis done for experimental data or for stochastic models, is equilibrium analysis. We also analyze the local effect of parameter perturbations on the different variables through local sensitivity analysis. Let (M * , D * , R * , E * , U * ) denote an equilibrium solution to (1). The only possible equilibrium solution occurs when S(R * ) = 0. Briefly, if S(R * ) = 0, M * = D * = 0, which implies that R * = 0, but then S(R * ) = 0, a contradiction. Therefore, R * ≤ 1. Furthermore, since no new cells are produced in the scope of this model and the explant assays, M * + D * ≤ 1.

Parameter relationships
We assume that under homeostasis (undamaged cartilage), the values of cell with functional mitochondria, ROS, ATP, and GAG (M, R, E, and U respectively) remains 1, while the value for cells with dysfunctional mitochondria, D, remains 0. The only reason for changes is disruption of this equilibrium due to an impact. To ensure this equilibrium, the following relationships between parameters were assumed 1. We assumed that in the function f E that R 0 = 1/(1 + ) so as to produce the maximum amount of ATP when R = 1 and M + D = 1.

2.
We considered a level of 100% M to be optimal/normal for cartilage. This assumption requires that α M = δ R , since we seek an equilibrium R * = 1 when M * = 1.
3. With the assumptions above, in order to produce an equilibrium E * = 1 when R * = 1, we assume that The details of the equilibrium analysis are given in Appendix A. Briefly, the non-trivial equilibrium is stable and will be determined by the effect of the oxidative stress on the cell viability. In other words, we expect to reach a new homeostasis with lower levels of cell viability and appropriate levels of ROS, ATP, and GAG. No chaos is present in the system.

Numerical Results and Data Fitting
System (1) was solved using the Matlab TM function ode15s. The parameter values used for generating the results can be seen in Table 1. The data used for parametrization of our model can be seen in Table 2. We note that the data is modified. In the experimental results in [13], all explants had mean initial viability of 89%, including control. If 89% viability is normal for cartilage, we divided all the data by 89% to get the normal viability to be 100% (or 1 in the simulation calculations). In other words, the viability data was scaled. The initial conditions, in order (M (0), D(0), R(0), E(0), U (0)) for the 7 J/cm 2 impact were (1, 0, 1.0202, 1, 1), and (1, 0, 1.036, 1, 1) for the 14 J/cm 2 impact simulation. We assumed that undamaged cartilage only contains cells with functional mitochondria, that ATP, and GAG content are optimal, and the impact increased the initial amount of ROS above 1, depending on the impact's energy. We fit all parameters, besides R(0) for the 14 J/cm 2 impact using the 7 J/cm 2 data in Table 2 (cell viability, ATP, GAG). Then, using the parameters we found, we fit the initial amount of ROS after the 14 J/cm 2 impact to the cell viability data in that case. We used the Matlab particle swarm function for minimizing the error. The root mean square error for the fit to the 7 J/cm 2 data is 0.074, and to the 14 J/cm 2 data is 0.123. The results are presented in Figures  3 to 9. The total cell viability (functional plus dysfunctional mitochondria) fits well with the cellular viability presented in [13], as evident from Figures 3 and  4. The ATP simulation also fits well with the available data from [6] as seen in Figure 7. The GAG simulation also fit well with the given data ( Figure 9). Overall, the model seems to capture the biochemical dynamics of the impact site of the cartilage explant.  Table 2: Data used in parameter estimation. Standard deviation is given after the mean as ±.

Sensitivity Analysis
Let us denote by S par,var the effect of the parameter par on the variable var. Standard methods of local sensitivity analysis boil down to solving a set of differential equations with respect to S par,var , namely where S par is the vector of S par,var with respect to each variable, J is the Jacobian matrix, and F is a vector of partial derivatives of the corresponding variable with respect to the parameter of interest. The parameters we want to analyse are k S , δ D , δ R , s C , k E , λ E , k U , and λ U , as well as the initial conditions for each variable, M (0), D(0), R(0), E(0), U (0). The method is outlined in [2]. The relative local effect was measured by S par,var (t)/var(t). None of the parameters affect any of the variables locally. The initial conditions had some effect, although none of them had a local effect on U . E(0) and U (0) did not affect any of the variables. The effect of changes in M (0) on M (t) is constantly 1, on D is 0, and on R and E can be seen in Figures 10 and 11. D(0) does not affect M and D and its effect on R and E can be seen in Figures 12 and 13. Changes in R(0) affects R and E, as seen in Figures 14 and 15, respectively.

Discussion
We constructed a model of the effects of oxidative stress on the energy production and proteoglycan release of a cartilage explant after a blunt impact. The model considered the effect of the impact on the mitochondrial ROS release and the resulting disruption in ATP production, which in turn negatively affects GAG release and cartilage structure. The model's results fit well with the results of the laboratory experiments both qualitatively and quantitatively. The simulations seem to capture well the cell viability dynamics, the amount of ATP available at the impact site post injury, as well as the GAG content. Particularly, they capture the recovery in ATP production after the initial cell death and disruption. They capture, to an extent, the GAG production recovery as well.
The fact that the model's outcomes are not sensitive to the local perturbations of the variables implies that our system is stable, as suggested by the equilibrium analysis. While a lot of biological systems exhibit chaos, we expect certain outcomes in this phenomenon, in particular for a significant impact on a cartilage explant to reduce the cellular viability and disrupt the production of ATP. Our system is sensitive to the initial conditions, which is understandable. The amount of ROS gradually becomes independent of the initial burst of ROS due to the impact, and dependent on the current viability of both cell types. This effect is seen in Figures 10, 12, and 14. At the same time, the amount of initial ROS has a significant effect on the amount of cells with dysfunctional mitochondria, D, which is not surprising given that we assume that oxidative stress leads to apoptosis. Significantly less of an effect was observed in the sensitivity of cells with functional mitochondria to the initial amount of ROS, so it was assumed to be 0.
Several limitations of our modeling efforts should be addressed. Scarcity of longitudinal data for ATP production and GAG availability means the model results only vaguely support the underlying hypothesis about the effect of the post-impact oxidative stress on the biochemical functions in articular cartilage, but they do not allow us to conclude that the dynamics we describe are entirely accurate. Furthermore, the whole model is non-dimensional and the estimated parameters non-mechanistic, which makes the model only appropriate for estimating relative levels of the biochemical compounds. More data and measurements would be needed for addressing these issues.
The parameters themselves were estimated and the error found may be a local minimum, rather than a global one, so other parameter sets might give us a similar fit. However, the idea that the impact changes the initial amount of ROS released in the tissue seems to work to validate the viability of cells after the 14 J/cm 2 impact. More data at that stronger impact level would be needed to validate our ATP and GAG predictions. Our predictions for the amount of ROS present, both in the explant impact area, and per cell in the impact area seem to qualitatively capture expectations, namely high amounts of initial ROS, which level off eventually, as seen in [9].
A major result presented in [13] suggests the positive effect of antioxidants, N-acetylcysteine (NAC) particularly, on the post-impact cellular viability and overall cartilage stability, as measured by GAG content. The fact that treating the cartilage explant with NAC results in mitigating the effects of the blunt impact, leads to the conclusion that reducing oxidative stress and mitochondrial dysfunction post-impact is a viable option for preventing the development of OA. Modeling the effect of NAC and the timing of its application will be the subject of further work.
Work on creating and implementing in silico models like the one presented here may have a significant role in predicting the harmful effect of impact on cartilage explants and eventually translate to predicting post-impact patient outcomes. Using mathematical models to describe and quantify the biochemical reactions that lead to cartilage damage after an impact, may eventually remove the need to run a large portion of laboratory experiments. An ODE system such as (1) is easily encapsulated in code (e.g. Matlab) that can be used directly in the lab to predict experimental outcomes. The computations of the solutions the system take on the order of seconds on contemporary desktop equipment, which can save significant experimental time and resources relative to conducting expensive, time-consuming, and error-prone experiments. However, in the near term more experiments are needed to inform models and for creating an accurate map of the important biochemical interactions.   Figure 1: A diagram of the dynamics expressed in 1. External strain from the blunt impact causes cells with functional mitochondria to transition into cells with dysfunctional mitochondria and cells with dysfunctional mitochondria to go into apoptosis. Cells with dysfunctional mitochondria release twice the amount of reactive oxygen species (ROS) as normal cells, which further affects the oxidative stress. ROS is used in production of ATP, which in turn is utilized for the release of glycosaminoglycans (GAG), which strengthen the ECM.    Figure 7: Relative concentration of ATP after an impact of 7 J/cm 2 and its fit to available ATP data from [6].  Figure 9: Relative concentration of GAG after an impact of 7 J/cm 2 and its fit to available GAG data from [13]. for which 0 is the equilibrium (keep in mind R → 1 + is fixed). The eigenvalues of this system are −k S s C (R − 1) and −δ D s C (R − 1), both negative. Therefore, the equilibrium (M * , D * ) is asymptotically stable, which implies that the equilibrium of the original system is also asymptotically stable.