Macroscale mesenchymal condensation to study cytokine-driven cellular and matrix-related changes during cartilage degradation

Understanding the pathophysiological processes of cartilage degradation requires adequate model systems to develop therapeutic strategies towards osteoarthritis (OA). Although different in vitro or in vivo models have been described, further comprehensive approaches are needed to study specific disease aspects. This study aimed to combine in vitro and in silico modeling based on a tissue-engineering approach using mesenchymal condensation to mimic cytokine-induced cellular and matrix-related changes during cartilage degradation. Thus, scaffold-free cartilage-like constructs (SFCCs) were produced based on self-organization of mesenchymal stromal cells (mesenchymal condensation) and (i) characterized regarding their cellular and matrix composition or secondly (ii) treated with interleukin-1β (IL–1β) and tumor necrosis factor α (TNFα) for 3 weeks to simulate OA-related matrix degradation. In addition, an existing mathematical model based on partial differential equations was optimized and transferred to the underlying settings to simulate the distribution of IL–1β, type II collagen degradation and cell number reduction. By combining in vitro and in silico methods, we aimed to develop a valid, efficient alternative approach to examine and predict disease progression and effects of new therapeutics.


Introduction
Cartilage is a highly complex tissue that can be found in different forms at various locations in the body. The cartilage type is determined by the cell density and the composition of the extracellular matrix (ECM) resulting in function-adapted properties. Fibrous cartilage found in the menisci or intervertebral discs is cell-poor and rich of type I collagen, while elastic cartilage has a high cell density and elastic fibers. Articular cartilage, a subtype of hyaline cartilage, is organized in a unique structure consisting of different layers and distinct chondrocyte phenotypes enabling impact absorbance and pressure distribution within joints. The ECM of articular cartilage is defined by high concentrations of type II collagen and several proteoglycans, mainly aggrecan and glycosaminoglycans binding water within the tissue. Cartilage in general evolves during embryonic development through mesenchymal condensation [1].
Cartilage degradation within articular joints is a major feature of osteoarthritis (OA). OA ranks among the most common musculoskeletal disorders with no available cure except for joint replacement surgery. The progression of cartilage degradation is driven by an imbalance in anabolic and catabolic metabolic processes that normally ensure cartilage maintenance and homeostasis [2]. The inflammatory microenvironment results in increased expression of specialized matrix-degrading enzymes such as matrix-metalloproteinases (MMP) or a disintegrin and metalloproteinase with thrombospondin motifs (ADAMTS). MMPs, calcium-dependent zinccontaining endopeptidases, are involved in different physiological and pathological remodeling processes. MMP-1, MMP-3, MMP-8 and MMP-13 have been found to be increased in osteoarthritic cartilage and to be responsible for specifically degrading collagen fibers, but also proteoglycans [3].
The pathogenesis of OA has not been fully understood so far, partly due to the lack of an optimal model system which sufficiently incorporates all aspects of the disease [2,4]. A vast variety of in vivo, ex vivo, in vitro and (to some extent) in silico models for OA already exists, and they all mimic distinct features of OA pathophysiology including the degradation of ECM, inflammation and alterations in cell metabolism, viability and differentiation [5][6][7]. In vivo OA models are crucial for translational research, but mainly make use of small rodents, especially mice, although these species show significant differences in articular cartilage anatomy, loading conditions and life span [8,9]. Thus, in vitro models represent an important tool to investigate the mechanisms involved in the pathogenesis of OA and to examine possible therapeutic options. Model systems include monolayer cultures with cell lines or primary chondrocytes, co-cultures, 3D-cultures, cartilage explants from either humans or animals or even new cartilage-on-a-chip approaches [10]. However, a major challenge is the source of primary chondrocytes or explants, since healthy human cartilage samples are rare. Thus, sample collection is mainly undertaken from joint replacement surgeries. Independent of the cell source, it is a consensus that the 3D cultivation of primary chondrocytes resembles the in vivo situation more closely [11]. However, most 3D culture systems involve a scaffold to provide the cells with a predetermined structure, although their distinct effects on the cells are often not considered. Recent developments aim to establish scaffold-free tissue engineered cartilage which could serve as a promising approach for cartilage repair in vivo, but also as an excellent in vitro model, since ECM formation and degradation can be evaluated without the interference of a scaffold [12][13][14]. The increasing understanding on developmental biology resulted in improved tissue-engineering approaches using e.g. self-assembly and self-organization [12] processes that recapitulate mesenchymal condensation and chondrogenesis [13,[15][16][17][18]. Therefore, mainly mesenchymal stromal cells (MSCs) are used for scaffold-free cartilage tissue engineering owing to their benefits in terms of availability, proliferative capacity and phenotype stability when compared to primary chondrocytes. Several groups have proven the possibility to produce in vitro cartilage grown scaffold-free based on MSCs [13,15]. The main aim is to develop cartilage transplants as a therapeutic option to repair cartilage defects in joints such as the knee or ankle [15]. However, it has been also proposed that those scaffold-free cartilage analogous can be also used for experimental studies [13].
In addition to biological models, in silico modeling offers a powerful tool to bring clarity to the processes that evolve during in vitro or in vivo modeling [19,20]. Few studies have been performed so far to obtain parameters for in silico modeling, although a combination of biological data with mathematical modeling promises to accelerate translation in OA research [21][22][23]. Catt et al describe a partial differential equation (PDE) model of cartilage with a focus on the production of ECM and the growth of chondrocytes synergistically leading to tissue expansion [22]. Moreover, Kar et al investigated cartilage degradation induced by interleukin-1β (IL-1β) using a PDE model based on data from the literature together with existing experimental data to calibrate their model [23]. In contrast, Baker et al used an ordinary differential equation (ODE) system to describe the interaction of anti-inflammatory cytokines, proteinases and fibronectin in osteoarthritic cartilage [24]. Since biomechanical influences are of utmost interest, several approaches have already been used, and these include modelling arthritic cartilage under cyclic compressive loading or a combining of inflammation and biomechanics [25,26].
In this study, we tested whether scaffold-free cartilage-like constructs (SFCCs), engineered by macroscale mesenchymal condensation and biomechanical stimulation, can be used as an in vitro model to investigate cytokine-driven cellular and matrixrelated changes during cartilage degradation that occur e.g. during OA. In addition, we evaluated the feasibility to adapt an existing OA-related mathematical model to resemble the matrix degradation processes in these SFCCs (figure 1).

hMSC isolation and expansion
Human mesenchymal stromal cells (hMSCs) were obtained from human femoral bone marrow of patients undergoing hip arthroplasty (provided by the Centre for Musculoskeletal Surgery, Charité-Universitätsmedizin Berlin and distributed by the 'Tissue Harvesting' Core Facility of the Berlin Institute of Health Center for Regenerative Therapies (BCRT), Germany). All protocols were approved by the Charité-Universitätsmedizin Ethics Committee and performed according to the Helsinki Declaration (ethical approval EA1/012/13). Donor information is summarized in table 1. hMSCs were cultivated in DMEM GlutaMAX™ Medium (Thermo Fisher Scientific, MA) with 20% StemMACS™ (Miltenyi Biotech, Germany), 10% fetal calf serum (FCS) (Thermo Fisher Scientific, MA), 100 units ml −1 penicillin and 0.1 mg ml −1 streptomycin (Thermo Fisher Scientific, MA) at a temperature of 37 • C in 5% CO 2 atmosphere. hMSCs were cultured separately for each donor and characterized by flow cytometry (CD90 + , CD105 + , CD73 + , CD14 − , CD20 − , CD34 − , CD45 − , HLA-DR − ) as well as osteogenic, adipogenic and chondrogenic differentiation assays. Only cells successfully passing the characterization were further cultivated up to passages 3-5.

Generation of SFCCs and experimental setup
SFCCs were produced based on a patented protocol (patent no.: EP1550716B1) [27]. In short, around 10-15 million hMSCs from one donor (randomized) were transferred into a 3D state via centrifugation followed by maturation for 3 to 4 weeks applying biomechanical forces until reaching cartilage-like stiffness [28]. Due to the input of biomechanical forces, the SFCCs are formed by self-organization [12]. Before the experiments, SFCCs were cultured for 4 weeks in DMEM GlutaMAX ™ Medium, supplemented with 10% FCS, 100 units ml −1 penicillin, 0.1 mg ml −1 streptomycin (all from Thermo Fisher Scientific, MA) and 9.39 mg l −1 ascorbic acid (Sigma-Aldrich, MO); in the following, this is referred to as 'standard medium' . In a first step, SFCCs were cultivated under cytokine-free conditions and samples were taken on days 0, 7, 14 and 21 (n = 3 per time point). Secondly, to resemble the proinflammatory state that has been described in patients with early OA, stimulation was performed using standard medium supplemented with 50 ng ml −1 recombinant human IL-1β (2 × 10 8 IU mg −1 ) and 100 ng ml −1 recombinant human tumor necrosis factor alpha (TNFα) (2 × 10 7 IU mg −1 ; both from Immun-oTools, Germany). Samples of cytokine stimulated SFCCs (STIM) and controls (CTL) were taken on day 66 Male + 21 (n = 9 with triplicates of 3 different donors). A subset of cytokine stimulated SFCCs was further cultivated for 21 more days under cytokine free conditions to evaluate whether regeneration is possible (REG) (n = 9 with triplicates of 3 different donors). SFCCs from each donor were randomized into groups and analysis was performed in a blinded fashion by numbering samples at random.

Histological staining and immunofluorescence
Samples for histological analysis were first fixed in 4% paraformaldehyde (PFA) for 6 h and subsequently treated with 10%, 20% and then 30% sucrose solution, each for 24 h. After fixation, samples were cryo-embedded in SCEM embedding medium and cryo-sections of 8-µm thickness were prepared using cryofilms (Sectionlab, Japan). Prior to each histological and immunohistochemical staining procedure, slices were dried for 20 min at room temperature. Von Kossa and Alizarin red staining were performed as published previously [29].

Immunofluorescence staining-MMP-1, MMP1-3 and TUNEL
Immunofluorescence staining was used to quantify MMPs. First, the slides were air-dried at room temperature and then rehydrated with PBS for 10 min. Subsequently, unspecific binding sites were blocked with PBS/5% FCS for 30 min. Primary MMP-13 antibody (mouse anti-human; Invitrogen Thermo Fisher Scientific, monoclonal, MA5-14247) was diluted 1:200 in PBS/5% FCS/0.1% Tween ® 20 and primary MMP-1 antibody (mouse anti-human; Invitrogen Thermo Fisher Scientific, monoclonal, MA5-15872) was diluted 1:500 in PBS/5% FCS/0.1% Tween ® 20 and incubated according to the manufacturer's instruction for 3 h. After each incubation step, the preparation was washed 3 times with PBS/0.1% Tween ® 20. The secondary antibody (goat antimouse A546; Invitrogen Thermo Fisher Scientific, A-11003) was diluted 1:500 in PBS/5% FCS/0.1% Tween ® 20 and applied for 2 h. In the final staining step, core staining was performed using DAPI (1 µg ml −1 diluted in PBS/5% FCS/0.1% Tween ® 20) for 15 min. After air bubble-free covering with FluoroMount covering medium, microscopic evaluation was performed with the fluorescence microscope BZ-9000A (Keyence, Germany) using the DAPI and TRITC channels. Image analysis was performed using ImageJ. In order to determine the cell number, the Find Maxima tool was used for the DAPI image, whereas the area of MMP positive signals was measured using the Color Threshold tool. TUNEL staining (Sigma Aldrich, USA) was performed according to the manufacturer's instructions. The positive control was treated with desoxyribonuclease (DNase) I (0.34 Kunitz units, Qiagen Germany) for 10 min.

Histomorphometry
Histomorphometry was performed using FIJI ImageJ 1.52i [30,31]. H&E stained sections were used to analyze the cell count per area and cell distribution within the SFCCs. The cell count per tissue area (cells mm −2 ) was identified from H&E overview pictures of each SFCC with 50× magnification using a modified color deconvolution method [32] (for additional detail information see figure A1 and table A1). First, a free hand selection tool was used to define the region of interest (ROI) for the section outline representing the Total Area (Tt.A.) of the section. The Gap Area (Gp.A.) where no tissue was present was identified using the Color Threshold tool and subtracted from the Tt.A. to obtain the Total Tissue Area (Tt.T.A.). Next, the Color Deconvolution plugin of ImageJ with a vector to separate hematoxylin and eosin staining into each color layer was applied. Within the hematoxylin layer, the cell nuclei were identified by applying the Threshold tool of ImageJ based on their brightness within the layer. A binary image was created representing the nuclei of the cells within the section. Finally, cell count was performed from these binary images with a combination of the Particle Analysis tool in ImageJ and manual counting.
The cell count for the Tt.T.A. was performed identically for the experiments under normal and inflammatory conditions. A ROI for the Outer Area (Ot.A.) was identified with a manual selection tool for the experiment under normal conditions. For the stimulation experiment, ROIs were determined for the Ot.A. with a reduction of section diameter to 0.95 for X-and Y-axes, which then was further divided into an Outer Core Area (Ot.C.A.) and Inner Core Area (In.C.A.) by another reduction of the diameter by 0.5 for X-and Y-axes. Sections of 2 different levels per SFCC were analyzed respectively and the mean taken for statistical analysis. In order to analyze the immunohistochemically stained sections, the DAB coverage area was measured for type I collagen (Col-1) and type II collagen (Col-2), respectively.
Two to three pictures of 100× magnification per section were analyzed and the mean taken for statistical analysis. Tt.T.A was determined again by measuring the Tt.A. and subtracting Gp.A. which was identified with the Color Threshold tool in ImageJ. The area stained positive for Col-1 and Col-2 was defined by the Color Threshold tool also. Color Thresholds were determined for each set of staining separately.

RNA isolation, cDNA synthesis and qPCR
Total RNA was isolated from the SFCCs using the TissueRuptor II (QIAGEN, Germany) to homogenize the tissue and the RNeasy Fibrous Tissue Mini Kit (QIAGEN, Germany) was used to extract the RNA according to the provided protocols. RNA concentrations were measured via NanoDrop Fluorometer (Thermo Fisher Scientific, MA) and RNA integrity was confirmed via the 2100 Bioanalyzer (Agilent Technologies, CA). Sensiscript RT Kit (QIA-GEN, Germany) was used for cDNA synthesis with 50 ng per reaction according to the manufacturers' instructions. Primers were designed using Primer Blast (NCBI, MD) and sequence analysis of qPCR products was performed at LGC genomics (LGS genomics GmbH, Berlin, Germany) to confirm primer specify (for primer sequences see table 2). To analyze RNA expression, quantitative PCR (qPCR) was performed using the DyNAmo ColorFlash SYBR Green qPCR kit (Thermo Fisher Scientific, MA) at a Mx3000P qPCR System (Agilent Technologies, CA) with approximately 1.5 ng cDNA per 20 µl reaction and the following temperature profile: 7 min denaturation at 95 • C, 45 cycles of 5 s at 95 • C, 7 s at 57 • C and 9 s at 72 • C. Two technical replicas per sample and gene were performed. After each qPCR run, a melting curve analysis was performed to confirm primer specificity. In cases where no amplification curve reached the threshold before 45 cycles, the threshold-cycle value (Ct-value) was assumed to be 45. Gene expression data is shown as ∆Ct-value (=2 −∆Ct ) with normalization to the housekeeping gene elongation factor 1-alpha (EF1A).

RNA isolation from human cartilage
Human cartilage was collected from femoral condyles taken during total knee replacement surgeries (ethical approval EA1/012/13). Cartilage was harvested from areas as unaffected as possible and transferred to RNAlater (QIAGEN, Germany) for 1 h at 4 • C: Before cryo-conservation at −80 • C, the RNAlater was removed completely. Cartilage samples were cryo-pulverized (59012N, Biospec, Bartlesville, OK) and gently resuspended in TriFast ™ (VWR, Germany) mixed with 1-bromo-3-chloropropane (Sigma Aldrich, MO). Centrifugation was performed after 10 min of incubation for 10 min at 10 000× g. The top aqueous phase was further used for RNA isolation using the RNeasy Mini Kit (QIAGEN, Germany)

Gen symbol
Sequence of forward primer Sequence of reverse primer according to the protocols provided. cDNA synthesis and qPCR were performed as described above.

Statistical analysis
Statistical analysis was performed with the Graph-Pad Prism V.8 software. Quantitative data is shown as mean ± SEM for gene expression data and mean ± SD for all other data. Since the number of samples was small and Gaussian distribution was not assumed, we used the Kruskal-Wallis test with Dunn's multiple comparison as non-parametric statistical test for group differences. Each SFCC from one donor was assumed to be an individual replicate, and thus paired analysis was not performed. Statistical numbers and adjusted p-values are listed in supplementary results. A p-value of <0.05 was considered to be statistically significant. Image analysis was carried out blinded for treatment groups. Important numbers and adjusted p-values are stated either in the text or in the graphs. Due to the explorative character of this study, no sample size calculation was performed.

In silico model generation
To describe the temporal evolution of the spatial distribution of cellular and matrix-related processes in the in vitro model of OA, PDEs were used to interpret the experimental outcomes and to enable the identification of key steps within the progression of matrix degradation. The parameters of the in silico model were calibrated based on the histomorphometric in vitro data for distribution of type II collagen (Col-2) and cell numbers. Since each section of the SFCCs which were analyzed had a thickness of 8 µm, we transformed the dimension from area (mm 2 ) to volume (m 3 ). The geometrical shape of the SFCCs was approximated by a cylindrical form (radius r = x = 5.0 · 10 −3 m and height h = y = 4.5 · 10 −3 m) (figures 1 and 4A). SFCCs show a cylinder-like shape comparable to the disc form that has been described by others before [33,34]. This cylinder-like construct was surrounded by culture medium, so that the added IL-1β was able to diffuse into the SFCC from all surfaces. The cylindrical form of the SFCC was radially symmetric around its center. Therefore, the computations were performed on a two-dimensional rectangle instead of the full threedimensional domain (figures 1 and 4(A)). In this first modeling approach, we only focused on the effects of IL-1β, although the in vitro model was additionally stimulated with TNFα. Based on the mathematical model by Kar et al (table A2) the adapted model (table  3) was obtained as described in the following. For equation (1), it was possible to distinguish between the decrease due to the matrix formation and due to IL-1β by comparing the reduction of the cell density under non-inflammatory and inflammatory conditions in the in vitro models. Thus, we could identify the parameters p 1 (non-inflammatory conditions) and p 9 (inflammatory conditions) separately. Following Kar et al, we expected a small increase of Col-2 under non-inflammatory conditions (small value for p 2 ) which was negligible compared to the significant decrease induced by the MMPs [23]. Furthermore, we have omitted equation (3) for degCol-2. This component cannot be reliably measured and has no significant impact on other components. MMP measurements are only used for verification (figure 1), and therefore, parameter p 4 and the basal MMP activity p 6 -for which no estimates exist, have been omitted. Finally, in the IL-1β equation (5), the diffusion from outside is certainly the predominant source of IL-1β, which can also be verified e.g. by performing simulations of the model of Kar et al with and without a correspondingly simplified equation for IL-1β. The generated systems of time-dependent PDEs were solved using the state-of-the-art software package KARDOS [35,36]. Given a user-specified accuracy tolerance, the code automatically computes numerical approximations by means of a fully adaptive grid in space and time based on local error estimates with optimal computational complexity. The graphical outputs were produced using MATLAB (Version 2017b, The MathWorks, Inc., Natick, MA).  In addition, a layer-like structure was visible that led to a certain loose morphology as a result of the embedding procedure. Cell density decreased over the period of 3 weeks from a mean cell density of 1.378 ± 694 cells mm −2 at d0 to a mean cell density of 505 ± 268 cells mm −2 at d21 (figure 2(C)). Although this difference is not statistically significant, there is a clear biological relevance that should be considered. In addition, histomorphometry revealed a higher cell count in the Outer Area (Ot.A.) than in the Total Area (Tt.A.), indicating a spatial arrangement of cells within the construct (figure 2(C)). Type I collagen (Col-1) coverage was higher than type II collagen coverage (Col-2) although no difference could be identified between time points (figure 2(D)). qPCR analysis revealed the mRNA expression of the cartilage specific markers COL2A1 and ACAN, although COL1A1 expression was higher when compared to that of COL2A1 (figure 2(E)). COL10A1 expression was low in comparison with all other genes. No statistical differences were found between time points indicating a stable phenotype over time. With respect to the induction of chondrogenesis, we found a strong upregulation of COL2A1 at day 0 when compared to monolayer MSCs ( figure A3(A)). In addition, there is a more pronounced expression of COL2A1 than COL1A1 or COL10A1 ( figure A3(A)). Although significant differences do remain when comparing the results to native human cartilage ( figure A3(B)), we did not observe osteogenic induction/mineralization ( figure A4). Hence, SFCCs revealed a cartilage-like as evidenced by gene expression and histological analysis ( figure 2).

Stimulated SFCCs show cellular changes and matrix degradation
To model the inflammatory environment of OA, SFCCs were treated with IL-1β and TNFα. Both are considered to be the most important proinflammatory cytokines during the onset of OA. Cytokines were applied with concentrations of 50 ng ml −1 for IL-1β, and 100 ng ml −1 for TNFα, representing a highly aggressive proinflammatory stimulation in order to achieve maximal effects on the SFCCs. Tissue softening was observed macroscopically by volume increase of SFCCs stimulated for 3 weeks (STIM) compared to untreated controls (CTL) which was partly reversed after 3 additional weeks without cytokine treatment (REG) (figures 3(A) and A2). This observation was supported by histological H&E staining, microscopically indicating pronounced water retention, maceration of the superficial cell layer and morphological changes in cell phenotypes in the STIM group ( figure 3(B)). Immunohistological analysis showed a constant Col-1 coverage with no statistically significant differences between groups (figure 3(C)). Col-2 coverage however decreased after stimulation with IL-1β and TNFα to a mean coverage of 14.4 ± 5.5% (STIM) and moreover significantly to 8.8 ± 5.9% (REG), when compared to the untreated controls (25.3 ± 4.7%; figure 3(C)). Histomorphometry revealed statistically significant changes between groups in cell count per area for the Tt.A., Ot.A. and Outer Core Area (Ot.C.A.; figure  3(D)). Therefore, a lower cell density was observed in the STIM group compared to the CTL, while slight changes were also seen in the REG group compared to CTL. No significant differences were found in the In.C.A. (figure 3(D)). TUNEL staining showed a few TUNEL positive cells in the CTL group while more TUNEL positive cells were found in the STIM group ( figure A5). On mRNA level, COL1A1 was significantly downregulated upon cytokine stimulation in comparison with controls ( figure 4(E)). The gene expression of COL2A1 and ACAN was obviously lower in the STIM and REG group. Interestingly, COL10A1 was upregulated in the REG group (figure 3(E)). As expected, gene expression levels of the inflammatory markers IL1, IL6 and IL8 were significantly upregulated compared to untreated controls (figure 3(F)). Gene expression of TNF was not found to be different  between the experimental groups. However, IL1, IL6 and IL8 were numerically diminished in the REG compared to the STIM group although no statistical significance was detected ( figure 3(F)). MMP1 and MMP3 were significantly upregulated in the cytokinetreated group compared to the CTL group and numerically downregulated in the REG compared to the STIM group ( figure 3(G)). There were no significant differences in gene expression for MMP13 within the experimental groups. In summary, stimulation with IL-1β and TNFα leads to OA-like changes which can be observed in vivo during the early phase of the disease. Additional cultivation for 3 weeks without cytokines after stimulation did partially reverse those changes.

Refining an existing mathematical model result in a reduced PDE model of the in vitro observations
In order to focus on the main processes of cartilage degradation expected during OA onset, the in silico model was built on the assumption that under noninflammatory conditions, the chondrocytes (termed here as cells) produce ECM (here mainly Col-2) while under inflammatory conditions with the addition of IL-1β, cells release MMPs, which then degrade the ECM, and potentially go into apoptosis (see figure 1). Within the mathematical model, we focused on the specific Col-2 degrading enzymes MMP-1, -3 and -13. Table 3 shows the underlying equations derived from the pathway described above and our in vitro observations (adapted from [23]). Considering the underlying in vitro experiments and the complexity of the derived in silico model, it was necessary to derive a reduced model preserving our in vitro observations and the underlying biology.
First, we studied the stability of the homogeneous nonlinear dynamical system derived from our reduced order model by neglecting diffusive processes. It turned out that all states of equilibrium are Lyapunov stable, i.e. small perturbations of the equilibrium points stay small for all times. Next, applying a sensitivity analysis, we found out that all parameters exhibit sufficiently large sensitivities with respect to variations in cell and Col-2 concentrations. This formed the basis for our parameter calibration described in the following.
In table 4 the parameter values, initial values and the boundary conditions to solve the mathematical (PDE) system are described. Concerning the parameter values, the cell apoptosis rate p 1 was fitted by a least squares approach, utilizing measurements of the cell concentration (n = 3, d = 0, 7, 14, 21). Since the MMP-1/-3/-13 measurements have been neglected up to now, the diffusion coefficient is taken from the literature (MMP-1 [37];). The remaining parameters have been calibrated with our in vitro observations using a least squares fitting. For MMP-1/-3/-13, we assumed that there are homogeneous Neumann boundary conditions for all surfaces. Because the evolution of the cell number and Col-2 was described by ODEs, i.e. where no diffusion is present, there exist no boundary conditions for these components.
Due to the diffusion of IL-1β into the SFCC with a constant concentration of IL-1β in the surrounding media, inhomogeneous Dirichlet boundary conditions were assumed to be present on the outer surfaces. For the inner surface of the construct (center of the cartilage construct), homogeneous Neumann boundary conditions were applied ensuring the symmetry assumption on the SFCC. Col-2 decreased spatiotemporally assuming a homogeneous initial distribution (second row in figure  4(B)). Since IL-1β entered from the outside and stimulated the MMP-1/-3/-13 production, MMPs increased on the boundaries first and diffused further into cartilage (third row in figures 4(B) and (C)). Under non-inflammatory conditions, the IL-1β concentration was assumed to be zero. This also holds true for MMP-1/-3/-13 (see table 3, equations (3) and (4)). Thus, Col-2 was constant over time and cells decrease only due to matrix formation leading to a representation of the healthy state observed in vitro as well.

Modeling the in vitro findings
To predict cellular and matrix-related changes over a longer time period (5 weeks) and a 10-fold lower concentration of IL-1β stimulation (5 ng ml −1 ) -which is more comparable to the in vivo situation -we accordingly adapted the in silico model. We observed a full penetration of IL-1β and MMP-1/-3/-13, a fast decline in the cell number, but a rather slower Col-2 degradation when compared to the 3 weeks simulation which had higher IL-1β concentrations (supplementary video (available online at stacks.iop.org/BF/12/045016/mmedia) Appendix A).

MMP staining in the SFCCs verifies in silico simulations
In order to verify the mathematical model, MMP-1 and MMP-13 were stained via immunofluorescence for CTL, STIM and REG sections and analyzed for their spatial distribution and displayed as a sum of that. Therefore, MMP coverage normalized to the cell number was slightly higher in the STIM group ( figure  4(D)). The relative MMP coverage was lower in the inner core area (In.C.A.) than it was in that of the Ot.C.A. (figures 4(E) and (F)). This also matches the   results gained from the in silico model and nicely verifies our approach.

Discussion
Since multiple approaches and a variety of OA model systems are necessary to fully understand the complexity of this disease, we here report on a successful combination of in vitro and in silico modeling to simulate the main features of OA, matrix degradation and the upregulation of proinflammatory cytokines induced by stimulating scaffold-free cartilage-like constructs with IL-1β and TNFα. In most cases, cartilage explants or monolayer cell cultures are used as an in vitro model to study OA and its underlying mechanism or to test potential therapeutic options. However, there are some studies using tissue-engineered cartilage for their experiments such as those of Mohanraj et al reported similar outcomes using tissue-engineered cartilage compared to cartilage explants in a model of post-traumatic OA [10,38] However, most studies on tissue engineered cartilage still focus on the development of regenerative approaches for cartilage repair and do not consider using their tissues as in vitro models for the study of cartilage and OA pathophysiology [12,14,15,39,40]. SFCCs revealed a cartilage-like phenotype detected by gene expression and histological analysis (figure 2); however, significant differences do remain when comparing the results to native human cartilage ( figure A3(B)). Thus, the gene expression of COL1A1 and the presence of Col-1 on a protein level were still visible in all SFCCs examined. hMSCs have the potential to differentiate into different lineages including the chondrogenic, osteogenic and adipogenic lineage, although it has not been yet described that MSCs can fully form human articular cartilage in vitro or in vivo. The major challenges to achieve full phenotypic conversion are the profound induction of chondrogenesis as indicated by e.g. upregulated COL2A1 expression ( figure A3(A)) and the maintenance of the chondrogenic phenotype over a longer time period (figures 2(C)-(E)). The 3D structure of the SFCCs was gained through biomechanical loading [41,42]. However, it is clear that further improvements can be achieved by e.g. adding transforming growth factor beta (TGFβ) [13,15,41], using culture regimes that are specifically adapted to the timepoint of chondrogenesis [18,43] or using iPSCs as cell source [44,45]. Of note, we could not observe hypertrophic cells or mineralization which has been described for mesenchymal condensation in vitro [16,46]. By using a macroscale approach (constructs app. 0.5 cm diameter), we avoid the disadvantages of pellets regarding their physiologically irrelevant size, geometry, high cell number and density compared to the matrix and mechanical properties [15].
The in vitro simulation of the inflammatory environment in OA through a cytokine treatment with IL-1β at 50 ng ml −1 and TNFα at 100 ng ml −1 allowed us to observe significant changes in the SFCCs with regard to a decrease in cell density and matrix degradation (figure 3) [47,48]. Although the in vivo concentrations for IL-1β and TNFα are known to be much lower as demonstrated in synovial fluid of patients suffering from knee OA [49], cytokine concentrations of 1-100 ng ml −1 are usually applied when mimicking the proinflammatory environment in early OA in vitro. The latter situation is necessary in order to shorten the length for treatment by adapting the in vitro experiments, because the chronic degenerative joint disease OA in humans usually evolves over decades. IL-1β and TNFα do not only contribute to the upregulation of proteases but also inhibit the synthesis of ECM molecules, mainly Col-2 and aggrecan [50,51]. MMPs are one predominant group of enzymatic proteins which plays an important role in the pathogenesis of OA [3,52]. In both our in vitro studies and in the in silico model, we focused on MMP-1, −3, and −13 as main Col I and II degrading enzymes [53] since MMP-2 (gelatinase A) and MMP-9 (gelatinase B) are known to interact with Col IV (MMP-2) and are mainly involved in wound healing (both) or bone development (MMP-9) [54][55][56]. Furthermore, MMP-9 is induced by MMP-3 and −13 [57].
Mathematical models are very flexible, enabling adjustments and testing of various hypotheses in parallel. However, biological experiments are most often complex, expensive and resource-demanding, leading to a gain of less experimental data than model parameters can provide [58]. Thus, problems with model non-identifiability can occur, since model parameters cannot be estimated properly. Experimental design approaches aim to resolve these problems by identifying the data gap and proposing the required additional experimental data. The flexibility of mathematical models including experimental design can be used to structure future in vitro or in vivo models efficiently. However, the first considered mathematical model, which takes most of the prominent mechanisms of OA pathogenesis into account, has unidentifiable parameters (table A2). Its applicability is nevertheless limited, since neither enough experimental human data nor reliable literature-based parameters are available for model calibration and validation. This means that more quantitative measurements of the components at different time and spatial points are needed. The reduced model, however, does allow us to reproduce the experimental data with significantly fewer parameters (table 3). Validation concerning MMP-1/−13 measurements which were not used to calibrate the model indicates that the reduced in silico model is suited to investigate the influence of different targets for matrix degradation ( figure 4). In addition, the Supplementary video in this article shows the feasibility of the in silico model as a prediction tool for further in vitro experiments. The observed fast decline in the cell number accompanied by a slower Col-2 degradation is in accordance with a previously published study indicating a high chondrocyte death before matrix changes occur [59]. In comparison to the model described by Kar et al [23], our model still considers changes of the cell concentration in space as well as over time. A more detailed analysis of the model proposed by Kar et al revealed IL-1α as the driving force in the model, while MMP-1/−3/−13, Col-2 and degraded Col-2 only have a small influence on the remaining components. Thus, the simplified model in our study displays the same features but at the same time involves a smaller number of parameters and uncertainties. Other mathematical simulation approaches for OA cartilage degradation and cartilage focus on e.g. poroelastic models or coupling cellular phenotype and mechanics [26,60,61].
Nevertheless, there are many ways to foster the current approach. The proposed in vitro model displays only one specific part of the disease and includes only one tissue type, although OA has been described as a whole organ disease including several tissues such as the synovial membrane and the subchondral bone. Thus, we are currently working on a more complex whole joint model for both in vitro and in silico aspects. Furthermore, biomechanical loading has not yet been addressed within the first in vitro and in silico models presented here, although one's awareness of its tremendous impact has already been published [62]. Moreover, we shortened the experimental time window by using high concentrations of the proinflammatory cytokines IL-1β and TNFα in our experiments. That bridges the long-lasting cumulative effect of these cytokines over years and decades in the course of OA pathogenesis. Finally, an extension of the in silico model should include the treatment with TNFα due to its important role in OA. For a profound model validation and parameter fitting, more quantitative in vitro data will be necessary.

Conclusion
In this study, we describe a human 3D in vitro model based on SFCC (mesenchymal condensation) which simulates the main features of cartilage degradation-inflammation and upregulation of matrix degrading enzymes. Advantages of such a model include the 3D environment, the possibility for mid-throughput analysis, a sufficient sample volume (macroscale) for several analyses and a wide availability of tissue-engineered cartilage with more scaffold-free approaches emerging. With the combination of in vitro and in silico modeling, we aim to allow the immediate adaptation and modifications of different models. We strived to develop a mathematical model to refine and optimize the in vitro model and vice versa, especially regarding the development of a whole joint in vitro/in silico model.

Acknowledgments
The authors would like to thank Manuela Jakstadt for her excellent technical assistance. Bone-marrow was provided by the Tissue Harvesting Core Facility of the BIH Berlin. We also thank Kar et al for making their mathematical model available for our study. The scientific discussions with David Smith (The University of Western Australia) and Atte Eskelinen (University of Eastern Finland) were very helpful and greatly appreciated. AL, FB, AD, MP and TG are members of Berlin-Brandenburg research platform BB3R and Charité 3 R . This study was funded by the German Federal Ministry for Education and Research (BMBF) (project No. 031L0070) and partly by the Wolfgang Schulze Foundation ('ArthroMo'). The work of SR was partly funded by the Trond Mohn Foundation. The work of TG was funded by the Deutsche Forschungsgemeinschaft (353142848). AL is currently being supported by the Joachim Herz Foundation. Additionally, the project was distinguished with the Berlin-Brandenburg Award for Alternatives to Animal Testing 2019 (AL and RE as official awardees). Funding bodies did not have any role in designing the study, in collecting, analyzing and interpreting the data, in writing this manuscript, and in deciding to submit it for publication. Figure A1. Pipeline for histomorphometric analysis. (A) Histomorphometric analysis for cell count per area using a modified color deconvolution method. A 50× magnification overview image of an H&E stained section was taken. The section outline was defined using the Polygonal Selection tool in ImageJ (1). In order to subtract the background (or gap area) the Color Threshold tool was utilized (2). Next, the Color Deconvolution plugin (3) was used to obtain the hematoxylin color channel only in which the cell nuclei could be identified using a threshold for brightness (4). Particle analysis could then be performed based on the binary image attained in the previous step to obtain the cell count per area [cells mm −2 ]. (B) For further analysis, the 'Total Area' was further divided into an 'Outer Area' by scaling the region of interest (ROI) for 'Total Area' to 0.95 for X-and Y-axes. The remaining 'Core Area' was then divided into an 'Outer Core Area' and an 'Inner Core Area' by scaling the ROI for the 'Core Area' by 0.5 for the X-and Y-axes. Figure A2. Comparison of the wet weight after stimulation for 3 weeks. The increase in the wet weight indicates water retention induced by matrix degradation. Graph bar show with mean ± SEM. Statistical differences between groups were tested with the Kruskal-Wallis test and Dunn's multiple comparisons test. Inner Core Area Outer area subtracted from total area, diameter reduced by 0.5 for X-and Y-axes, outer part of the remaining area represents the outer core area, the inner part the inner core area. In.C.A.