The osteogenetic activities of mesenchymal stem cells in response to Mg2+ ions and inflammatory cytokines: a numerical approach using fuzzy logic controllers

Magnesium (Mg2+) ions are frequently reported to regulate osteogenic activities of mesenchymal stem cells (MSCs). In this study, we propose a numerical model to study the regulatory importance of Mg2+ ions on MSCs osteoblastic differentiation in the presence of an inflammatory response. A fuzzy logic controller was formulated to receive the concentrations of Mg2+ ions and the inflammatory cytokines of TNF-α, IL-10, IL-1β, and IL-8 as cellular inputs and predict the cells’ early and late differentiation rates. Five sets of empirical data obtained from published cell culture experiments were used to calibrate the model. The model successfully reproduced the empirical data regarding the concentration- and phase-dependent effect of Mg2+ ions on the differentiation process. In agreement with the experiments, the model showed the stimulatory role of Mg2+ ions on the early differentiation phase, once administered at low concentration, and their inhibitory role on the late differentiation phase. The numerical approach used in this study suggested 6–8 mM as the most effective concentration of Mg2+ ions in promoting the early differentiation process. Also, the proposed model sheds light on the fundamental differences in the behavioral properties of cells cultured in different experiments, e.g. differentiation rate and the sensitivity of the cultured cells to stimulatory signals such as Mg2+ ions. Thus, it can be used to interpret and compare different empirical findings. Moreover, the model successfully reproduced the nonlinearities in the concentration-dependent role of the inflammatory cytokines in early and late differentiation rates. Overall, the proposed model can be employed in studying the osteogenic properties of Mg-based implants in the presence of an inflammatory response.


Introduction
Magnesium (Mg)-based biomaterials are an attractive choice in the orthopedic industry due to their biodegradability and superior osteogenic capacity compared to non-degradable metallic implants e.g., titanium alloys [1][2][3]. Mg 2+ ions released in the implantation site due to the degradation process modulate a wide range of physiological processes involved in bone fracture healing, in particular osteogenesis [1], [4]. The mediatory effects of Mg 2+ ions on osteogenic differentiation can be evaluated in two aspects (see Fig 1A). First, Mg 2+ ions directly modulate multiple signaling pathways associated with the differentiation of mesenchymal stem cells (MSCs) to osteoblasts [5], [6]; Mg 2+ ions are shown to activate MAPK/ERK and Wnt/ β-catenin signaling pathways which are associated with osteogenic differentiation [7]. Secondly, Mg 2+ ions mediate the inflammatory response and thereby indirectly guide osteogenesis [8][9][10].
The inflammatory response is known to play a pivotal role in all stages of bone tissue regeneration [11]- [13]. In particular, macrophages are repeatedly shown to guide osteogenesis [12], [14], [11]. During the healing process, macrophages obtain different functional roles ranging from pro-inflammatory M1 to anti-inflammatory M2 and secrete a multitude of cytokines essential for the osteogenic ability of MSCs [14], [15]. Empirical data have shown the significance of Mg 2+ ions in modulating macrophage polarization and cytokines production [9], [10], [16], [17]. For instance, Qiao et al. [17] applied Mg 2+ ions to macrophages in a cell culture environment and observed an upregulation in the production of interleukin (IL)-8, which is an important factor in the osteogenic differentiation of MSCs. Such an indirect influence of Mg 2+ ions on osteogenesis is reported to even overweight the direct regulatory effect of Mg 2+ ions on MSCs [17].
Although Mg 2+ ions are generally shown to promote osteogenesis, multiple studies have reported their detrimental effects on tissue regeneration [8], [18]. Such a conflicting finding implies an incomplete understanding of the role of Mg 2+ ions in the complex process of fracture healing. In order to effectively design a Mg-based implant, it is paramount to fully understand the bioregulatory role of Mg 2+ ions on tissue regeneration in the presence of other signaling factors, in particular inflammatory reactions [4], [19] [20]. So far, the experimental approach has been the primary method in investigating the physiological roles of Mg 2+ ions. In recent years, numerous studies have shown the importance of numerical investigation in the study of biological systems and bone-implant design [21][22][23][24][25]. Fuzzy logic (FL) has been given a special attention in quantitative modelling of biological systems with uncertain kinetic data [26] [27]. Aldridge et al [28] used FL in describing the dynamics of intracellular pathways of human colon carcinoma cells associated with different growth factors and insulin receptors. Their FL-based simulations successfully produced several predictions of pathways crosstalk and regulation [28]. They also proposed a relationship between MK2 and ERK pathways which was unknown previously [28]. Wang et al [29] and Niemeyer et al [30] employed FL to simulate the transitional formation of different tissues during bone regeneration. They have successfully shown the patterns of bone tissue regeneration under various mechanical environment [29], [30].
Recently, we proposed a fuzzy agent-based model to numerically study the proliferation and osteogenic differentiation of MSCs in response to Mg 2+ ions [31]. We showed that Mg 2+ ions within 3-6mM concentration have the highest stimulation effect on cell population growth. The model also captured the stimulatory role of Mg 2+ ions on early differentiation and its inhibitory effect on the late differentiation process. In this study, as an extension of our previous model, we propose a FL-based model to investigate the osteogenic response of MSCs to Mg 2+ ions in tandem with the inflammatory cytokines of tumor necrosis factor alpha (TNFα), interleukin 10 (IL-10), interleukin 1 beta (IL-1β), and IL-8. In the present model, the previously proposed controller for the Mg 2+ ions is further extended to include the inhibitory role of Mg 2+ ions on the early osteogenic response once the concentrations is below the physiological level. In addition, five more controllers are presented to account for the regulatory roles of the given inflammatory cytokines.

The overview of the model formulation
In this study, the proposed cell model receives five cellular inputs and predicts the osteogenic differentiation rates (see Fig 1B). The choice of cytokines is based on their importance in osteogenesis and their involvement in the interplay of Mg 2+ ions with the inflammatory cells, as elaborated in section 4.1. Since the cellular inputs can regulate the early and late differentiation process differently (Fig 1C), two distinguished processes of early-and late differentiation rates are calculated by the model. A FL controller is designed to act as the core calculator of the cell model. In order to define fuzzy rules, we gathered the available information in the literature regarding the regulatory role of the cellular inputs (see section 4.1). A complete introduction to the development of the FL controller is given in section 5. The formulation of the model generated 30 unknown parameters given in S3 Table. We estimated the values of these parameters using a calibration process by employing differential evolution (DE) [32]. The published data of five cell-culture experiments were used for this purpose (see section 5.5). These experiments reported the bioregulatory effect of Mg2+ ions and different inflammatory cytokines on the osteoblastic differentiation of MSC. In these reports, alkaline phosphate (ALP) was used as the early differentiation marker, while osteocalcin (OC) and alkaline red staining (ARS) were the late differentiation markers. The parameter estimation process was carried out to maximize the fitness value, which is defined as the normalized absolute difference between the simulation results (S) and the empirical observations (E) as | R 2 = |E-S|/E. The calibration process was conducted on the dataset of each experiment individually, encoded as C1 to C5, as well as on the combined data of all experiments, encoded as C1-5. Once the model was calibrated, we investigated the sensitivity of the model to the parameters using large-scale simultaneous perturbations (LSSP) and small-scale individual perturbations (SSIP). The details of the calibration process and the sensitivity analysis are given in section 5.6.

Simulation results versus empirical observations
2.2.1 Study 1. Study 1 reports the regulatory effect of Mg 2+ ions on the early differentiation rate, measured by ALP [17]. The fits of the simulation results to the data of study 1 are presented in Fig 2. The empirical data suggest that the application of Mg 2+ ions at the concentration of 8 mM increases the early differentiation rate at both measurement days of 3 and 7 compared to the physiological concentration of 0.8 mM (see Fig 2). Contrarily, the application of Mg 2+ ions at the concentration of 0.08 mM decreases the early differentiation rate (see Fig  2). These observations were correctly captured by the model calibrated by C1-5 with a fitness value of 0.84. Once calibrated by C1, the model was able to closely reproduce the reported measurements with the average fitness value of 0.92. The variations in the simulation results of study 1 due to SSIP are given in Fig 2 for [6]. The fits of the model to the data of study 2 are given in Fig 3. The data suggests that by increasing Mg 2+ ions from 0.8 to 5 mM, the early differentiation rate increases while the late differentiation rate decreases. The model calibrated by C1-5 was able to capture the former but not the latter. The mean fitness values obtained during Effects of magnesium degradation products on mesenchymal stem cell fate and osteoblastogenesis C1-5 were 0.94 for ALP and 0.49 for OC (see Fig 3). Once calibrated by C2, the model was capable of correctly reproducing both trends and exact values with the fitness values of 1 for both markers (see Fig 3). The variations in the simulation results of study 2 due to SSIP are given in Fig 3 for different measurement items.

Study 3.
Study 3 reports the regulatory influence of IL-10 and TNF-α on the early and late differentiation rates, measured by ALP and ARS at day 14 and 21, respectively, for the application period of 48 hours [13]. The fits of the model to the data of study 3 are given in Fig  4 for IL-10 and in Fig 5 for TNF-α. In consistent with the empirical data, the model calibrated by C1-5 showed the stimulatory effect of IL-10 on both early and late differentiation rates within the concentration range of 0 to 10 ng/ml (see Fig 4). The obtained fitness values were 0.88 and 0.97 for ALP and ARS, respectively (see Fig 4). Once calibrated by C3, the fits of the model to the data experienced only a slight improvement (see Fig 4). The variations in the simulation results of study 3 due to SSIP are given in Fig 4 for the case of IL-10.
In agreement with the experiments, the model calibrated by C1-5 successfully captured the stimulatory effect of TNF-α on the differentiation rates at the applied concentration of 1ng/ml (see Fig 5). C1-5 produced the fitness values of 0.91 and 0.93 for ALP and ARS, respectively, for the case of TNF-α (see Fig 5). Once calibrated by C3, the fits of the model improved to 0.95 and 0.96 for ALP and ARS, respectively (see Fig 5). The variations in the simulation results of study 3 due to SSIP are given in Fig 5 for the case of TNF-α.

Study 4.
Study 4 reports the effect of IL-10, administrated over 48 hours, on the early and late differentiation rates, measured by ALP and ARS, respectively [33]. The fits of the model to the data of study 4 are given in Fig 6. The data shows the stimulatory effect of IL-10 on the early and later differentiation rates once applied in the concentration range of 0 to 10 ng/ml with a maximum effect at the concentration of 0.1 ng/ml (see Fig 6). However, IL-10   shows an inhibitory effect on the differentiation process once applied at 10 ng/ml or more (see Fig 6). The model calibrated by C1-3 successfully captured these observations with the fitness values of 0.85 and 0.82 for ALP and ARS, respectively (see Fig 6). However, the simulation results notably deviate from the empirical data for the applied concentration of 1 ng/ml (see Fig 6). Once calibrated by C4, the simulation results were improved with the fitness values of 0.95 and 0.96 for ALP and ARS, respectively (see Fig 6). The variations in the simulation results of study 4 due to SSIP are given in Fig 6 for different measurement items.

Study 5.
Study 5 measures the regulatory effects of IL-8 and IL-1β on the early differentiation rate, measured by ALP [17]. The fits of the model to the data of study 5 are presented in Fig 7. Based on the data, IL-8 stimulates the early differentiation rate in a dose-dependent fashion increasing from 0 to 100 ng/ml (see Fig 7). The model calibrated by C1-5 was not capable of reproducing the increase in ALP from 10 to 100 ng/ml (see Fig 7). Once calibrated by C5, the model reproduced the data throughout all applied concentrations. The fitness value was significantly improved from 0.66 to 0.98 from the case of C1-5 to C5 (see Fig 7). In the case of IL-1β, the empirical data show a significant stimulatory effect at the applied concentration of 10 ng/ml. The simulation results reproduced this observation with a fitness value of 0.87 (see Fig 7). Once calibrated by C5, the match of the simulation results to the data improved resulting in a fitness value of 1 (see Fig 7). By simultaneous application of IL-8 and IL-1β, the data showed a significant upregulation in the measured ALP (see Fig 7). The model calibrated by C1-5 reproduced this observation with a fitness value of 0.80 (see Fig 7). Once calibrated by C5, the simulation results perfectly matched the data with the fitness value of 1 (see Fig 7). The variations in the simulation results of study 5 due to SSIP are given in Fig 7 for different measurement items.

Results of the calibration process
The inferred parameter values obtained during different runs of C1-5 are given in Fig 8A. The results of the rest of the calibration scenarios, i.e. C1, C2, C3, C4, and C5 can be found in S1

Results of the sensitivity analysis
The significance of the parameters obtained during different calibration scenarios is provided in Fig 3. It can be seen that LSSP and SSIP suggest different significance orders for the parameters. For instance, in the case of C5, α es (sensitivity of the early differentiation rate to the stimulatory signals) is recognized as the top important parameter using SSIP, while this parameter is not among the top five parameters suggested by LSSP. Overall, the parameters of n ALP and β ALP , which map the simulated maturity to the measured quantity of ALP in experiments, are among the first two significant parameters across different studies. In the second rank are the parameters of M t (early maturity threshold), T d (differentiation time), p ef (the Fast membership level of early differentiation), and α es (sensitivity of the early differentiation rate to the stimulatory signals).

The overall goodness of the model
In summary, the present model calibrated by the accumulated data of all studies, i.e. C1-5, was able to successfully reproduce most of the important empirical observations reported in studies 1 to 5 (Figs 2-7). In agreement with the experiments, the results of the simulations correctly demonstrated that Mg 2+ ions in low concentration stimulates the early differentiation process, while at a concentration below the physiological level downregulates the differentiation phase (see Figs 2 and 3). The calibration process suggested 6-8 mM as the optimal concentration of Mg 2+ ions in promoting the early differentiation rate (see the estimated value of p ms in S3 Table). In addition, the model was capable of reproducing the following observations: • The stimulatory effect of IL-10, once administrated for 48 hours, on the early and late differentiation processes within the concentration range of 0 to 10 ng/ml (see Fig 4) • The stimulatory role of TNF-α once applied at the low concentration of 1 ng/ml (see Fig 5) • The dose-dependent role of IL-10, once administrated over 48 hours, on early and late differentiation processes-an increasing trend in the stimulatory role within the concentration of 0 to 0.1 ng/ml and an inhibitory role at the concentration of above 1 ng/ml (see Fig 9) • The increasing trend in the stimulatory effect of IL-8 on the early differentiation process within the concentration range of 0 to 10 ng/ml (see Fig 9)  • The stimulatory role of IL-1β on the early differentiation rate once applied at the concentration of 10 ng/ml (see Fig 9) • The cumulative stimulatory effect of IL-8 and IL-1β on the early differentiation process once administrated together (see Fig 9).
However, the simulation results deviated from the experiments in capturing the inhibitory effect of Mg 2+ ions on the late differentiation process (see Fig 3). In addition, the model was unable to reproduce the increasing trend in the stimulatory effect of IL-8 on the early differentiation rate once the applied concentration increased from 10 to 100 ng/ml (see Fig 7). Moreover, the simulation results were noticeably different from the empirical data for several observations; for instance, the quantity of ALP in study 3 for the applied concentration of 0.1 ng/ml (see Fig 4) and the quantity of ALP in study 4 for the applied concentration of 1 ng/ml (see Fig 6). In the next step, we conducted the calibration process for each study individually, i.e. C1, C2, C3, C4, and C5. The overall fits of the model to the empirical data improved significantly with an average accuracy of over 95 percent (see Figs 2-7).

Investigating the sources of discrepancy between the data and the simulation
In order to investigate the possible sources of discrepancy among the empirical data of different studies, we examined the estimated parameters obtained during different calibration scenarios (see Fig 8B). The results showed that the values of several parameters were noticeably different between different calibration schemes (see section 2.3). Such discrepancy can stem from either the different exploration of the calibration algorithm in finding the global minimums or the fundamental differences between different studies [31].
As previously noted, the model calibrated by C1-5 could not explain the inhibitory effect of Mg 2+ ions on the late differentiation process, measured by OC in study 2 (see Fig 3). According to Fig 8B, the parameters of M t , T d , β OC , n OC , α li , and α es show divergence between C1-5 and C2 (see Fig 2). β OC and n OC , which map the simulated maturity to OC (see Eq (6)), are not associated with the calculation of the simulated maturity (see Eq (1) in section 4.4). Therefore, we exclude them from further evaluation. To examine the remaining parameters, we first elaborate on the regulatory effect of Mg 2+ ions on the differentiation process. In the present model, the differentiation process is divided into the early and late phases (see Eq (1)). According to the fuzzy rules, Mg 2+ ions stimulate the early differentiation rate while inhibiting the late differentiation rate (see Table 1). Therefore, the outcome of the differentiation progress is partly upregulated and partly downregulated by Mg 2+ ions. M t and T d control the degree of contribution of each part in the final outcome of the differentiation process (see Eqs (1) and (2)). The lower quantity of M t , i.e. maturation threshold, obtained during C2 compared to C1-5 (see Fig  2B) shortens the early maturation phase and consequently constrains the stimulatory role of Mg 2+ ions on the overall differentiation process. The lower quantity of T d , i.e. differentiation time, obtained during C2 compared to C1-5 (see Fig 2B) serves a similar purpose; a lower T d results in an earlier maturation along the differentiation line, consequently limiting the stimulatory influence of Mg 2+ ions. The evaluation of M t and T d suggests that the cells used in the experiments of study 2 possess a lesser differentiation capacity compared to the average of the rest of the studies.
Additionally, α es and α li contribute to the deviation of the model calibrated by C1-5 from the empirical data. α es and α li control the regulatory impact of Mg 2+ ions on early and late differentiation rates, respectively (see Eq (3)). C2 reports a lower quantity of α es and a higher quantity of α li , respectively, compared to C1-5 (see Fig 2B). These parameters, while possessing lesser significance on the simulation results compared to M t and T d (see Fig 9), contribute to reducing the overall effect of Mg 2+ ions on the late differentiation process. This observation implies that cells used in the experiments of study 2 are more sensitive to stimulatory signals during early differentiation phase and less sensitive during late differentiation phase.
We previously noted that the model calibrated by C1-5 was not capable of capturing the increasing trend in the stimulatory effect of IL-8 on the early differentiation rate moving from 10 to 100 ng/ml, reported in study 5 (see Fig 7). According to Fig 8B, the values of the three parameters of n ALP , p 1bie , and α es were notably different between C1-5 and C5. As mentioned earlier, the parameter of n ALP is not relevant in the calculation of the differentiation rate and is, therefore, excluded from further discussion. The parameter of p 1bie , which is associated with IL-1β, is also irrelevant in the formulation of IL-8. Thus, we further examine the parameter of α es , i.e. the sensitivity of the differentiation process to the stimulatory factors, which is directly involved in the calculation of the early differentiation rate (see Eq (3)). It can be seen that C1-5 reports a higher quantity of α es compared to C5 (see Fig 2B). According to Eq (3), a higher quantity of α es increases the overall differentiation rate and results in a faster maturation process. Further evaluation of the model calibrated by C1-5 showed that the simulated maturity already reaches 1, i.e. ultimate value, for the applied concentration of 10 ng/ml, leaving no space for further improvement for higher stimulatory signals, e.g. 100 ng/ml (these results were not provided). This observation implies that the cells used in study 5 are less sensitive to stimulatory signals.

Calibration process
We showed that the employed DE approach results in different sets of parameter values at each run of the calibration process (see Fig 2A). In order to obtain consistent values, we repeated the tuning process at least two hundred times, depending on the study, and used the mean parameter values (see section 3.1). This approach is computationally expensive but results in consistency in the obtained parameter values, which can be used to investigate the differences among different models. Such discrepancy can be due to inherent differences Table 1. Fuzzy logic rules define the osteogenic reaction in response to stimulatory inputs. For the qualitative definition of each linguistic term, refer to Fig 10A. The symbol � indicates that the given conditions must occur simultaneously to produce the given intensity. The symbol~indicates that any choice of one or more from the chosen inputs produces the same intensity. The term 'Not', preceded by a linguistic level, indicates that the rule applied for all except the given level. The terms 'ED' and 'LD' stand for early differentiation and late differentiation, respectively. The rules are given in IF/THEN format in S2 among different experiments such as donor characteristics, culture medium characteristics, and experimental protocol [11], [17], [34].

Sensitivity analysis
In order to study the sensitivity of the model to its parameters, we employed two different approaches of LSSP and SSIP. LSSP captures the effect of large and simultaneous variations in the parameter values, while SSIP examines the impact of small and individual perturbations (see section 4.6). Together, these approaches provide a superior understanding of the system compared to the single approach utilized in our previous study [31]. In the formulation of the present model, we assumed non-linear relationships between the simulated maturity and the differentiation markers (see Eq (6)), which was assumed linear in our previous model [31]. The results of the sensitivity analysis showed n ALP , n ARS , and n OC , which control the degree of non-linearity in these relationships, are among the most significant parameters in the model (see Fig 9). The inferred values of these parameters were higher for the case of C1 to C5 compared to C1-5 (see Fig 8 and S3 Table). This is primarily due to the overfitting phenomenon as the calibrated models of C1 to C5 are too closely aligned to the data of individual studies. For the case of C1-5, the values of these parameters are close to 1 (see S3  Table), showing a linear relationship between the simulated maturity and the differentiation markers. Another significant parameter in the model is β ALP , according to the results of the sensitivity analysis (see Fig 3). This parameter represents the baseline value of ALP secretion when the simulated maturity is zero (see Eq (6)). The estimated values of this parameter lie within 0.4 and 1.4 for different calibration scenarios (see S3 Table). The addition of this parameter to the relationship between the simulated maturity and the differentiation markers is an improvement to the formulation of our previous model, in which the baseline value is assumed zero [31]. The parameters of M t , T d , and α es were also shown as the influential parameters according to the sensitivity analysis (see Fig 3). As elaborated earlier, these parameters, which are directly involved in the formulation of the simulated maturity, were the determinant factors in explaining the discrepancy in the simulation results.

Fuzzy-logic controller as the decision-making center of cell
In the present study, we used FL-based controllers as the deciding center of cells to simulate the process of osteogenic differentiation. FL-based simulations use plain language to describe a system which can potentially help in dissolving technical barriers between simulation and experimental experts easing their involvement in the rapid development of a computer model [26], [35]. Moreover, FL-based models are tractable and interpretable, which makes them a suitable option in investigating and incorporating the experimental datasets which are not in agreement with one another, similar to this study and our previous report [31]. Since the present model receives the signals from the environment and predicts cellular actions, similar to an agent, this descriptive model can serve as a natural basis for the predictive reinforcement model in the future [36].

Limitations, future developments, and concluding remarks
The proposed computer model in this study has several important limitations. Firstly, we were not able to find sufficient information in the literature to define the quality of interactions between different signaling factors in the FL controller, except for the two factors of IL-8 and IL-1β (see Table 1). A similar limitation was reported in our previous model [31] as well as in other computer models in the literature [37][38][39]. Further experiments are required to investigate the potential synergic effects between different cellular inputs, in particular, Mg 2+ ions with the inflammatory cytokines [17]. Secondly, the proposed model in this study, which simulates the behavior of a single cell, is used to explain the average behavior of the cells cultured in the experiments. We defined several parameters to scale the outputs of the cell model to match the empirical observation (see section 4.5). However, the current model is not able to capture the heterogeneity in the behavior of individual cells. In the next step, the proposed cell model will be incorporated into an agent-based model to address this issue. This can potentially improve the simulation results but will require significantly higher computational costs, i.e. in the order of thousand times. Thirdly, the regulatory effect of substrate stiffness as an important factor in guiding osteogenesis is not incorporated in the present model [40]. This parameter will be included in our future studies. Similarly, transforming growth factor-beta (TGF-β) and bone morphogenic protein (BMP), which are two crucial growth factors in the osteogenic differentiation of MSCs, will be adopted from our previous model [31]. To do so, the quality of interaction between Mg 2+ ions and these growth factors need to be studied first as Mg 2+ ions are also reported to influence the activation of transforming growth factor-beta (TGF-β) and bone morphogenic protein (BMP) signaling [8]. Lastly, in this study, we simulated the influence of the inflammatory cytokines, primarily produced by macrophages, on MSC osteogenic differentiation. However, the cross-talk between MSC and macrophage is a double-sided process; macrophage polarization and cytokine production are also regulated by MSC. In the next steps, we will address this problem and complete the circle of MSC-macrophage interaction.

Materials and methods
In this section, we first explain the complete process of the development of the FL controller (see Fig 10). The bioregulatory roles of the cellular inputs in the differentiation activities are explained in the following subsection together with the process of converting the qualitative knowledge into the machine-readable algorithm, termed fuzzification. After, the fuzzification process is explained for the cellular outputs. Then, the process of fuzzy inference and defuzzification is elaborated. The outputs of the FL controller require post-processing in order to define the actual rates of early and late differentiation rates. This process is further elaborated in the subsequent section. Finally, we elaborate on the empirical data used for the calibration and sensitivity analysis with the technical details of each process in the following sections. We used several software and packages to develop the current model, which can be found S1 Text. The source code of the present model can be found online [41].

Cellular inputs and fuzzification process
The first step in the development of the FL controller is to define and fuzzify the cellular inputs (see Fig 10). In this section, we elaborate on the bioregulatory roles of the five signaling factors in the early and late differentiation processes using the information available in the literature. Then, FL membership functions are defined according to this information.
Mg 2+ ions are shown to affect osteogenic differentiation depending on the concentration and the state of cell differentiation [42][43][44][45]. A Mg 2+ ion concentration of 0.8 mM, which is commonly used as the minimal essential Medium-MEM in cell culture experiments, is considered as the physiological concentration. A concentration below 0.8 mM is shown to delay the early differentiation process of MSCs [17]. Mg 2+ ions within the concentration range of 2-10 mM is reported to promote early differentiation rate [6], [43], [45][46][47], while the concentration of Mg 2+ ions above 1.8 mM is reported inhibitory for late differentiation process [6], [43], [45]- [47]. In addition, a concentration over 20 mM is shown to compromise cell viability [47][48][49]. To account for these observations, we create six membership functions to define the cellular input of Mg 2+ ions as shown in Fig 10A. The parameter of p ms is defined to mark the peak occurrence of the Stimulatory effect. The parameter of p md marks the beginning of the Inhibitory effect. The membership function of Ineffective is defined as the transition state from Stimulatory to Inhibitory state.
IL-10 is a key anti-inflammatory cytokine in the bone regenerative process including osteogenic differentiation [33]. The mediatory effects of IL-10 are shown to significantly depend on the applied concentration as well as the length of application [13], [33] [50]. Chen et al. [33] showed that IL-10 has a dual effect on osteogenic activities; IL 10 in low concentration activates p38/MAPK signaling pathway and stimulates osteogenic activities, while higher concentrations of IL-10 inhibit p38/MAPK signaling by activating NF-kB and consequently downregulating osteogenesis. They found that IL-10 has the highest stimulatory effect at the concentration of 0.1 ng/ml, while the concentration of 10 ng/ml showed a significantly inhibitory effect on osteogenesis [33]. However, Valles et al. [13] showed an increasing stimulatory effect of IL-10 on osteogenesis by increasing the concentration from 0.1 to 10 ng/ml. They also experimented the importance of application time and concluded that while IL-10 can be proosteogenic after short-term treatment (48 hours), its continuous application can produce inhibitory effects [13]. Therefore, in the definition of the FL controller, we set the application time of 48 hours to differentiate short-and long-term treatment periods. Accordingly, we define two separate sets of functions to fuzzify the cellular input of IL-10 depending on the application time. In the controller that simulates IL-10 for the application of less or equal to 48 hours, we define a constant increase in the stimulatory effect of IL-10 from 0 to 10 ng/ml (see Fig 10A). However, for the case where the application time exceeds 48 hours, the promotory effect of IL-10 decreases by exceeding 0.1 ng/ml (see Fig 10A), in compliance with the findings of Chen et al. [33].
TNF-α is a primary cytokine in the inflammatory reaction which plays an important role in osteogenic differentiation [51]. The regulatory effects of TNF-α on osteogenesis have been shown to be concentration-dependent; while low concentration of TNF-α is assumed stimulatory, its high concentration has an inhibitory effect on osteoblastic differentiation [13], [52]. According to the results of Valles et al. [13] and Glass et al. [52], we define the stimulatory concentration of TNF-α as 1 ng/ml and the inhibitory concentration as 100 ng/ml [52]. TNF-α at a concentration of 10 ng/ml is shown to have neither stimulatory nor inhibitory effect on osteogenesis [13], [52]. Thus, the concentration of 10 ng/ml is defined as the neutral state with an ineffective role in osteogenesis (see Fig 10A).
IL-8 is traditionally classified as a pro-inflammatory cytokine with the main role of recruiting the inflammatory cells to the injury site [17] [53]. In addition, IL-8 is shown to play an important role in the commitment of MSCs to bone cells [54] [55]. Qiao et al. [17] showed that IL-8 is especially important when the inflammatory response is mediated by Mg 2+ ions; IL-8 is the primary player in promoting osteogenesis when macrophages are treated by Mg 2+ ions [17]. They also showed that the regulatory effect of IL-8 on osteogenesis is dose-dependent and substantially increases within the concentration of 0 to 100 ng/ml [17]. To account for these observations, we define three membership functions to fuzzify the cellular input of IL-8 (see Fig 10A), where the stimulatory role of IL-8 rapidly increases by an increase in its concentration. The parameter of p 8f marks the peak occurrence of Favorable condition which is an intermediate state toward Stimulatory level.
IL-1β is a major pro-inflammatory cytokine with an important role in the simulation of osteoclastogenesis [56]. IL-1β is also shown to increase osteogenic activities of MSCs, especially when it is administrated in a low concentration within 1 to 10 ng/ml [17]. To account for these observations, we define three membership functions to fuzzify IL-1β (see Fig 10A), where two parameters of p 1bs and p 1bie mark the peak and end of the stimulatory effect. In addition, IL-1β is reported to be an antagonizer of IL-8 osteogenic effects [17]; Qiao et al. [17] showed that IL-1β substantially reduces the secretion of the osteogenic biomarkers stimulated by IL-8. Thus, in the present FL controller, IL-1β has a dual role in osteogenetic response; on the one hand, it stimulates osteogenesis; on the other hand, it hampers the stimulatory role of IL-8.

Cellular outputs and fuzzification process
The outputs of the fuzzy controller are early and late differentiation rates. These values are defined in linguistic formats using the set of membership functions given in Fig 10B. The early differentiation process is assumed to occur in six different rates within the range of 0 and 1, with 0.5 representing the physiological rate. The parameters of p es , p ef , and p evf are defined to mark the intensities of Slow, Fast, and Very fast levels, accordingly. The late differentiation process is fuzzified using five membership functions shown in Fig 10B, with 0.5 representing the physiological rate (see Fig 10B). The parameter of p ls and p lf are defined to mark the Slow and Fast membership levels. Gaussian membership functions (GMFs) with the activation value of 1 and the sigma value of 0.05 are employed for the fuzzification of the output variables. GMFs offer multiple advantageous such as smoothness and concise notation as well as superior reliability and robustness of the system [57][58].

Fuzzy inference and defuzzification process
The set of rules given in Table 1 is used to determine the differentiation rates in response to the cellular inputs. A Mamdani-type FL controller is used for the calculations. Since multiple rules can be triggered simultaneously, the centroid technique is used to determine the final output of the controller. In this technique, the area bound by the activation degree is summed over different fuzzy sets, and the center of the sums is calculated to represent the final outcome (see Fig 10C) [59].

Simulation of the differentiation process
Similar to our previous study [31], we define the factor of maturity to mark the degree of progression along the line of osteogenic differentiation. Maturity holds a value in the range of 0 and 1 and is calculated for a given time (T) using early and late differentiation rates (r e and r l , respectively) as, where T e is time required for the early maturation process.
where M t is the early maturation threshold, and T d is the time required for MSC to fully differentiate to osteoblast. M t is a value between 0 and 1, marking the end of early differentiation process.
In Eq (1), r e and r l are the scaled versions of f e and f l (fuzzy outputs), respectively, where r 0 is the physiological rate of osteogenic differentiation, calculated as 1/T d , where S is a function that scales the outputs of fuzzy controller, which are initially between 0 and 1, where x is the input from the FL controller. α s,k and α i,k are the scaling coefficients for the stimulatory and inhibitory effects, respectively, k 2 {e, f}.

Empirical data
The empirical data obtained from 5 sets of published experiments are used to estimate the free parameters of the present model (see S1 Table for the summary of the experiments). These experiments evaluate the bioregulatory effect of different Mg 2+ ions as well as different inflammatory cytokines on the osteogenic differentiation of MSCs. Three markers of ALP, OC, and ARS are used to study the differentiation process. ALP is commonly used as an early differentiation marker, while OC and ARS are recognized as the markers of the late differentiation phase [40], [60]. These factors are correlated with maturity, which is the simulated indicator of differentiation progress in the present study, where x is the conditioned value of maturity, y i is the quantity of the biomarker with i 2 {ALP, OC, ARS}, β i is the baseline of the marker, n i is the degree of nonlinearity, and k i,j is the correction factor applied for each experiment with j 2 1:5. For ALP, x is calculated as, In Eq (7), we assume that ALP correlates with maturity until the early maturation threshold and stays constant afterward, accounting for the fact that ALP is the early differentiation marker. OC and ARS, as the late marker of differentiation, correlate with maturity in its whole range, i.e. x = maturity. In Eq (6), β i takes into account that the differentiation markers can be detected even at the beginning of an experiment where the simulated maturity is zero [33]. The parameter of n i is defined to catch the non-linear correlation between maturity and the differentiation markers. Lastly, k i,j takes into account the variations in the measurements of the markers across different experiments, e.g. differences in the reported units.

Calibration process and sensitivity analysis
The overview of the calibration process is given in Fig 11A. DE is a metaheuristic stochastic search algorithm based on an evolutionary process which improves the candidate solution by an iterative approach [61]. Although DE is considered a global optimization algorithm, the estimated parameter values can be different at each run. To overcome this issue, we execute multiple DE runs and use the posterior distributions to infer the parameter values (see Fig 11A).
The overview of the sensitivity analysis is given in Fig 11B. Two approaches of LSSP and SSIP are employed to study the sensitivity of the model to its parameters. In LSSP, one or several parameters are simultaneously perturbed with a magnitude of 50 percent around the inferred values. Fractional factorial design (FFD) and the analysis of variance (ANOVA) are used for sampling and analysis of this approach, which are explained in detail in our previous publication [31]. In SSIP, the parameters of the model are perturbed one at a time, with a magnitude of 15 percent around the inferred values.