Discerning the spatio-temporal disease patterns of surgically induced OA mouse models

Osteoarthritis (OA) is the most common cause of disability in ageing societies, with no effective therapies available to date. Two preclinical models are widely used to validate novel OA interventions (MCL-MM and DMM). Our aim is to discern disease dynamics in these models to provide a clear timeline in which various pathological changes occur. OA was surgically induced in mice by destabilisation of the medial meniscus. Analysis of OA progression revealed that the intensity and duration of chondrocyte loss and cartilage lesion formation were significantly different in MCL-MM vs DMM. Firstly, apoptosis was seen prior to week two and was narrowly restricted to the weight bearing area. Four weeks post injury the magnitude of apoptosis led to a 40–60% reduction of chondrocytes in the non-calcified zone. Secondly, the progression of cell loss preceded the structural changes of the cartilage spatio-temporally. Lastly, while proteoglycan loss was similar in both models, collagen type II degradation only occurred more prominently in MCL-MM. Dynamics of chondrocyte loss and lesion formation in preclinical models has important implications for validating new therapeutic strategies. Our work could be helpful in assessing the feasibility and expected response of the DMM- and the MCL-MM models to chondrocyte mediated therapies.


Text B. Zoning
The superficial zone was chosen as the area of interest as it was observed that the naive mice had near identical chondrocyte populations through time. Immunofluorescence stainings with specific antibodies directed against collagen type II were used to identify the region of articular cartilage. The clear and consistent sparsity of collagen type II staining in the upper cartilage layer was used to outline the superficial zone of the tibial plateau and femoral condyle. The superficial zone was considered the region of interest, and was used in the following for the automated analysis of chondrocyte and apoptosis populations.

Text C. Thresholding and Contouring
First, the image was thresholded on RGB such that any noise and light scattering due to the extracellularmatrix (ECM) in the image was removed. The DAPI and TUNEL signals were processed within their respective channels only. The thresholding of the respective channel was kept fixed for all images. Then any pixel which had an intensity above the threshold was amplified to the maximum intensity. Following this procedure, the image should only contain regions where there is a possibility of a signal and regions of no signal. Then a contour was drawn around each connected set of pixels. Each contour was classified as a signal. The number of pixels within the contour was considered the area of the signal. The signal was then classified to be accepted, rejected or recounted (Fig 1).

Text D. Classification
The area of a true positive signal in the superficial zone in all images was assumed to follow the same contour area distribution. This distribution was empirically calculated by manually measuring the areas of 200 cells that are considered to be a true positive signal from randomly selected images. These sampled areas were then fitted to a gamma distribution. The accepted signal area begins within the 95 percentile confidence interval of the fitted gamma distribution. The extracted signals were then automatically assessed if they should be accepted, rejected or recounted (Fig 1). Within the processing of tissue sections it occurs that sometimes cell nuclei leave their natural position in the lacunae and are found in a region of collagen type II stained extracellular matrix. Since these nuclei cannot be clearly assigned to a specific lacunae or region, these signals were rejected and excluded from the analysis. These signals can be identified as they show an overlap of the DAPI or TUNEL staining and the collagen type II staining. In the case that separate nuclei are close to each other, the signal appears as the union of more than one nuclei. In such cases, if the signal area is greater than the 95th percentile, this signal is considered for recounting, by cutting them using the WaterShed method. All mentioned rules were applied to the images, so no human adjustments were needed to the data set generated by the automated process. To test the sensitivity and specificity of the data extraction method, ten slices were randomly selected from the pool of all images (excluding the images used for constructing the empirical distributions of the signal area) to compare the automated classification and human classification.
The automated method had a sensitivity of 91% and a specificity of 95%.

Text E. Test Change Point Analysis
The aim is to split the time course into two phases: a transient phase and a stationary phase. The transient phase captures the studied effect changing, whereas, the stationary phase captures the studied effect stabilising. We distinguish these phases by detecting changes in the mean and variance through time. Specifically, in a transient phase we expect the means between adjacent time points to vary dramatically, whereas in the stationary phase, we expect the means among adjacent time points to be very similar. There may be more than two phases which the studied effect undergoes, however, in this work we consider modelling only two phases (we seek at most one change point). Finding a change point in a time series is referred to as a binary segmentation-or At Most One Change (AMOC) analysis in the Change Point Analysis literature [2]. We now formulate the problem statistically as done in [3]: There are N time points which are ordered chronologically, and in each time point there are I replicates.
Let µ n , σ 2 n be the mean and the variance of the studied effect at time n. Then the following hypotheses are considered: versus the alternative hypothesis, where µ and σ 2 are unknown common parameters if there is no change, and k is the possible change point location. We use the Schwarz Information Criterion (SIC) as the test statistic for testing the hypotheses above.
Before presenting the computational aspects of the hypothesis test, we describe the steps and principles of the analysis. The null hypothesis states that fitting a single Gaussian to all the samples is a better model than aggregating points on either side any k and fitting two Gaussians respectively. We denote SIC(N ) to be the Schwarz Information Criterion for fitting all samples to one Gaussian (H 0 ), and SIC(k) the Schwarz Information Criterion of splitting the samples at time k and fitting two Gaussians respectively (H 1 for k). Hence, we have a series of alternative hypotheses; then to reject H 0 implies that there is no k, no splitting, which describes the samples better than a single Gaussian. More formally, the hypothesis test is given as follows: The term c α is the critical value from the distribution of SIC(N ), which was empirically calculated and given in [3, Table 1]. In our case c α = 10.317 for α = 0.05 and the total number of samples is equal to 18. If H 0 is rejected, then the positionk is a good estimate for the change point in the interval. We now present the details for calculating the test statistics described above. For k ∈ {2, . . . , N − 1} we define: where L 0 is the likelihood, that all samples are from a Gaussian distribution with meanμ 0 and varianceσ 2 0 . Similarly, the SIC for the alternative hypotheses has four unknowns for each k and is given by: where L 1 μ 1,k ,μ 2,k+1 ,σ 2 1,k ,σ 2 2,k+1 is the likelihood, that samples up to and at time t = k are from a Gaussian distribution with meanμ 1,k and varianceσ 2 1,k , and samples from time greater than k are observed from a Gaussian distribution with meanμ 2,k+1 and varianceσ 2 2,k+1 .