Mathematical modeling and application of IL-1β/TNF signaling pathway in regulating chondrocyte apoptosis

Introduction: Mathematical model can be used to model complex biological processes, and have shown potential in describing apoptosis in chondrocytes.

Method: In order to investigate the regulatory mechanisms of TNF signaling pathway in regulating chondrocyte apoptosis, a fractional-order differential equation model is proposed to describe the dynamic behavior and mutual interaction of apoptosis-related genes under the activation of TNF signaling pathway. Compared with the traditional molecular biology techniques, the proposed mathematical modeling has advantages to providing a more comprehensive understanding of the regulatory mechanisms of TNF signaling pathway in chondrocyte apoptosis.

Result: In this paper, differentially expressed genes induced by IL-1β in human chondrocyte apoptosis are screened using high-throughput sequencing. It is found that they were significantly enriched in the TNF signaling pathway. Therefore, a mathematical model of the TNF signaling pathway is built. Using real-time PCR experiments, mRNA data is measured and used to identify the model parameters, as well as the correlation coefficient. Finally, the sensitivity of the model parameters is discussed by using numerical simulation methods, which can be used to predict the effects of different interventions and explore the optimal intervention strategies for regulating chondrocyte apoptosis.

Discussion: Therefore, fractional-order differential equation modeling plays an important role in understanding the regulatory mechanisms of TNF signaling pathway in chondrocyte apoptosis and its potential clinical applications.

1 Introduction

Osteoarthritis (OA) is a chronic degenerative joint disease characterized by articular cartilage destruction, subchondral bone remodeling, osteophyte formation, and inflammatory changes in periarticular tissues, which seriously affects the quality of life of patients. According to statistics, the incidence of osteoarthritis is 50% in people over 50 years old, 80% in people over 60 years old, and the disability rate is as high as 53% (Wang et al., 2019). Osteoarthritis is the first chronic disease that causes disability in adults, and it is also one of the most common diseases in the world (Frank et al., 2015; Yuan et al., 2016). With the aggravation of population aging in China, the prevalence of OA will show an increasing trend, so the research on OA is becoming more and more urgent. In the process of OA, the apoptosis of chondrocytes is the main pathological feature of OA. Studies have found that the content of IL-1β in OA patients is significantly increased. When the content of IL-1β is increased, it will activate the TNF signaling pathway, and TNF-α and IL-1β in the TNF signaling pathway can promote each other to absorb and degrade articular cartilage, aggravating the degree of OA. The higher the content of IL-1β, the more severe the OA, and there is a positive correlation between them.

However, for the treatment of OA, the current goal is mainly to relieve pain and improve joint function, and there is no specific drug. Therefore, exploring an intervention that can effectively delay the progress of OA and improve the OA condition in the early stage is considered to be a very potential treatment strategy that can effectively reduce the development process of OA patients (Yuan et al., 2016). Studies have shown that chondrocyte apoptosis plays a key role in the occurrence and development of osteoarthritis, and its number is positively correlated with the severity of osteoarthritis (Gu et al., 2016; Toshihisa, 2016). Therefore, preventing and reducing chondrocyte apoptosis is an effective method for the treatment of osteoarthritis (Zahoor et al., 2017).

The signal transduction pathway of chondrocyte apoptosis is very complex, and the signal pathways involved mainly include TNF signaling pathway (Wang et al., 2005; Qiu et al., 2007). MAPK signaling pathway (Yoon and Kim, 2004; Ma and Ren, 2012; Gao et al., 2014; Cai and Li, 2019; Liu et al., 2020), JAKS/STAT signaling pathway (Li et al., 2016; Liu et al., 2016), Wnt/β-catenin signaling pathway (Hu et al., 2019; Kalamakis et al., 2019; Liu et al., 2020) and NF-κB signaling pathway (Rigoglou and Papavassiliou, 2013; Ji et al., 2017; Zhao and Gong, 2019) are the hot topics at home and abroad. Tumor necrosis factor-alpha (TNF-α) in the TNF signaling pathway is a cytokine with pleiotropic biological effects. The biological effects of TNF are triggered by two TNF receptors (TNFR) on the cell surface, and its signal transduction pathways mainly include caspase family-mediated apoptosis, adaptor protein TRAF-mediated activation of transcription factor NF-κB and JNK protein kinase (Wang et al., 2005). When TNF-α level is increased, JNK signaling pathway is activated and starts to participate in chondrocyte apoptosis, which also leads to a significant decrease in the expression of apoptosis inhibitor protein Bcl-2 (Yoon and Kim, 2004). In addition, JNK is also involved in inhibiting the expression of Sox-9 and blocking the apoptosis of chondrocytes induced by NO, while SP600125, an inhibitor of JNK, can significantly inhibit the pathological damage of cartilage (Gao et al., 2014), providing another direction for the treatment of OA.

In the IL-1β-induced primary human chondrocytes, baicalin downregulated the mRNA and protein expression of MMP1 and MMP13, and promoted the expression of collagen type II and Aggrecan. Baicalin significantly inhibits cartilage degradation in DMM-induced OA mice, suggesting that breviscapine may be a potential drug for OA treatment (Liu et al., 2020). In addition, breviscapine inhibited the migration of β-catenin and the phosphorylation of p38 into the nucleus, which is associated with the regulation of MAPK signaling pathway.

Studies have found that (Hu et al., 2019), OA chondrocytes were transfected with miR-320c and its inhibitor β-catenin-siRNA, and the results suggested that MiR-320c was decreased and β-catenin was increased in the late stage of OA chondrocytes formation. Overexpression of miR-320c and knockdown of β-catenin had similar effects on cartilage specific gene expression and hypertrophy related gene expression in OA chondrocytes. Injection of mmu-miR-320-3p attenuated OA progression in a mouse model of OA, indicating that miR-320c inhibits osteoarthritis chondrocytes degeneration by inhibiting the canonical Wnt signaling pathway, and miR-320c has the potential as a new drug for osteoarthritis treatment. Hu et al. (2020) reported that pretreatment with Loureirin A significantly inhibited IL-1β-induced production of NO, PGE2, COX-2, TNF-α, iNOS, and IL-6 in mouse articular chondrocytes. In addition, Loureirin A significantly inhibited IL-1β-mediated AKT phosphorylation and NF-κB entry into the nucleus in chondrocytes in signaling pathway studies; therefore, Loureirin A may be a potential therapeutic candidate for OA. Wu et al. (2018) found that after the activation of Notch signaling pathway by the specific activator Jagged1 protein in rat knee joint, Bax protein was upregulated and Bcl-2 protein was downregulated through the apoptotic pathway, thereby promoting chondrocyte apoptosis and aggravating OA. After the Notch signaling pathway was inhibited by the injection of γ-secretase inhibitor DAPT (GSI-IX), the Bax protein was downregulated and the Bcl-2 protein was upregulated through the apoptotic pathway, thereby inhibiting the apoptosis of chondrocytes and alleviating the development of OA.

In summary, the signaling pathways involved in the regulation of articular chondrocyte apoptosis are diverse and extremely complex, but the most significant signaling pathway regulating chondrocyte apoptosis is still unclear, and the most critical signaling molecules in the signaling pathway are also unclear. Therefore, it is urgent to screen the signaling pathways that most significantly affect chondrocyte apoptosis and discover the most critical signaling molecules in this signaling pathway. The mathematical modeling of signaling pathways through high-throughput sequencing is the key to solve the above problems.

Dynamic modeling and analysis of biological systems can simplify the biological system into a mathematical model for analysis and numerical simulation, thus replacing the actual complex, long-term, expensive and even impossible experiments, greatly improving the research efficiency, and studying the influence of artificially imposed control conditions on the operation process of biological systems. For example, infectious disease modeling (Sun and Liu, 2017; Sun and Liu, 2018; Sun et al., 2020) and signaling pathway modeling (Liu et al., 2017), etc. Kinetic modeling based on signaling pathways by relevant scholars mainly includes:

In 2003, Matsuno et al. (2003), based on the Notch signaling pathway to regulate the formation mechanism of Drosophila large intestine boundary cells, used Hybrid functional Petri net (HFPN) for modeling. In 2009, Manu et al. (2009) used gene regulatory networks to elucidate interactions between four target genes in the early embryonic period. The equilibrium point, stability and branching of the system were analyzed by using the dynamics theory of ordinary differential equations. Based on the qualitative characteristics of the dynamic system, the complex mechanism of the channel formation and pattern formation of Drosophila germ layer was further understood. Immediately following this, Nakamasu et al. (2009) proposed a system of diffusion differential equations with three factors. This model takes into account cell spreading and computationally simulates models of the wild-type pigment production mode as well as other models of mutant production modes.

In the above study, although Manu et al. (2009) calculated the equilibrium point, stability and Hopf branching, they only studied the model theoretically and did not combine the experimental phenomenon with the theoretical results. Although some progress has been made in the system of diffusion differential equations with three factors proposed by Nakamasu et al. (2009), it does not consider that model organism development is susceptible to the influence of external environment such as temperature, especially the interference of random factors. In addition, fluctuations in the level of gene products can generate noise at the molecular level, which not only affects the accuracy of the signal gradient but also reduces the output of the target. Fully considering the above reasons and current situation, the applicant established a deterministic model of Drosophila large intestine border cell formation relying on Delta-Notch signaling pathway (Liu et al., 2017), and added white noise to study the effect of random factors on Drosophila large intestine border cell formation. The model is as follows:

$\left\{\begin{array}{l}\frac{d{D}_i}{dt}=\frac{\lambda }{1+\Delta \cdot A_i}-d_1{D}_i-{\displaystyle \sum _{NG\left(\mathrm{i}\right)}}{f}_1\cdot D_i,1\le i\le NC,\\ \frac{d{N}_i}{dt}={\lambda }_N-d_2{N}_i+{\displaystyle \sum _{j\in NG\left(i\right)}}{f}_2\cdot D_j-\frac{a{N}_i}{b{D}_i+N_i},\\ \frac{d{A}_i}{dt}=-d_3{A}_i+\frac{a{N}_i}{b{D}_i+N_i}.\end{array}\right.(1)$

Differential equation dynamics theory was used to study the equilibrium point and its stability when the Delta gene was not expressed and overexpressed. The different phenotypes of the deterministic system are calculated by means of numerical simulations, which are in good agreement with the experimental results.

In addition to the Notch signaling pathway mentioned above, in 2018, Frederik et al. established a corresponding mathematical model based on the mechanism of WNT signaling pathway regulating adult hippocampal nerves, and revealed age-related changes in the nervous system of adult hippocampus (Ziebell et al., 2018; Alahdal et al., 2021). In 2019, Georgios et al. used mathematical kinetic modeling to study that inflammatory signaling and WNT antagonist sFRP5 induce quiescent neural stem cells to regulate the maintenance and regeneration ability of stem cells in the aging brain (Kalamakis et al., 2019).

However, there are few studies on the mathematical modeling of chondrocyte apoptosis based on signal pathways. Accurate mathematical models can quantitatively analyze the mechanism of each signal molecule in the signaling pathway regulating chondrocyte apoptosis at the system level. Therefore, it is particularly urgent to screen the signal pathways that most significantly affect chondrocyte apoptosis, and to establish a kinetic model of chondrocyte apoptosis regulated by this signal pathway from the perspective of mathematical biology.

2 High-throughput sequencing and enriched signaling pathway

Knee Human chondrocyte was collected from the Shenzhen Second People’s Hospital (The First Affiliated Hospital of Shenzhen University). Written informed consent was obtained from all patients. The study has been approved by the Shenzhen Second People’s Hospital (The First Affiliated Hospital of Shenzhen University), China. The method to deal with the cartilage in order to obtain the proper chondrocytes is as follows:

After scraping the cartilage tissue with a blade, rinse it three times with physiological saline containing 3% PS. Use a two-step digestion method to separate the articular chondrocytes. Transfer the cartilage tissue to a sterile culture dish and wash it three times with PBS containing 3% PS. Carefully remove any attached soft tissue or synovial fluid, and use ophthalmic scissors to crush it into pieces smaller than 1mm3. Transfer the crushed tissue to a new 50 mL centrifuge tube and add 5 times the volume of 0.25% trypsin. Digest for 20 min and then wash the remaining tissue three times with low glucose DMEM medium. Add 5 times the volume of a solution containing 0.2% type II collagenase and digest gently on a shaking bed at 37°C for 4-5 times, 25 min each time. After each digestion, collect the cells in the supernatant and culture the cells using low glucose DMEM medium containing 20% fetal bovine serum.

In the early stage, normal chondrocytes were stimulated with IL-1β at a concentration of 10 ng/mL. The results showed that chondrocytes began to apoptosis at 6h, reached the peak at 18h, and then the degree of apoptosis decreased (Figure 1).

FIGURE 1

FIGURE 1. The mRNA expression levels of Caspase-3, BCL-2 and Bax were detected by real-time PCR.

Therefore, in this study, the experiments are divided into three groups: 1) normal chondrocytes served as blank control, 2) normal chondrocytes were stimulated with IL-1β (10 ng/mL) for 6 h, and 3) normal chondrocytes were stimulated with IL-1β(10 ng/mL) for 18 h. The cells were collected and subjected to high-throughput sequencing.

The results of sequencing showed that the significant difference value was p = 3.21e-09 for TNF signaling pathway, p = 5.76e-05 for NF-κβ signaling pathway, p = 1.23e-04 for Cytokine-cytokine receptor interaction signaling pathway. The significant difference value for T-cell receptor signaling was p = 1.77e-04. The detailed results are shown in Table 1.

TABLE 1

TABLE 1. High-throughput sequencing Significance table of each signaling pathway affecting chondrocyte apoptosis.

The results of the above high-throughput sequencing suggest that the most significant signaling pathway regulating chondrocyte apoptosis under IL-1β-stimulated conditions is the TNF signaling pathway. Therefore, this project proposes a scientific hypothesis: Through the mathematical modeling of IL-1β/TNF signaling pathway screened by high-throughput sequencing, the regulatory effect of signal molecules in the IL-1β/TNF signaling pathway on chondrocyte apoptosis can be systematically analyzed, and new targets for the treatment of osteoarthritis can be further explored through the sensitivity analysis of parameters.

3 Mathematical model

According to the signaling molecular transduction mechanism of TNF signaling pathway regulating chondrocyte apoptosis, we use differential equations to model the dynamics of chondrocyte apoptosis. Firstly, the signal molecule transduction map of TNF signaling pathway regulating chondrocyte apoptosis was drawn by Portable Pathway Builder Tool 2.0 (as shown in Figure 2).

FIGURE 2

FIGURE 2. Signal molecule transduction diagram of TNF signaling pathway regulating chondrocyte apoptosis.

Here, $x_1\left(t\right),x_2\left(t\right),...,x_n\left(t\right)$ is utilized to denote the relative concentrations of mRNA for each of the signaling molecule at time $t$. The differential equation governing changes in signal molecule concentration is expressed as follows.

$\frac{d{x}_i\left(t\right)}{dt}=v_{i,production}-v_{i,\mathrm{deg}radation}(2)$

where $v_{i,production}$ and $v_{i,\mathrm{deg}radation}$ denote the production and degradation rates of the ith signaling molecule, respectively. The IL-1β/TNF signaling pathway can be modeled as a system of differential equations, with each signaling molecule’s dynamic changes represented by its own equation. In this paper, the variables $x_1\left(t\right)$, $x_2\left(t\right)$, $x_3\left(t\right)$, $x_4\left(t\right)$, $x_5\left(t\right)$, $x_6\left(t\right)$, $x_7\left(t\right)$, $x_8\left(t\right)$, $x_9\left(t\right)$, $x_{10}\left(t\right)$, $x_{11}\left(t\right)$, $x_{12}\left(t\right)$, $x_{13}\left(t\right)$, $x_{14}\left(t\right)$, $x_{15}\left(t\right)$, $x_{16}\left(t\right)$, $x_{17}\left(t\right)$, $x_{18}\left(t\right)$ and $x_{19}\left(t\right)$ represents the relative concentration of mRNA expressed by the signaling molecule IL-1β, IRAK1/4, TNFα, TNFR1, TRADD, TRAF2/5, RIP1, NIK, TAK1, IKKs, MKK4/7, FADD, JNK1/2, NF-κB, ITCH, Caspase10, c-FLIP, Caspase8 and Caspase3/7 at time $t$, respectively. Parameter $a_1$ is used to denote the growth rate of IL-1β expressed mRNA and $d_1$ the degradation rate of IL-1β. Therefore, the differential equation about the rate of mRNA change of IL-1β at time $t$ is given as $\frac{d{x}_1\left(t\right)}{dt}=a_1-d_1{x}_1\left(t\right)$. Similarly, parameter $a_2$ denotes the growth rate of IRAK1/4 mRNA expression, while $b_2$ represents the activation rate of IL-1β on IRAK1/4. Furthermore, considering the time delay effect of IL-1β on IRAK1/4 activation, a delay term denoted by $\tau$ is incorporated in this study. Additionally, $d_2$ signifies the degradation rate of IRAK1/4. Therefore, the differential equation governing the change in mRNA level of IRAK1/4 at given time $t$ is expressed as: $\frac{d{x}_2\left(t\right)}{dt}=a_2+b_2{x}_1\left(t-\tau \right)x_2\left(t-\tau \right)-d_2{x}_2\left(t\right)$. $a_3$ and $d_3$ represent the growth rate and degradation rate of TNFα, respectively. $a_4$ and $d_4$ represent the growth rate and degradation rate of TNFR1, $b_4$ represents the activation rate of TNFα on TNFR1. $a_5$ and $d_5$ represent the growth rate and degradation rate of TRADD, $b_5$ represents the activation rate of TNFR1 on TRADD. $a_6$ and $d_6$ represent the growth rate and degradation rate of TRAF2/5, $b_6$ represents the activation rate of TRADD on TRAF2/5. $a_7$ and $d_7$ represent the growth rate and degradation rate of RIP1, $b_7$ represents the activation rate of TRAF2/5 on RIP1. $a_8$ and $d_8$ represent the growth rate and degradation rate of NIK, $b_8$ represents the activation rate of TRAF2/5 on NIK. $a_9$ and $d_9$ represent the growth rate and degradation rate of TAK1, $b_9$ and $c_9$ represents the activation rate of IRAK1/4 and TRAF2/5 on TAK1. $a_{10}$ and $d_{10}$ represent the growth rate and degradation rate of IKKs, $b_{10}$, $c_{10}$ and $f_{10}$ represents the activation rate of RIP1, NIK and TAK1 on IKKs, respectively. $a_{11}$ and $d_{11}$ represent the growth rate and degradation rate of MKK4/7, $b_{11}$ and $c_{11}$ represents the activation rate of TAK1 and NF-κB on IKKs. $a_{12}$ and $d_{12}$ represent the growth rate and degradation rate of FADD, $b_{12}$ represents the activation rate of TRADD on FADD. $a_{13}$ and $d_{13}$ represent the growth rate and degradation rate of JNK1/2, $b_{13}$ represents the activation rate of MKK4/7 on JNK1/2. $a_{14}$ and $d_{14}$ represent the growth rate and degradation rate of NF-κB, $b_{14}$ represents the activation rate of IKKs on NF-κB. $a_{15}$ and $d_{15}$ represent the growth rate and degradation rate of ITCH, $b_{15}$ represents the activation rate of JNK1/2 on ITCH. $a_{16}$ and $d_{16}$ represent the growth rate and degradation rate of Caspase10, $b_{16}$ represents the activation rate of FADD on Caspase10. $a_{17}$ and $d_{17}$ represent the growth rate and degradation rate of c-FLIP, $b_{17}$ represents the activation rate of IKKs on c-FLIP and $c_{17}$ represents the inhibiting rate of ITCH on c-FLIP. $a_{18}$ and $d_{18}$ represent the growth rate and degradation rate of Caspase8, $b_{18}$ represents the activation rate of c-FLIP on Caspase8 and $c_{18}$ represents the inhibiting rate of FADD on Caspase8. $a_{19}$ and $d_{19}$ represent the growth rate and degradation rate of Caspase3/7, $b_{19}$ and $c_{19}$ represents the activation rate of Caspase10 and Caspase8 on Caspase3/7. Subsequently, a comprehensive mathematical model for the entire system is established. The time delay of activation or inhibition of downstream signaling molecules varies among different signaling molecules. However, for the sake of convenience in calculation and simulation, a uniform time delay $\tau$ is adopted in this study. In subsequent mathematical modeling, different time delays will be utilized. In this paper, the Eq. 2 was operated by Matlab 2020b.

$\left\{\begin{array}{l}\frac{d^{\alpha }{x}_1\left(t\right)}{dt}=a_1-d_1{x}_1\left(t\right),\\ \frac{d^{\alpha }{x}_2\left(t\right)}{dt}=a_2+b_2{x}_1\left(t-\tau \right)x_2\left(t-\tau \right)-d_2{x}_2\left(t\right),\\ \frac{d^{\alpha }{x}_3\left(t\right)}{dt}=a_3-d_3{x}_3\left(t\right),\\ \frac{d^{\alpha }{x}_4\left(t\right)}{dt}=a_4+b_4{x}_3\left(t-\tau \right)x_4\left(t-\tau \right)-d_4{x}_4\left(t\right),\\ \frac{d^{\alpha }{x}_5\left(t\right)}{dt}=a_5+b_5{x}_4\left(t-\tau \right)x_5\left(t-\tau \right)-d_5{x}_5\left(t\right),\\ \frac{d^{\alpha }{x}_6\left(t\right)}{dt}=a_6+b_6{x}_5\left(t-\tau \right)x_6\left(t-\tau \right)-d_6{x}_6\left(t\right),\\ \frac{d^{\alpha }{x}_7\left(t\right)}{dt}=a_7+b_7{x}_6\left(t-\tau \right)x_7\left(t-\tau \right)-d_7{x}_7\left(t\right),\\ \frac{d^{\alpha }{x}_8\left(t\right)}{dt}=a_8+b_8{x}_6\left(t-\tau \right)x_8\left(t-\tau \right)-d_8{x}_8\left(t\right),\\ \frac{d^{\alpha }{x}_9\left(t\right)}{dt}=a_9+b_9{x}_2\left(t-\tau \right)x_9\left(t-\tau \right)+c_9{x}_6\left(t-\tau \right)x_9\left(t-\tau \right)-d_9{x}_9\left(t\right),\\ \frac{d^{\alpha }{x}_{10}\left(t\right)}{dt}=a_{10}+\left[b_{10}{x}_7\left(t-\tau \right)+c_{10}{x}_8\left(t-\tau \right)+f_{10}{x}_9\left(t-\tau \right)\right]x_{10}\left(t-\tau \right)-d_{10}{x}_{10}\left(t\right),\\ \frac{d^{\alpha }{x}_{11}\left(t\right)}{dt}=a_{11}+b_{11}{x}_9\left(t-\tau \right)x_{11}\left(t-\tau \right)+c_{11}{x}_{14}\left(t-\tau \right)x_{11}\left(t-\tau \right)-d_{11}{x}_{11}\left(t\right),\\ \frac{d^{\alpha }{x}_{12}\left(t\right)}{dt}=a_{12}+b_{12}{x}_5\left(t-\tau \right)x_{12}\left(t-\tau \right)-d_{12}{x}_{12}\left(t\right),\\ \frac{d^{\alpha }{x}_{13}\left(t\right)}{dt}=a_{13}+b_{13}{x}_{11}\left(t-\tau \right)x_{13}\left(t-\tau \right)-d_{13}{x}_{13}\left(t\right),\\ \frac{d^{\alpha }{x}_{14}\left(t\right)}{dt}=a_{14}+b_{14}{x}_{10}\left(t-\tau \right)x_{14}\left(t-\tau \right)-d_{14}{x}_{14}\left(t\right),\\ \frac{d^{\alpha }{x}_{15}\left(t\right)}{dt}=a_{15}+b_{15}{x}_{13}\left(t-\tau \right)x_{15}\left(t-\tau \right)-d_{15}{x}_{15}\left(t\right),\\ \frac{d^{\alpha }{x}_{16}\left(t\right)}{dt}=a_{16}+b_{16}{x}_{12}\left(t-\tau \right)x_{16}\left(t-\tau \right)-d_{16}{x}_{16}\left(t\right),\\ \frac{d^{\alpha }{x}_{17}\left(t\right)}{dt}=a_{17}+b_{17}{x}_{10}\left(t-\tau \right)x_{17}\left(t-\tau \right)-c_{17}{x}_{15}\left(t-\tau \right)x_{17}\left(t-\tau \right)-d_{17}{x}_{17}\left(t\right),\\ \frac{d^{\alpha }{x}_{18}\left(t\right)}{dt}=a_{18}+b_{18}{x}_{17}\left(t-\tau \right)x_{18}\left(t-\tau \right)-c_{18}{x}_{12}\left(t-\tau \right)x_{18}\left(t-\tau \right)-d_{18}{x}_{18}\left(t\right),\\ \frac{d^{\alpha }{x}_{19}\left(t\right)}{dt}=a_{19}+b_{19}{x}_{16}\left(t-\tau \right)x_{19}\left(t-\tau \right)+c_{19}{x}_{18}\left(t-\tau \right)x_{19}\left(t-\tau \right)-d_{19}{x}_{19}\left(t\right).\end{array}\right.(3)$

4 Real-time PCR experiment

4.1 Extract total RNA by Trizol method

(1) Pre-processing of samples

Cell: Add 1 mL Trizol reagent, mix thoroughly by pipetting, transfer to a RNase-free 1.5 mL EP tube, and incubate for 10 min to lyse the cells.

(2) Add 200 μL chloroform, vigorously shake the tube several times, and let it stand at room temperature for 5 min.

(3) Centrifuge at 4°C and 12,000 rpm for 15 min to separate the sample into three phases: the upper phase (RNA), the middle phase (protein), and the lower phase (DNA).

(4) Transfer the upper aqueous phase (approximately 400 μL) to a new 1.5 mL EP tube, add 400 μL isopropanol, mix thoroughly, and let it stand at room temperature for 10 min.

(5) Centrifuge at 4°C and 12,000 rpm for 10 min, and a white RNA pellet should be visible at the bottom of the tube.

(6) Discard the supernatant, add 1 mL 75% ethanol without RNase, vortex, and centrifuge at 4°C and 10,000 rpm for 5 min.

(7) Repeat step 6 once.

(8) Discard the supernatant, air-dry the RNA pellet for 5–10 min, and dissolve the pellet in 20 μL DEPC water.

(9) Measure the OD260, OD280, and OD260/OD280 values of the RNA using a micro spectrophotometer, and calculate the purity and concentration of the RNA. Estimate the RNA quality based on the OD260/OD280 ratio, which should be between 1.8 and 2.0 for experimental requirements. Calculate the total RNA concentration (μg/μL) using the following formula: Total RNA concentration = OD260 × 40×10⁻³.

(10) Store the total RNA at −80°C for future use.

4.2 Reverse transcribed into cDNA

4.2.1 Genomic DNN removal reaction

Melt the RNA template, NASe-free Water, and 10× gDNA Remover Buffer on ice. In a nuclease-free micro centrifuge tube, prepare a reaction system (10 μL) on ice according to the table below. The reaction system is shown in Supplementary Table S1.

Subsequently, the mixed system was blown and briefly centrifuged using a pipettor, and the reaction was carried out on a PCR instrument. The procedure is given as follows:

The cells were incubated at 42°C for 2 min and at 60°C for 5 min.

4.2.2 Reverse transcription reaction

Following the preceding reaction, the system was rapidly chilled on ice and subsequently subjected to a brief centrifugation. The addition of specific components was then carried out based on the desired detection index (Supplementary Tables S2, S3).

Subsequently, employing a liquid motion beat blending system and brief centrifugation, the sample was added to the PCR reaction with the following program: incubation at 25°C for 10 min, followed by incubation at 50°C for 15 min and then further incubated at 85°C for an additional 5 min. The reverse transcription product was subsequently placed on ice or refrigerated as required.

4.3 Real-time PCR detection

The Real-time PCR reaction system is shown in Table 2. Subsequently, the mixed system was blown with a pipette and briefly centrifuged, and then placed on a fluorescent quantitative PCR instrument for amplification detection, the procedure is given in Supplementary Table S5.

TABLE 2

TABLE 2. Real-time PCR reaction system.

The operation data was obtained, and the final data was analyzed by 2^-△△Ct method (as shown in Table 3). A more intuitive bar graph is shown in Supplementary Figure S1.

TABLE 3

TABLE 3. Gene relative expression of IRAK1 (x₂), TNF-α (x₃), TRADD (x₅), RIP1(x₇), NIK (x₈), TAK1 (x₉), IKK-β (x₁₀), MKK4(x₁₁), FADD (x₁₂), ITCH (x₁₅), Caspase-8 (x₁₈) and Caspase-3 (x₁₉) with the interval was 6 h from 0 h to 72 h.

5 Parameter estimate and numerical simulations

According to the expression of various signaling molecules in the TNF signaling pathway, we divided parameter estimation into three stages. The first stage was from 0^th h to 36th h, the second stage was from 36th h to 60th h, and the third stage was from 60th h to 72nd h. The nonlinear least square method was programmed to estimate the 62 parameters by Matlab 2020b and the parameters of stage 1 are as follows (see Table 4):

TABLE 4

TABLE 4. The estimated parameters of stage 1 for numerical simulation.

The parameters of stage 2 are as follows (see Table 5):

TABLE 5

TABLE 5. The estimated parameters of stage 2 for numerical simulation.

The parameters of stage 3 are as follows (see Table 6):

TABLE 6

TABLE 6. The estimated parameters of stage 3 for numerical simulation.

Based on the model parameters obtained from three stages and combined with mathematical model (2), we obtained simulation results of 19 genes, as shown in Figures 3–7. The red dots with error bars represent the real experimental data, while the deep blue curves are numerical simulations of the model, and the light blue areas are three times variance of fitted curves.

FIGURE 3

FIGURE 3. The gene relative expression of IL-1β, IRAK1, TNF-α and TNFR1 from 0 h to 72 h.

FIGURE 4

FIGURE 4. The gene relative expression of TRADD, TRAF2/5, RIP1 and NIK from 0 h to 72 h.

FIGURE 5

FIGURE 5. The gene relative expression of TAK1, IKK-β, MKK4 and FADD from 0 h to 72 h.

FIGURE 6

FIGURE 6. The gene relative expression of JNK1/2, NF-kB, ITCH and Caspase-10 from 0 h to 72 h.

FIGURE 7

FIGURE 7. The gene relative expression of c-FLIP, Caspase-8 and Caspase-3 from 0 h to 72 h.

The correlation coefficient (R) was calculated to compare the experimental findings with model stimulation. The results, as presented in Table 7, demonstrate that all tested variables exhibited relatively strong correlations, indicating the efficiency of the model.

TABLE 7

TABLE 7. The correlation coefficient between experimental data and numerical simulation.

6 The sensitivity analysis of molecules based on the mathematical model

Based on the accurate mathematical model obtained, the influence of model parameter changes on the system will be studied next, that is, the sensitivity of parameters in the model will be analyzed. Given the positive correlation between chondrocyte apoptosis and osteoarthritis severity, this project aims to employ gene interference experiments to silence related genes while making subtle adjustments to model parameters that significantly impact output. The ultimate goal is to reduce chondrocyte apoptosis and identify new targets and strategies for treating osteoarthritis. In this section, sensitivity of molecules based on the mathematical model was analyzed by numerical simulations. Because the severity of chondrocyte apoptosis is positively correlated with the severity of osteoarthritis, this project aims to greatly change the output of the model by making subtle adjustments to the sensitive parameters of the model, that is, to reduce chondrocyte apoptosis, and to provide new targets and strategies for the treatment of osteoarthritis. Therefore, in this study, the degree of change in the mRNA expression of the gene Caspase-3 was used as an indicator of parameter sensitivity. The simulation results show that $b_4$, $a_5$, $d_5$, $a_6$, $a_{10}$, $b_{10}$, $g_{10}$, $f_{10}$, $d_{10}$, $b_{12}$, $d_{12}$, $a_{16}$, $b_{16}$, $d_{16}$, $a_{17}$, $a_{18}$, $b_{18}$, $a_{19}$, $b_{19}$, $g_{19}$, $d_{19}$ are sensitive parameters and others are insensitive parameters. Below we randomly selected a sensitive parameter ($a_5$) and demonstrated its simulation results. The results are as follows:

Similarly, an insensitive parameter ($b_7$) was randomly selected and the simulation results were demonstrated. The results are shown in Supplementary Figures S2–S6.

In this section, all parameters in the model were numerically simulated (see Supplementary Material). We randomly selected parameter $a_5$ to demonstrate the sensitivity of parameters on Caspase-3. The numerical simulation results show that when parameter $a_5$ increases from 0.6 to 3.8 with a step size of 0.8, the signaling molecules IL-1β, IRAK1, TNF-α, and TNFR1 show no significant changes (as shown in Figure 8). This is because parameter $a_5$ represents the production rate of the signaling molecule TRADD, which does not have a feedback regulatory effect on upstream signaling molecules. From Figure 9, it can be clearly seen that as parameter $a_5$ increases from 0.6 to 3.8, the mRNA expression level of TRADD gradually increases, and the growth amplitude is large. The black, green, red, light blue, and dark blue curves represent the dynamic trends of TRADD over time when $a_5$ is 0.8, 1.4, 2.0, 2.6, 3.2, and 3.8, respectively. Similarly, with the increase of $a_5$, the mRNA expression levels of TRAF2/5 and NIK also gradually increase, and the increase is significant. Especially when $a_5$ is 3.8, the mRNA expression level of NIK is approximately 10 times that of $a_5=0.8$. As for RIP1, although the expression level increases with the increase of $a_5$, the increase is relatively small. From Figure 10, it can be seen that when $a_5$ increases from 0.8 to 3.2, the expression levels of TAK1, IKK-β, and MKK4 only increase slightly, but when $a_5$ increases to 3.8, the expression levels of TAK1, IKK-β, and MKK4 increase significantly. The expression level of FADD also increases gradually with the gradient increase of $a_5$. Despite the large increase in $a_5$ from 0.8 to 3.8, JNK1/2, NF-kB, and ITCH show no significant changes. However, with the increase of $a_5$, the mRNA expression level of Caspase-10 gradually increases, and when $a_5$ is 3.8, the mRNA expression level of Caspase-10 is significantly higher than the other cases (as shown in Figure 11). From Figure 12, it can be seen that when $a_5$ increases from 0.8 to 3.2, the mRNA expression levels of c-FLIP and Caspase-8 increase but not significantly, and the mRNA expression level of Caspase-3 also increases with the gradient increase of $a_5$. However, when $a_5$ increases to 3.8, c-FLIP, Caspase-8, and Caspase-3 all undergo a stepwise increase. For example, when $a_5$ is 3.8, c-FLIP increases by 1.38 times, Caspase-8 increases by 1.8 times, and Caspase-3 increases by 1.27 times. Therefore, from the numerical simulation results, it can be seen that small changes in parameter $a_5$ can significantly affect the mRNA expression of downstream signaling molecules in the signaling pathway, especially the indicator signaling molecule Caspase-3, which is indicative of the severity of osteoarthritis. The mRNA expression level also undergoes significant changes. Therefore, parameter $a_5$ is a sensitive parameter that affects the progression of osteoarthritis and is a potential therapeutic target for treating osteoarthritis.

FIGURE 8

FIGURE 8. The dynamic trend of IL-1β, IRAK1, TNF-α and TNFR1 from 0 h to 72 h when the parameter $a_5$ increases from 0.6 to 3.8, and the interval is 0.8.

FIGURE 9

FIGURE 9. The dynamic trend of TRADD, TRAF2/5, RIP1 and NIK from 0 h to 72 h when the parameter $a_5$ increases from 0.6 to 3.8, and the interval is 0.8.

FIGURE 10

FIGURE 10. The dynamic trend of TAK1, IKK-β, MKK4 and FADD from 0 h to 72 h when the parameter $a_5$ increases from 0.6 to 3.8, and the interval is 0.8.

FIGURE 11

FIGURE 11. The dynamic trend of JNK1/2, NF-kB, ITCH and Caspase-10 from 0 h to 72 h when the parameter $a_5$ increases from 0.6 to 3.8, and the interval is 0.8.

FIGURE 12

FIGURE 12. The dynamic trend of c-FLIP, Caspase-8 and Caspase-3 from 0 h to 72 h when the parameter $a_5$ increases from 0.6 to 3.8, and the interval is 0.8.

On the other hand, when parameter $b_7$ increases from 0.04 to 0.12 with a step size of 0.02, only the mRNA expression level of the signal molecule RIP1 in the entire signal pathway shows a significant change, while the other signal molecules show no significant changes. The numerical simulation results are shown in Supplementary Figures S2–S6. Therefore, parameter $b_7$ has no significant impact on osteoarthritis and is considered an insensitive parameter.

7 Conclusion

Since most of the current research on the regulation of signal pathways in chondrocyte apoptosis is based on molecular biology methods, it can only focus on studying the signaling molecules in a specific pathway, making it difficult to understand the mechanism of signal transduction at a systemic level. In this study, high-throughput sequencing is firstly used to screen for differentially expressed genes involved in IL-1β-induced apoptosis in human chondrocytes. The results show significant enrichment in the TNF signaling pathway. Therefore, mathematical modeling approach is adopted to quantitatively analyze the mechanism of TNF signaling pathway in regulating chondrocyte apoptosis. Next, qPCR experiments are performed to measure the mRNA expression levels of IRAK1 (x₂), TNF-α (x₃), TRADD (x₅), RIP1(x₇), NIK (x₈), TAK1 (x₉), IKK-β (x₁₀), MKK4(x₁₁), FADD (x₁₂), ITCH (x₁₅), Caspase-8 (x₁₈) and Caspase-3 (x₁₉) at 0 h, 6 h, 12 h, 18 h, 24 h, 30 h, 36 h, 42 h, 48 h, 54 h, 60 h, 66 h, and 72 h. The 62 parameters of the model are estimated by using the collected experimental data. Since the process of chondrocyte apoptosis is time-dependent, the parameter estimation was divided into three stages. The effectiveness of parameter estimation is evaluated by calculating the correlation coefficient between the experimental data and the mathematical model, with a maximum correlation coefficient of 0.969 (See Table 7). Finally, numerical simulation are used to calculate the sensitivity of the model parameters, and the results show that parameters $b_4$, $a_5$, $d_5$, $a_6$, $a_{10}$, $b_{10}$, $g_{10}$, $f_{10}$, $d_{10}$, $b_{12}$, $d_{12}$, $a_{16}$, $b_{16}$, $d_{16}$, $a_{17}$, $a_{18}$, $b_{18}$, $a_{19}$, $b_{19}$, $g_{19}$ and $d_{19}$ are sensitive parameters. These sensitive parameters significantly affect chondrocyte apoptosis and further influence the severity of osteoarthritis. Therefore, in future studies, we will design siRNA interference experiments to validate the sensitive and insensitive parameters at the cellular level, further explore the sensitive factors that affect chondrocyte apoptosis, and potentially identify new therapeutic targets for the treatment of osteoarthritis at a systemic level.

Data availability statement

The platform for the high-throughput sequencing is Health Time Gene and the sequencing data were uploaded to NCBI, the number is PRJNA1019114.

Ethics statement

The studies involving humans were approved by the first affiliated hospital of shenzhen university (shenzhen second people’s hospital) clinical research ethics committee. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Author contributions

YW: Data curation, Methodology, Project administration, Writing–original draft. JL: Conceptualization, Investigation, Software, Writing–original draft. BH: Formal Analysis, Funding acquisition, Project administration, Resources, Writing–original draft. XL: Supervision, Validation, Visualization, Writing–review and editing. XS: Software, Supervision, Writing–review and editing. DS: Conceptualization, Funding acquisition, Supervision, Writing–review and editing, Software, Validation.

Funding

The authors declare financial support was received for the research, authorship, and/or publication of this article. The work is supported by the National Natural Science Foundation of China (62103287), Guangdong Provincial Natural Science Foundation General Project (2023A1515010518), Basic Research General Project of Shenzhen (JCYJ20210324103209026), PhD Basic Research Initiation Project (RCBS20200714114856171, RCBS20200714114909099), Shenzhen Nanshan District Excellent Youth Fund Project (NSZD2023061), President's Fund of Southern University of Science and Technology Hospital (2022-A3) and Clinical Research Project of Shenzhen Second People’s Hospital (20213357007).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2023.1288431/full#supplementary-material

References

Alahdal, M., Sun, D. S., Duan, L., Ouyang, H., Wang, M., Xiong, J., et al. (2021). Forecasting sensitive targets of the kynurenine pathway in pancreatic adenocarcinoma using mathematical modeling. Cancer Sci. 00, 1481–1494. doi:10.1111/cas.14832

CrossRef Full Text | Google Scholar

Cai, W., and Li, G. (2019). Research progress of signaling pathway in the pathogenesis of osteoarthritis. Guangxi Med. J. 41 (06), 93–96.

Google Scholar

Frank, M., Bwemero, J., Kalunga, D., Sangu, W., Semeni, S., Hamisi, M., et al. (2015). OA60 Public health and palliative care mix; a ccpmedicine approach to reverse the overgrowing burden of non-communicable diseases in Tanzania. Support. Palliat. Care 5 (1), A19. doi:10.1136/bmjspcare-2015-000906.60

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, S., Yin, H., Liu, H., and Sui, Y. H. (2014). Research progress on MAPK signal pathway in the pathogenesis of osteoarthritis. China J. Orthop. Trauma 27 (5), 441–444.

PubMed Abstract | Google Scholar

Gu, Y., Liu, S., Xia, S., Liao, Q., and Wang, W. (2016). Research progress in pathogenesis and treatment of osteoarthritis. Chin. J. Bone Jt. 5 (10), 770–774.

Google Scholar

Hu, S., Mao, G., Zhang, Z., Wu, P., Wen, X., Liao, W., et al. (2019). MicroRNA-320c inhibits development of osteoarthritis through downregulation of canonical Wnt signaling pathway. Life Sci. 228, 242–250. doi:10.1016/j.lfs.2019.05.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, S. L., Wang, K., Shi, Y. F., Shao, Z. X., Zhang, C. X., Sheng, K. W., et al. (2020). Downregulating Akt/NF-κB signaling and its antioxidant activity with Loureirin A for alleviating the progression of osteoarthritis: in vitro and vivo studies. Int. Immunopharmacol. 78, 105953. doi:10.1016/j.intimp.2019.105953

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, B., Guo, W., Ma, H., Xu, B., Mu, W., Zhang, Z., et al. (2017). Isoliquiritigenin suppresses IL-1β induced apoptosis and inflammation in chondrocyte-like ATDC5 cells by inhibiting NF-κB and exerts chondroprotective effects on a mouse model of anterior cruciate ligament transection. Int. J. Mol. Med. 40 (6), 1709–1718. doi:10.3892/ijmm.2017.3177

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalamakis, G., Brüne, D., Ravichandran, S., Bolz, J., Fan, W., Ziebell, F., et al. (2019). Quiescence modulates stem cell maintenance and regenerative capacity in the aging brain. Cell. 176 (6), 1407–1419. doi:10.1016/j.cell.2019.01.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Hui, C., Ping, Z., Zhou, S. H., Tian, Q., Shi, J., et al. (2016). JAK2/STAT3 signal pathway mediating curcumin in cartilage cell metabolism of osteoarthritis. China J. Orthop. Trauma 12, 1104–1109. doi:10.3969/j.issn.1003-0034.2016.12.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, F., Li, L., Lu, W., Ding, Z., Huang, W., Li, Y. T., et al. (2020). Scutellarin ameliorates cartilage degeneration in osteoarthritis by inhibiting the Wnt/β-catenin and MAPK signaling pathways. Int. Immunopharmacol. 78, 105954. doi:10.1016/j.intimp.2019.105954

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, F., Sun, D. S., Murakami, R., and Matsuno, H. (2017). Modeling and analysis of the Delta-Notch dependent boundary formation in the Drosophila large intestine. BMC Syst. Biol. 11 (4), 80–60. doi:10.1186/s12918-017-0455-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., He, X., Ping, Z., Zhou, S., and Li, X. (2016). Protective effect of diosgenin on chondrocytes mediated by JAK2/STAT3 signaling pathway in mice with osteoarthritis. J. Zhejiang Univ. Med. Sci. 45 (5), 452–459.

PubMed Abstract | Google Scholar

Ma, G., and Ren, M. (2012). Reserch on apoptosis of chondrocytes in osteoarthritis. Chin. J. Histochem. Cytochem. 21 (2), 196–200.

Google Scholar

Manu, S. S., Spirov, A. V., Gursky, V. V., Janssens, H., Kim, A. R., Ovidiu, R., et al. (2009). Canalization of gene expression and domain shifts in the Drosophila blastoderm by dynamical attractors. PLoS Comput. Biol. 5 (3), e1000303. doi:10.1371/journal.pcbi.1000303

PubMed Abstract | CrossRef Full Text | Google Scholar

Matsuno, H., Murakami, R., Yamane, R., Yamasaki, N., Fujita, S., Yoshimori, H., et al. (2003). Boundary formation by notch signaling in Drosophila multicellular systems: experimental observations and gene network modeling by Genomic Object Net. Pac Symp. Biocomput 8, 152–163.

PubMed Abstract | Google Scholar

Nakamasu, A., Takahashi, G., Kanbe, A., and Kondo, S. (2009). Interactions between zebrafish pigment cells responsible for the generation of Turing patterns. Proc. Natl. Acad. Sci. U. S. A. 106 (21), 8429–8434. doi:10.1073/pnas.0808622106

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, C., Gan, H., and Huang, D. (2007). Molecular mechanism of TNF-α signal transduction. Chin. J. Biochem. Mol. Biol. 23 (6), 430–435.

Google Scholar

Rigoglou, S., and Papavassiliou, A. G. (2013). The NF-κB signalling pathway in osteoarthritis. Int. J. Biochem. Cell. Biol. 45 (11), 2580–2584. doi:10.1016/j.biocel.2013.08.018

CrossRef Full Text | Google Scholar

Sun, D., Duan, Li, Wang, D., and Xiong, J. (2020). Mathematical modeling and infection control strategy of COVID-19. Chin. J. Dis. Control Prev. 24 (5), 523–528.

Google Scholar

Sun, D. S., and Liu, F. (2017). Analysis of a new delayed HBV model with exposed state and immune response to infected cells and viruses. Biomed Res. Int. 2017, 7805675. doi:10.1155/2017/7805675

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, D. S., and Liu, F. (2018). Modeling and control of a delayed hepatitis B virus model with incubation period and combination treatment. Interdiscip. Sci. Comput. Life Sci. 10, 375–389. doi:10.1007/s12539-017-0275-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Toshihisa, K. (2016). Cell death in chondrocytes, osteoblasts, and osteocytes. Int. J. Mol. Sci. 17 (12), 2045–2071.

PubMed Abstract | Google Scholar

Wang, J., Yu, T., Zheng, Z., Lv, H., Chen, W., and Zhang, Y. (2019). Research progress of pathological mechanism of knee osteoarthritis. J. hebei Med. Univ. 40 (10), 1237–1238.

Google Scholar

Wang, W., Fu, L., and Ye, J. (2005). Research progress of TNF-α signaling pathway. J. Fujian Med. Univ. 39, 27–31.

Google Scholar

Wu, S., Liu, J., Zuo, Y., and Li, Z. (2018). The effect of Notch signaling pathway on apoptosis of articular chondrocytes in knee osteoarthritis. West China Med. J. 33 (09), 115–120.

Google Scholar

Yoon, H. S., and Kim, H. A. (2004). Prologation of c-jun N-terminal kinase is associated with cell death induced by tumor necrosis factor alpha in human chondrocytes. J. Korean Med. Sci. 19 (4), 567–573. doi:10.3346/jkms.2004.19.4.567

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, P., Kang, W., Li, X., Dong, B., Liu, D., and Wang, W. (2016). Research progress on the pathogenesis of osteoarthritis and related cytokines. Orthop. J. China 11, 1010–1015.

Google Scholar

Zahoor, T., Mitchell, R., Bhasin, P., Guo, Y., Paudel, S., Schon, L., et al. (2017). Effect of low-intensity pulsed ultrasound on joint injury and post-traumatic osteoarthritis: an animal study. Ultrasound Med. Biol. 44 (1), 234–242. doi:10.1016/j.ultrasmedbio.2017.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, H., and Gong, N. (2019). miR-20a regulates inflammatory in osteoarthritis by targeting the IκBβ and regulates NK-κB signaling pathway activation. Biochem. Biophysical Res. Commun. 518 (4), 632–637. doi:10.1016/j.bbrc.2019.08.109

PubMed Abstract | CrossRef Full Text | Google Scholar

Ziebell, F., Dehler, S., Martin-Villalba, A., and Marciniak-Czochra, A. (2018). Revealing age-related changes of adult hippocampal neurogenesis using mathematical models. Development 145, dev153544–12. doi:10.1242/dev.153544

PubMed Abstract | CrossRef Full Text | Google Scholar