Modeling the effect of blunt impact on mitochondrial function in cartilage: implications for development of osteoarthritis

Laboratory experiments

Our modeling efforts revolve around laboratory experiments outlined in Martin et al. (2009) and Coleman, Buckwalter & Martin (2015). Briefly, osteochondral explants (pieces of articular cartilage with subchondral bone underneath, harvested from cattle), are secured at the bottom of a drop tower and subjected to high-energy blunt impacts (from a 5 mm diameter brass rod dropped onto the explant from different heights) of different magnitudes, comparable to those estimated to occur in serious joint injuries (7 J/cm² and 14 J/cm²). The effects of impacts on cell viability within 72 h post-impact were recorded by putting the explants into media containing calcein acetoxymethylester (which dyes live cells green) and ethidium homodimer-2 (which dyes dead cells red), taking images with confocal microscopy and analyzing the images from six areas of three explants (eighteen images total per time point) (Martin et al., 2009). The viability was recorder as % live cells (green) to total cells (green plus red).

Cartilage is a tissue that comprises of mainly extracellular matrix with a dispersion of chondrocytes throughout. The extracellular matrix is composed of water, collagen, and proteoglycans, in different proportions depending on the cartilage depth (more proteoglycans toward the bottom and less toward the top, which makes it heterogeneous in stiffness—stiffer as one transitions from top to bottom). The effects of 7 J/cm² impacts on proteoglycans was measured by glycosaminoglycan (GAG) assay with dimethyl methylene blue. Proteoglycans contains GAG molecules, so GAG content is an indicator of cartilage stiffness and stability (Coleman, Buckwalter & Martin, 2015). The relative GAG content (GAG at impact versus GAG at nearby undamaged site) was measured at 7 and 14 days post-impact and averaged over six explants. Metabolic activity and energy metabolism, as revealed by adenosine triphosphate (ATP) content, was assessed at three different sites (impact, near impact, remote) at 24 h and 48 h post 7 J/cm² impact. No ATP or GAG data was collected for the 14 J/cm² impact. In this model we only consider the results at the impact site because that is consistent with the viability and GAG data available.

Biological hypothesis

Cartilage is a hypoxic tissue—while most human cells’ mitochondria use an oxidant-rich process of phosphorylation to produce ATP, chondrocytes’ mitochondria produce ATP through glycolisis, which requires small amounts of oxygen. Therefore, in articular cartilage there is a fine balance between the oxidants needed for normal cellular metabolism and an excess that can cause metabolic damage. High energy impacts to cartilage cause local oxidative chondrocyte death (Martin et al., 2009), as well as a decline in ATP production (Figure 17.1 in Coleman, Buckwalter & Martin, 2015). 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 resulting loss of ATP affects many cellular activities, but most notably it diminishes the production of GAGs, which undermines the stability of the cartilage matrix. The subject of the present work is describing these processes mathematically and creating a mechanistic model that qualitatively and quantitatively describes their effects on the progression of PTOA.

Mathematical model

Based on the biological hypothesis and the available data, we constructed a system of differential equations that describes the interactions between the mechanical impact from the drop tower, the mitochondrial function of the chondrocytes in the impact area, and the resulting concentrations of ROS, ATP, and GAG. The model is unitless to reflect the available data (% of cell viability for example).

We include the following variables:

1. M(t): proportion of live cells with functional (normal, undamaged) mitochondria in the cartilage explant impact area.

2. D(t): proportion of live cells with dysfunctional (damaged, abnormal) mitochondria in the cartilage explant impact area. These cells are characterized by their release of double the amount of ROS as cells with functional mitochondria.

3. R(t): relative concentration of reactive oxygen species (ROS) in the cartilage explant impact area.

4. E(t): relative concentration of ATP (produced through glycolisis) in the cartilage explant impact area.

5. U(t): relative concentration of GAG (measure of cartilage stability) in the cartilage explant impact area.

“Relative” in the descriptions above refers to comparison with normal, undamaged cartilage explant (control in the experiments). A brief description of the models dynamical system of equations below: the blunt impact causes a burst of ROS, which creates an environment of oxidative stress. The stress causes the mitochondria of normal chondrocytes to become dysfunctional, which makes them release additional ROS, and can cause cellular apoptosis. The amount of ROS that is available controls the ATP formation in the cell if the concentration of available oxidants is too little or too high, the ATP production is decreased (Coleman, Buckwalter & Martin, 2015). Decreased ATP production negatively affects GAG production and hence concentration in the impact area. A diagram of the interactions is presented in Fig. 1. 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. We consider 1 to be a level of the relevant variable that is considered optimal for cartilage function. This idea translates to assuming that when 100% (or a fraction of 1 of the total number of cells) of the viable cells have functional mitochondria, the released amount of ROS is 1, which translates to an optimal ATP production (assumed to be 1), and overall cartilage integrity (relative GAG concentration of 1). In other words the control case is assumed to be the case where M(t) = 1, D(t) = 0, R(t) = 1, E(t) = 1, and U(t) = 1. The control case is trivial, hence not pictured in the figures that accompany the model. (1)${\displaystyle \begin{array}{cc}\frac{dM}{dt}=\hfill & \underset{\text{mito damage due to ox. stress}}{\underbrace{-k_{S}MS\left(R\right)}}\hfill \\ \frac{dD}{dt}=\hfill & \underset{\text{mito damage due to ox. stress}}{\underbrace{k_{S}MS\left(R\right)}}-\underset{\text{apoptosis due to ox. stress}}{\underbrace{{\delta }_{D}DS\left(R\right)}}\hfill \\ \frac{dR}{dt}=\hfill & \underset{\text{mito ROS release}}{\underbrace{{\alpha }_M\left(M+k_{D}D\right)}}-\underset{\text{ROS clearance}}{\underbrace{{\delta }_{R}R}}\hfill \\ \frac{dE}{dt}=\hfill & \underset{\text{ATP production}}{\underbrace{f_E\left(\frac{R}{M+D+ϵ}\right)}}-\underset{\text{utilization}}{\underbrace{{\delta }_{E}E}}\hfill \\ \frac{dU}{dt}=\hfill & \underset{\text{GAG through ATP}}{\underbrace{k_{U}U\left(1-\frac{1+{\lambda }_U}{1+{\lambda }_{U}E}U\right)}}.\hfill \end{array}}$ The function S(R) represents the effect of oxidative stress on the system. It only triggers when an excessive amount of ROS is present. (2)${\displaystyle S\left(R\right)=\left\{\begin{array}{c}0\text{if}R\le 1\hfill \\ s_C{\left(R-1\right)}^{\alpha }\text{if}R>1.\hfill \end{array}\right.}$

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. (3)${\displaystyle f_E\left(x\right)=\left\{\begin{array}{c}0\text{if}x\le 0\text{or}x\ge 2R_0\hfill \\ \frac{k_E}{{\left(x-R_0\right)}^2+{\lambda }_E}-\frac{k_E}{R_0^2+{\lambda }_E}\text{if}x\in \left(0,2R_0\right).\hfill \end{array}\right.}$

What the function f_E describes is that if the relative amount of ROS is below some optimal level R₀ or above it, then the energy production is lower than optimal. Furthermore, no available ROS (R = 0) or too much ROS (R > 2R₀) shuts down ATP production. This idea is presented in Coleman, Buckwalter & Martin (2015). A plot of the function can be seen in Fig. 2.

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₀ = 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 ${\delta }_E=\frac{k_E}{{\lambda }_E\left(1+{\lambda }_E\right)}$.

Numerical solutions and data fitting

System Eq. (1) was solved using the MATLAB 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 Martin et al. (2009), 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. 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² impact, using the 7 J/cm² 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² impact to the cell viability data in that case. We used the MATLAB particle swarm function for fitting the parameters.

Local 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 ${\displaystyle \frac{d{\overrightarrow{S}}_{\text{par}}}{dt}=J\cdot {\overrightarrow{S}}_{\text{par}}+F,}$ where ${\overrightarrow{S}}_{\text{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 Atherton, Schainker & Ducot (1975). The relative local effect was measured by S_par,var(t)∕var(t).

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.

The details of the equilibrium analysis are given in Appendix. 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

The initial conditions (cell viability and chemical concentrations at time = 0 h), in order (M(0), D(0), R(0), E(0), U(0)) for the 7 J/cm² impact were (1, 0, 1.0202, 1, 1), and (1, 0, 1.036, 1, 1) for the 14 J/cm² impact simulation. The root mean square error for the fit to the 7 J/cm² data is 0.074, and to the 14 J/cm² data is 0.123. The results are presented in Figs. 3–7. The total cell viability (functional plus dysfunctional mitochondria) fits well with the cellular viability presented in Martin et al. (2009), as evident from Figs. 3 and 4. The ATP simulation also fits well with the available data from Coleman, Buckwalter & Martin (2015) as seen in Fig. 6. The GAG simulation also fit well with the given data (Fig. 7). Overall, the model seems to capture the biochemical dynamics of the impact site of the cartilage explant.

Sensitivity analysis results

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. The values of 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, which implies that the relative value of M(t) is aways dependent on M(0). The effect of M(0) on D is 0, and on R and E can be seen in Figs. 8 and 9. D(0) does not affect M and D and its effect on R and E can be seen in Figs. 8 and 9. Changes in R(0) affects R and E, as seen in Figs. 8 and 9, respectively.