1 Introduction

Osteoarthritis (OA) is an increasingly common disease in an aging society and is also strongly associated with past injuries and joint overload, genetic factors, diet and also obesity as a quality increasingly encountered. It is a disturbance of the balance between the degradation and regeneration processes of joint cartilage, leading to chronic inflammation of the musculoskeletal system, causing motor disability. Disease factors lead to chondrocyte apoptosis, where apoptosis is programmed cell death with distinct morphological signs. It is, therefore, a strictly regulated process that maintains homeostasis in the body, and its intensification is observed in many pathological conditions. Apoptosis of chondrocytes causes mechanical damage manifested by disruption of the integrity of the extracellular matrix. Unfortunately, the phenomenon of OA is still full of unknown, so it is essential to search for grounds and mechanisms which control OA development, where many biochemical factors, such as protein biomarkers measured in human body fluids and metabolic changes, influence the development of osteophytes [1,2,3].

However, current knowledge confirms the significant role of mechanical loading and blood vessels network during osteophytes growth [4,5,6,7], where exactly osteophytes are crucial discomfort and pain factor. Therefore, in this study, an attempt was made to identify a biomechanical parameter using the theoretical model describing the development of the disease, where earlier works had already extensively characterized the theoretical model [8, 9].

2 Mathematical modeling

The mathematical model describes bone spurs development that is formulated based on the similar theoretical descriptions of bone remodeling [10,11,12,13,14]. The character of the signals is close to the one introduced in [15, 16]. The most important effects considered in the theoretical model are shown in Diagram 1.

figure a

Diagram 1. Cause and effect diagram of osteophyte onset.

The relationships between the variables are presented below where:

Evolution of blood vessels density

$$\begin{aligned} \frac{\partial \rho _{v}\left( \textbf{x}{, }t \right) }{\partial t}{=}f_{1}{(}\rho _{v}\left( \textbf{x}{, }t \right) {, }\rho _{B}\left( \textbf{x}{, }t \right) {, }\rho _{N}\left( \textbf{x}{, }t \right) {, }U\left( \textbf{x}{, }t \right) {, }E\left( \textbf{x}{, }t \right) {)} \end{aligned}$$

Evolution of nutrients density

$$\begin{aligned} \frac{\partial \rho _{N}\left( \textbf{x}{, }t \right) }{\partial t}{=}f_{2}{(}\rho _{v}\left( \textbf{x}{, }t \right) {, }\rho _{B}\left( \textbf{x}{, }t \right) {, }\rho _{N}\left( \textbf{x}{, }t \right) {)} \end{aligned}$$

Evolution of bone cells density

$$\begin{aligned} \frac{\partial \rho _{B}\left( \textbf{x}{, }t \right) }{\partial t}{=}f_{3}{(}\rho _{B}\left( \textbf{x}{, }t \right) {, }\rho _{N}\left( \textbf{x}{, }t \right) {)} \end{aligned}$$

Evolution of Young modulus

$$\begin{aligned} \frac{\partial E\left( \textbf{x}{, }t \right) }{\partial t}{=}f_{4}{( }\rho _{B}\left( \textbf{x}{, }t \right) {, }U\left( \textbf{x}{, }t \right) {)} \end{aligned}$$

The equations of nutrient and bone cells density as well as Young modulus evolution were not changed in this work. Their specific description was done in previous paper [8].

figure b

Diagram 2. Relations between the variables in the theoretical model of OA.

The system of equations is complex and mathematically challenging to solve (Diagram 2). The weight parameters stabilize the process of the calculations, but simultaneously they deepen the qualitative character of the model. In order to get the model closer to the quantitative results, an effort was made to identify a first parameter related to the mechanical transduction of the cartilage cells. The equation of evolution of blood vessels density (Eq. 1) considers an extra-mechanical signal of nutrients demand \(S{(}x{, }t{)}\) (marked as I).

Other parts of the equation represent a biomechanical signal defined as a difference between nutrient density required by the cells to survive, current nutrient density (marked as II) and the density of blood vessels already present in the surrounding area (marked as III) and were described in [8, 9].

Evolution of blood vessels density

$$\begin{aligned} \frac{\partial \rho _{v}\left( \textbf{x}{, }t \right) }{\partial t}{=}A_{1} \underbrace{S{(x,t)}}_{I} \underbrace{\int _{\varOmega } H(\Delta N(\xi ,t))e^{\frac{R(\xi ,x)}{\beta }} d \xi _{1} \xi _{2} \xi _{3}}_{II} \underbrace{\int _\varOmega \rho {v}(\xi ,t)e^{\frac{R(\xi ,x)}{\gamma }} d \xi _{1} \xi _{2} \xi _{3}}_{III} \end{aligned}$$
(1)

Where:

\(\rho _{v}\left( \textbf{x}{, }t \right) \)

blood vessels density

\(A_{{1}}\)

amplification parameter

\(S\left( \textbf{x}{,}t \right) \)

extra-mechanical signal of nutrients demand

\(H{(\Delta }N\left( \xi {, }t \right) {)}\)

Heavisade’s function

\({\Delta }N\left( \xi {, }t \right) \)

physiological signal of nutrients demand

\(R\left( \xi {, }x \right) \)

distance between arbitrary point x in the domain\(\Omega \) and the location of the signal source

\(\beta {, }\gamma \)

range parameter

In the previous paper, an additional mechanical signal S(x, t) was generated by mechanically overloaded cartilage cells [8]. Essentially, the signal was based on the assumption that the physiological mechanical load could be increased by a load approaching twice the physiological signal. However, it should be discussed or explained why the physiological mechanical signal should be doubled. Presumably the multiplication of the signal depends on a lot of biophysical factors, but it should be related to the loading at least [17]. Additionally, the signal should be calculated over the domain with the decreasing range over the distance. Therefore, a new form of the signal was proposed in this paper (Eq. 2):

$$\begin{aligned} S\left( \textbf{x}{,}t \right) {=}P_{M}{(}\textbf{x}{,}t{)}\int _\Omega {A_{{2}}{U{(}\textbf{x}{,}t{)e}}^{{-}\frac{R\left( \xi {, }x \right) }{\lambda }}d\xi _{{1}}\xi _{{2}}\xi _{{3}}} \end{aligned}$$
(2)

Where:

\(P_{M}{(}\textbf{x}{,}t{)}\)

dimensionless microstructure parameter

\(A_{{2}}\)

amplification parameter

\(U{(}\textbf{x}{,}t{)}\)

current elastic strain density

\(\lambda \)

range parameter

To identify the parameter A2, that amplifies the current elastic strain density, it was assumed that the parameter is strictly related to the applied pressure and viability of the chondrocytes. It should be remembered that hypertrophic (dying) chondrocytes release essential factors of nutrients demand and the blood vessels respond to these factors and enable osteophytes onset.

Hence, the nutrient demand factor should be amplified as chondrocytes’ viability decrease. The parameter A2 thereupon was assumed to have a linear behavior (\(y{=-}ax{+}b)\) obtained from in vitro tests (Graph 1). It was assumed that variable x in the function is related to the value of applied pressure [MPa] wherein the process determines the range of amplification of the current elastic strain density influence on a whole signal S(x, t):

For \(ax{\ge }b\), the evolution of blood vessels density will be negative or equal to 0, which means that the cells are dying immediately and there is no angiogenesis.

For \(ax{<}b\), the evolution of blood vessels density will be positive but cannot be higher than b, which means that it will tend to the amplification value intended for slight loading. Based on the literature, the temporal compression is greater than 1 MPa and goes up to 20 MPa without killing the chondrocytes [18, 19]. Nevertheless, the loading changes during dynamic gait. In this model, we assumed that continuous pressure should not go above 1 MPa.

The applied pressure was tuned to impact on the cells but not to crush them. These cells were secured by the corn-starch-based material for cell culture described in Sect. 2.

3 Material and methods

3.1 Materials

The Maize starch (CAS 9005-25-8) and glycerin (purity \(\ge \) 99%, CAS 56-81-5) from Sigma-Aldrich were used as the feedstock material to develop the starch film.

The powder material for 3D printing was Polyamide 12 (PA2200) (modulus of elasticity 1650 MPa, tensile strength 48 MPa, deformation at break 18%, melting point 176 \(^{{\circ }}\)C, density 930 kg/m\(^{{3}})\) purchased from EOS Store. The 3D printing process was carried out using a Formiga P100 printer, with the following standard parameters of the laser process, respectively: beam speed of 1200 mm/s and laser power of 16W.

Rat chondrocytes (cryo-preserved chondrocytes, chondrocyte medium kit (basal medium, fetal bovine serum, chondrocyte growth supplement, penicillin/streptomycin solution), poly-L-lysine and Dulbecco’s phosphate-buffered saline) were purchased from Innoprot and were prepared as recommended.

3.2 Methods

3.2.1 Quasi-3D base surface preparation

To mimic the mechanism of cartilage overload, a biocompatible material based on corn-starch hydrogel was prepared, and the preparation of the starch material was described in detail in our previous works [20,21,22].

The film solution formulation is shown in Table 1. The mixture of starch and glycerol was dissolved in water and heated up to 95 °C while continuously mixing with a mechanical stirrer. After cooling, the solution was cast on Teflon®plates using a K Paint applicator (RK Print, Royston, UK) equipped with a micrometer gauge set to 3 mm to ensure a constant film thickness. The films were dried for seven days in standard conditions (23 °C, 50% RH). Next, they were carefully peeled off from the Teflon®plates by hand. Obtained film sheets have dimensions of approx. 100 mm \(\times \) 150 mm. 

Table 1 Hydrogel solution formulation
Full size table
Fig. 1
figure 1

Image of the developed base surface using KEYENCE VHX 950 optical microscope at magnification 500\(\times \)

Full size image

A digital microscope (Keyence VHX-7000, Japan) was used to observe the surface of the developed material. The roughness profile surface of the base material as a candidate scaffold for cell culture was determined based on 3D microscopic images (see Fig. 1). The Sa = 1.6 µm roughness parameter defined the arithmetic average profile deviation of the surface and confirmed the possibility of using a starch-based hydrogel as a potential scaffold material for chondrocyte culture.

The corn-starch-based material was primary developed as a non-permeability film. After 4 h to precisely adhere to the bottom of the wells of sliders (vessel for cell culture) dipped in medium with cells, the structure rippled and formed a quasi-3D structure. The rippled structure with the chondrocytes culture mimics closely a cartilage tissue [23, 24].

3.2.2 In vitro tests preparation

The cells used for the experiments are rat chondrocytes (MDF) after four population doublings. These cells were cryopreserved and frozen in a 1 ml probe containing at least \({5\times }{{10}}^{{5}}{ }\)living cells. First, the culture was established in a \(\phi \) 60 mm in Petri dish, which was covered with an aqueous solution of poly-L-lysine (1 mg/ml) for 16 h and rinsed twice with sterile water before cells were seeded.

Then, the cells were thawed by immersing the probe in a water bath, and then the content of the test liquid was pipetted onto a Petri dish, adding 3 ml of previously prepared warm medium.

To ensure an even distribution of cells on the surface, the Petri dish was gently rocked. The cells were incubated overnight in a Thermoscience incubator.

After 24 h, medium was added, and cells were incubated for a further 48 h. The medium was changed until cells covered 75–80% of dish, and then the medium was changed daily until 90% of the dish was covered with cells, after which they were subcultured.

EDT trypsin solution, trypsin neutralization solution, chondrocyte medium and DPBS were used for passage, where all necessary fluids were heated to 37 °C. Twenty-four proper slides were prepared with corn-starch hydrogel material on the bottom for the new culture, which (as before in a Petri dish) were covered with a membrane on the bottom, and then the cells were rinsed with DPBS. Then, enough trypsin and DPBS solution (2 ml trypsin and 8 ml DPBS) were added to cover the entire surface. The cells were incubated for 2 min, and then the solution was pipetted into a 15 ml Falcon, leaving the Petri dish in the incubator for another minute, where at the end of trypsinization the side of the dish was gently taped to separate the remaining cells. For neutralization, about 2 ml of trypsin was added and transferred to the Falcon probe. Cells were centrifuged for 5 min at 1000 rpm, resuspended in fresh medium and counted using trypan blue. To count cells were added equal parts of trypan blue and cells suspension to 1 ml probe and incubated for 2–3 min. The cells were then pipetted into the counting chambers and counted in an EVE automatic cell counter. 10⌃4 cell suspension was used per well for slide culture, and fresh medium was added to 2 ml. Cells were cultured for two days before performing static loading tests.

After series of mechanical loading, the samples viability of cartilage cells was measured by using trypan blue and Brucker’s chamber. The cells were trypsinysied, centrifuged and resuspended in fresh medium. 200 µl sample was taken from each well and stained with 200 µl of trypan blue in separate probe. The staining took 3 min in incubator. Then, 20 µl of stained cells was transferred to Brucker’s chamber (10 µl to each chamber) and counted on inverted microscope. We have counted cells in four squares in corners and then viability—number of dead and living cells were calculated.

3.2.3 Development of a stand for cell cultures mechanical stimulation

For the purpose of mechanical stimulation, a stand for cells was developed using selective laser sintering (SLS) additive manufacture. The computer-aided design (CAD) model of the apparatus is shown in Fig. 2. The stand was printed by a 3D printer with standard process parameters laser beam speed 1200 mm/s and laser power 16 W. In SLS, laser beam moves in horizontal directions (XY plane) and cures paths on subsequent layers of material [22,23,24,25]. The stand was printed from PA2200 material, due to its biocompatibility, availability, chemical resistance, significant durability and stability of the decomposition time of the material [25]. A stand with chondrocyte culture after seven days of culture, prepared for mechanical testing, is shown in Fig. 3.

Fig. 2
figure 2

Stand for laboratory tests enables simultaneous mechanical loading of four samples, a the model of a stand in an axonometric projection, b the schematic construction of a stand with component parts (dimensions: \(75\times 100\times 110\) mm), where F is a force applied to the arm, and Ft is a calculated force acting on the press

Full size image
Fig. 3
figure 3

Four holes slider with chondrocytes culture after 7 days of culturing prepared for mechanical tests

Full size image

The investigations for each sample with certain loading were conducted 6–8 times. The difference depended on the quality of the culture—some of them were contaminated.

The mechanical loading Ft force was applied to the presses and afterward calculated as a pressure p, considering the stamp area at the end of the press. The area of the every stamp was 60 mm\(^{{2}}\). Values of pressure p for every series of mechanical loads are presented in Table 2. The stamp was cut out from a culture flask and fit to the bottom of the press. The material of the culture flask was flat and biocompatible.

Table 2 In vitro tests results data
Full size table

3.2.4 Numerical simulations

The system of mathematical equations was implemented using the COMSOL Multiphysics software [26].

Figure 4 shows the example of the initial-boundary conditions, and in Fig. 5, finite element meshes used during numerical calculations are presented. To find the finest mesh density, the mesh study was conducted. The geometry model mimics part of the joint related to the end of the bone with cartilage under the local mechanical loading p.

Fig. 4
figure 4

Initial-boundary conditions. The unit of length is mm, and since the geometry is symmetrical along the short left edge, only the results shown in the red rectangle are presented

Full size image
Fig. 5
figure 5

Representation of a mesh with 2973 triangular elements and 47904 numbers of degrees of freedom

Full size image
Table 3 Values of parameters and initial values of variables
Full size table

Table 3 presents only parameters whose values were changed as well as changed initial values of variables compared to previous work [8].

4 Results and discussion

4.1 In vitro tests

The in vitro tests were conducted on a stand dedicated to mechanical stressing of the cell culture. Cells were sided on the corn-starch-based hydrogel which likely absorbed part of the stress. The cell cultures were not directly under a pressure. It can be said that the hydrogel takes a role of extra cellular matrix (ECM) in a real-life cartilage tissue. After series of mechanical loading, the samples viability of cartilage cells was measured. Obtained results are shown in Fig. 6. It can be observed that the cell viability decreases along with increase in pressure which is much more realistic. The equation of the linear trend line is in form \(y =-221,24x {+60,207}\) and is given by a parameter \(A_{{2}}\).

Fig. 6
figure 6

Viability of cartilage cells after 5 min loading

Full size image

Even though the results align with expectations and in our opinion sufficient for theoretical deliberate, we are willful to admit that there should be more samples from a biological point of view.

4.2 Numerical tests

The proposed model with a new form of A2 parameter has been tested through numerical simulations. Figure 7 shows the distribution of von Mises stress in the first time step, and around the maximum value, there is the area of an expected pathological changes.

Fig. 7
figure 7

Distribution of von Mises stresses (a) and distribution of elastic strain energy density (b) in the first time step

Full size image

The results presented in Table 4 show the changes of bone degeneration in the expected region. The degeneration of tissue begins around where the elastic strain energy reaches its greatest value. Calculated variables that determine the bone development, such as blood vessel density, bone cell density and Young’s modulus, increase over time.

Table 4 Results of numerical analysis of blood vessel density distribution (kg/m⌃3), bone cell density distribution (kg/m⌃3) and Young’s modulus distribution (Pa) in bone—cartilage area at three different points in time
Full size table

The qualitative theoretical model reflects the general character of bone formation based on real-world medical data. However, the results, on the one hand, provide information of the predicted trend of osteophyte development but also does not provide quantitative information on the rate of tissue changes. It should be noted that, since the interaction between the biomechanical variables—osteophytes onset was investigated, the result’s scale is irrelevant.

5 Conclusions

An in vitro experiment revealed the information about the character of the relation between mechanical loading and viability of cartilage cells. As the value of mechanical loading gets higher, the cells die more. According to our assumptions, the dying cells release more signals about nutrient demand, leading to angiogenesis and osteophytes formation. This information was adapted to the theoretical model of osteophytes development during OA. Although the character of the results is as expected, the tilt angle could be slightly different. The difference could be related to used hydrogel material as a quasi-3D scaffold which was soft and degraded in the medium. The degradation could have caused some interference during automatic counting.

The numerical results show the correct character of bone spurs evolution. Three variables, which define bone tissue, evolve in the predicted direction of increased von Mises stress (Fig. 5a) and elastic strain energy density (Fig. 5b) from the bone tissue domain.

Table 4 presents three time steps of variables blood vessels density, bone cells density and Young modulus. The first time step is the initial time step, and the other two were liberally selected. The evolution of the variables is changing with time from the domain representing healthy bone tissue. The observed changes are qualitative and more correlated to the real-life observations; however, there still is a requirement to identify the rest of the amplification parameters. Only then the theoretical model will be fully applicable for medical use.