Modeling the Effects of Macrophages on Bone Fracture Healing

A new mathematical model is presented to study the effects of macrophages on the 1 bone fracture healing process. The model consists of a system of nonlinear ordinary differential 2 equations that represents the interactions among classically and alternatively activated macrophages, 3 mesenchymal stem cells, osteoblasts, and proand anti-inflammatory cytokines. A qualitative 4 analysis of the model is performed to determine the equilibria and their corresponding stability 5 properties. Numerical simulations are also presented to support the theoretical results and to monitor 6 the evolution of a broken bone for different types of fractures under various medical interventions. 7 The model can be used to guide clinical experiments and to explore possible medical treatments that 8 accelerate the bone fracture healing process either by surgical interventions or drug administrations. 9


Introduction
Bone fractures are becoming a big worldwide problem due to their high frequency and surgical complications.Globally, osteoporosis causes more than 8.9 million fractures every year affecting 50% of women and 25% of men over age 50 [1][2][3].Furthermore, 10-15% of bone fractures do not heal or take longer to heal [4,5].Surgical complications, disabilities, and high morbidity rates often occurs from severe traumas and immune-compromised-fractured people [1,5].The most dangerous trauma is the hip fracture with mortality rates up to 20-24% in the first year after the fracture and greater risk of dying may persist for at least 5 years afterwards [6].Traffic accidents, the number one killer of young people and the major causes of fractures, are expected to be one of the top three causes of disabilities by 2020 [7,8].In addition, medical care costs for osteoporotic bone fractures are expected to be over US$25 billions by 2025; due, in part, to the expensive treatments and the prolonged hospitalization and rehabilitations [1,9].It is essential to have a better understanding of the bone fracture healing process in order to prevent unsuccessful healing and to develop optimal fracture union treatments.
Although significant improvements have been made in the experimental and mathematical modeling of the bone fracture healing process over the last twenty years, however the optimal conditions for bone repair are still unclear [10][11][12].Treatments based on anti-inflammatory cytokines, such as the cytokine-specific agents that block the pro-inflammatory cytokines productions, have exhibited promising clinical results and have led to intense orthopedic research activities [4,5,[12][13][14][15][16].
Recently, a mathematical model based on the interactions among macrophages, mesenchymal stem cells (MSCs), and osteoblast was developed to study the regulatory effects of two generic pro-and anti-inflammatory cytokines during the early stages of bone fracture healing [17].To our knowledge, it was the first attempt to incorporate the macrophages interactions in the modeling of the bone fracture healing process.The mathematical model revealed that high concentrations of pro-inflammatory cytokines negatively affect the healing time of a fracture and that the administration of anti-inflammatory cytokines can accelerate the healing time in a dose-dependent manner.However, the model assumed that the only source of anti-inflammatory cytokines is given by the MSCs, which may not be enough to promote and correctly represent the complex pattern of bone fracture healing formation.Therefore, it is better to incorporate the other sources of anti-inflammatory cytokines during the bone fracture healing process, such as the delivered from the macrophages [1,5,18].
In this paper, the mathematical model developed in [17] is extended to separately incorporate the two different phenotypes of macrophages: classically and alternatively activated macrophages, as they have distinct functions during the tissue healing [18,19].Classically activated macrophages release high levels of pro-inflammatory cytokines, including the TNF-α and IL-1β which exhibit inhibitory and destructive properties in high concentrations [13,19].In contrast, alternatively activated macrophages are characterized with the secretion of the anti-inflammatory cytokines such as the IL-10 which increase their phagocytic activities, mitigate the inflammatory responses, promote growth and accelerate fracture healing [1,5,18,19].This extension leads to a more realistic model by incorporating the different phagocytic rates and the separate production of the pro-and anti-inflammatory cytokines by the two types of macrophages [18,20].The model can be used to investigate potential therapeutic treatments based on the use of anti-inflammatory cytokines, stem cells, and macrophages, suggesting possible ways to guide clinical experiments and bone tissue engineering strategies [18,19].
The organization of the paper is as follows: Section 2 discuses the cellular and molecular interactions that occur during the bone fracture healing process.The simplifying assumptions are presented in Section 3. In Section 4, a system of nonlinear ordinary differential equations is introduced to describe the fundamental aspects of the bone fracture healing process during the resolution of inflammation and bone repair.The stability analysis of the system is presented in Section 5. Section 6 demonstrates the functionality of the model by numerically simulating the progression of the bone fracture healing process under normal and pathological conditions.The discussion and future work are presented in Section 7.

Biological Background
Bone fracture healing is a complex regenerative process that involves the participation of different cell types including the immune system and mesenchymal lineage cells [21].Their interactions and functions are strongly regulated by molecular and mechanical stimuli [19,22].Particularly, at the beginning of the healing process, cytokines either have positive or negative effects on the cellular functions depending on the influence of other cytokines, concentration, and exposed time [23][24][25].
Cytokines are functionally classified into pro-inflammatory and anti-inflammatory families.
Pro-inflammatory cytokines such as the tumor necrotic factor-α (TNF-α) activate the immune system defense to kill bacteria and fight infections, while anti-inflammatory cytokines inhibit pro-inflammatory synthesis and activate the mesenchymal lineage cellular functions [5].The interleukin-10 (IL-10) is one of the most potent anti-inflammatory molecules that inhibit the pro-inflammatory production [5,26].
The pro-inflammatory cytokines concentration during the bone fracture healing process is mainly delivered by necroses of cells and by the inflammatory immune cells, such as monocytes and neutrophils, that arrive to the injury site in response to the trauma [4].These pro-inflammatory profiles, which include the TNF-α, lead to an acute inflammation observed in the first 24 hrs [11,27,28] after injury.Monocytes migration mostly occurs during the beginning of the acute period, when monocytes also differentiate into macrophages to resolve the inflammation.Once this differentiation starts, the influx of the inflammatory cells ceases and they die out [29].
During the resolution of inflammation, macrophages increase their population by migration and they activate to their classical and alternative phenotypes accordingly to the cytokines stimuli [19,30].
The two phenotypes can also shift between each other during this process [31,32].Macrophages have the capability to release both pro-and anti-inflammatory cytokines through their different activation [31].Classically activated macrophages release high concentration of pro-inflammatory cytokines, such as the TNF-α, and low levels of anti-inflammatory cytokines [31] in responses to their engulfing functions.Alternative macrophages secrete high level of the Il-10 and low levels of TNF-α as they continue with the clearance of debris and the modulation of inflammation [31].The correct balance of TNF-α during bone repair is necessary for successful fracture healing.High levels of TNF-α induce a chronic inflammation and gradual destruction of cartilage and bone tissue [25], while the absence of TNF-α impairs fracture healing [13,15].
In addition, during the resolution of inflammation, MSCs arrive to the injury site, activate their immune-modulation functions by releasing the IL-10, and proliferate or differentiate into fibroblasts, chondrocytes, and osteoblasts [5,33].Fibroblasts and chondrocytes proliferate and release the fibrinous/cartilagenous extracellular matrix, which fills up the fracture gap [33,34].Osteoblast cells proliferate, synthesize the new woven bone, and differentiate into osteocytes or die out [22,33].
During the last stage of the bone fracture healing process, the fibrocartilage and the woven bone are constantly removed and replaced by a functional bone [35].This process is referred to as bone remodeling and consists of a systematic tissue degradation and production by osteoclasts and osteoblasts, respectively.Bone remodeling is a slow process that can take months to years until the bone recovers to its pre-injury state [11].
The inflammation is considered resolved when debris are eliminated, activated macrophages emigrate to the lymphatic nods to die, and inactivated macrophages return to their normal density [29].These evens are observed after two weeks from the beginning of the healing process [35,36].
Fibrinous/cartilaginous tissue production is observed in the first 3 days, it peaks in about 10 to 12 days, and its removal starts as early as 21 days [33].Approximately at 28 to 35 days, osteoclasts populate the tissue and the removal of the fibrocartilage is substantially observed [35].The fracture healing outcome is considered a delayed union if the fibrous/cartilaginous tissue is not removed completely in about 3 to 4 months after the injury, while it is considered a nonunion if no functional bone is obtained in 6 months after the trauma [37].

Modeling Assumptions
The most important effects of macrophages on the bone fracture healing process are observed during the inflammatory and repair phases of the bone fracture healing process [17].During the inflammatory phase, macrophages modulate and resolve the inflammation while during the repair phase macrophages provide an optimal environment for the cellular proliferation, differentiation, and tissue production.The primary cells during the inflammatory and repair phases of the bone fracture healing process are debris (D), unactivated macrophages (M 0 ), classical macrophages (M 1 ), alternative macrophages (M 2 ), MSCs (C m ), and osteoblasts (C b ).It is assumed that the cellular functions are regulated by the tumor necrotic factor-α (c 1 ) and the interleukin-10 (c 2 ) cytokines.It is also assumed that the regenerative process is given by the production of two extracellular matrices: the fibrocartilage (m c ), and the woven bone (m b ).These biological system interactions are depicted in Figure 1.The variables represent homogeneous quantities in a given volume.
In Figure 1 It is assumed that unactivated macrophages M 0 do not release cytokines and do not engulf debris.
It is also assumed that the population of M 0 increases in size proportionally to the debris concentration up to a maximal value of M max [30].The only source of activated macrophages, M 1 and M 2 , is M 0 .
Even though both phenotypes of activated macrophages have the abilities to release both pro-and anti-inflammatory cytokines, it is assumed that only M 1 deliver c 1 and M 2 deliver c 2 , as those are the major cytokines for each phenotype [38].M 0 activate to M 1 under the c 1 stimulus, while they activate to M 2 under the c 2 stimulus.M 1 and M 2 macrophages do not de-differentiate back to the M 0 macrophages [39]; and are able to switch phenotypes at a constant rate [32].The accumulation of macrophages at the injury site is modeled by its recruitment due to inflammation, which is assumed to be proportional to the debris density.It is assumed that MSCs differentiate into osteoblasts at a constant rate.MSCs synthesize the fibrocartilage, while osteoblasts synthesize the woven bone.It is also assumed that only the fibrocartilage is constantly removed by the osteoclasts.The density of the osteoclasts is assumed to be proportional to the density of the osteoblasts.The two types of tissue cells increase their populations by proliferation in a logistic growth fashion [33].It is also assumed that there is no recruitment of MSCs and osteoblasts.

Model Formulation
The process of bone fracture healing is modeled with a mass-action system of nonlinear ordinary differential equations.Following the outlined biological assumptions and the flow diagram given in Peer-reviewed version available at Math.Comput.Appl.2019, 24, 12; doi:10.3390/mca24010012 Equation ( 1) describes the rate of change with respect to time of the debris density, which decreases proportionally to M 1 and M 2 .The engulfing rate R D is modeled by a Hill Type II function to represent the saturation of the phagocyte rate of macrophages [38,40]: Equation ( 2) describes the rate of change with respect to time of the undifferentiated macrophages density.It increases because of migration and decreases by differentiating into M 1 and M 2 or by a constant emigration rate.It is assumed that M 0 migrate to the injury site proportionally to D up to a maximal constant rate, k max , [26,31]: where The differentiation rates of M 0 into M 1 and M 2 are stimulated by the cytokines accordingly to a Hill Type II equations, respectively [32]: Equation ( 3) describes the rate of change with respect to time of M 1 , which increases when M 0 activate to M 1 and M 2 shift phenotype; and decreases by emigration and when M 1 shift phenotype.Similarly, Equation (4) describes the rate of change with respect to time of M 2 .Equations ( 5) and ( 6) describes the rate of change with respect to time of c 1 and c 2 .Here, k 0 , k 1 , k 2 , and k 3 are the constant rates of the cytokine productions and d c 1 and d c 2 are the cytokine constant decay rates.The inhibitory effects of the anti-inflammatory cytokines are modeled by the following functions [32]:  7) describes the rate of change with respect to time of C m , which increases by cellular division up to a constant-maximal carrying capacity, K lm , and decreases by differentiation.The total MSCs proliferation rate is modeled by [28]:

Preprints
where in the absent of inflammation, c 1 = 0, MSCs proliferate at a constant rate k pm .However, when there is inflammation, c 1 > 0, the proliferation rate of MSCs increases or decreases according to the concentration of c 1 , i.e., high concentration levels of c 1 inhibit C m proliferation while low concentration levels of c 1 accelerate C m proliferation.The differentiation rate of C m is inhibited by c 1 , which is modeled by the following function: Equation ( 8) describes the rate of change with respect to time of C b .It increases when MSCs differentiate into osteoblasts or when osteoblasts proliferate.It decreases at a constant rate d b when osteoblasts differentiate into osteocytes.The osteoblasts proliferation rate is inhibited by c 1 , which is modeled by the following function: Equations ( 9) and ( 10) describe the rate of change with respect to time of the fibrocartilage and woven bone, where p cs and p bs are the tissue constant synthesis rates and q cd1 , q cd2 , and q bd are the tissue degradation rates, respectively [33].

Analysis of the Model
The analysis of Model ( 1)-( 10) is done by finding the equilibria and their corresponding stability properties.An equilibrium is a state of the system where the variables do not change over time [41].
Once the equilibria are identified, it is important to determine the behavior of the model near equilibria by analyzing their local stability properties.An equilibrium is locally stable if the system moves toward it when it is near the equilibrium, otherwise it is unstable [41].Therefore, the equilibria provide the possible outcomes of the bone fracture healing process and their corresponding stability properties define the conditions under which a particular healing result occurs.
System (1)-( 10) has the following three biologically meaningful equilibria of the vector form The existence conditions for the three equilibria are summarized in Table 1 and their stability conditions are summarized in Table 2 and proved in Appendix A.
The existence conditions listed in Table 1 arise from the fact that all biologically meaningful variables are nonnegative.Therefore, the existence condition for E 0 requires the steady state tissue densities to be either zero or any positive number.For E 1 , the existence condition arises from the requirement that the steady state density of C b must be greater than zero, which implies that the proliferation rate of osteoblasts must be greater than their differentiation rate, i.e., k pb > d b .
Similarly for E 2 , the existence condition arises from the requirement that the steady state density for C m must be greater than zero, which implies that the proliferation rate of MSCs must be greater

Equilibrium Points Existence Conditions Meaning
E 0 = (0, 0, 0, 0, 0, 0, 0, 0, The stability conditions of each biologically feasible equilibrium are listed in Table 2 and determined from the eigenvalues of its associated Jacobian matrix, see Appendix A, as follows: ) which implies that the differentiation rates of the MSCS and osteoblasts are greater than or equal to their proliferation rates, respectively.
The steady-state E 0 represents a nonunion.In this case, the inflammation is resolved since the first five entries of E 0 are zero; however, the repair process has failed since the osteoblasts and osteoclasts have died out before the beginning of the remodeling process.Hence, the tissue densities, m * c 0 and m * b 0 , can be any two positive values smaller than their maximal densities, p cs /q cd1 and p bs /q bd , respectively (see Theorem A1).

Numerical Results
The proposed new model ( 1)-( 10) is used to investigate the evolution of a broken bone under normal and pathological conditions during the first 21 days after trauma.Table 3 summarizes the baseline parameter values and units for the numerical simulations.These values are estimated in a qualitative manner from data in other studies [30,38,39] and are based on murine experiments with healthy mice having a moderate fracture (a broken bone with a gap size less than 3mm) [33,42].The bone fracture healing process for humans involves the same cells, cytokines, and qualitative dynamics, differing only in the number of cells, concentrations, and the length of time it takes for a full recovery.
First, a set of numerical simulation results is presented to compare two mathematical models of the bone fracture healing process that incorporate macrophages: the model developed in [17]

Comparison of existing models
The model developed in [17] and the present mathematical model ( 1)-( 10) are compared when D(0) < a ed = 4.71 × 10 6 , i.e., the initial debris concentration is below the half-saturation of debris.In this case, the macrophages' digestion rate increases approximately linearly with respect to the debris population, as it is assumed in model [17].The same parameter values are used in both models (Table 3), with k e 1 = 11, k e 2 = 48, k 2 = 3.72 × 10 −6 , and k 3 = 8 × 10 −6 .
Figure 2 shows the numerical evolutions of the tissues' production when D(0) = 2 × 10 6 .In all simulations, we refer to fibrocartilage and woven bone as cartilage and bone, respectively.The production of cartilage m c and bone m b given by the present model is much more realistic than the production given by the model developed [17], since, according to the experimental data, the cartilage production peaks to its maximal density of around 1g/mL about 10-12 days after trauma and a significant bone tissue production is observed after the second week [44].Figure 3 shows the qualitative behavior of E 1 for the macrophages, debris, TNF-α, and IL-10 densities, with the inflammation being resolved in about 40 days.The top-left plot of Figure 3 shows the temporal evolution of M 0 (dashed lines), M 1 (dotted lines), and M 2 (solid lines).It can be observed that M 1 first peaks to its maximum value, which is then followed by M 2 .Similar sequences of transitions of first M 1 and then M 2 are commonly observed in normal tissue healing conditions [5,39].Figure 4 shows the qualitative behaviors of E 1 for the MSCs, osteoblasts, cartilage, and bone densities.Here, the MSCs density decays to zero over the time, while the osteoblasts maintain a constant density below their carrying capacity K lb = 1 × 10 6 .In addition, the bottom plots of Figure 4 show that the cartilage is eventually degraded by the osteoclasts and the bone achieves its maximum density of 1 ng/mL.Therefore, E 1 exhibits the temporal progression of a successful bone fracture healing.Figure 5 shows the qualitative evolution for the MSCs, osteoblasts, cartilage, and bone densities for E 0 (solid lines) and E 2 (dotted lines).Since the temporal evolution of macrophages, debris, and cytokines densities in E 0 and E 2 are similar to those for E 1 showed in Figure 3, then they are omitted here.It can be observed in Figure 5 that the two cellular densities in E 0 , MSCs and osteoblasts, decay to zero over the time, with the osteoclasts failing to degrade the cartilage; which results in nonunion.
Mathematically, this case occurs when osteoblasts proliferate at a rate lower than their differentiation rate, i.e., k pb < d b .In practice, this scenario is commonly observed in advanced-age patients whose MSCs and osteoblast cells decrease their capability to proliferate and differentiate [1].On the other hand, the two cells and the two tissues in E 2 remain at positive constant values (Figure 5), but the final fracture healing outcome is still a nonunion.Here, the osteoclasts again fail to degrade the cartilage [1], even though the bone has achieved its maximum density of 1 ng/mL.Therefore, in this case, migration of osteoclasts must be enhanced through surgical interventions in order to achieve a successful bone repair [33].

Evolution of the healing process for different types of fractures
In this section, the model is used to monitor the evolution of a successful repair (Table 3) for different types of fractures.In healthy individuals, the simple, moderate, and severe fractures are correlated with the debris densities [45,46].Therefore, the initial debris concentration is set to D(0) = 3 × 10 5 , D(0) = 5 × 10 7 , and D(0) = 2 × 10 8 , for a simple, moderate, and severe fracture, respectively.Figure 6 shows that the tissues production is a slow process for a simple fracture, since the cartilage and bone densities are less than the corresponding tissue densities for moderate and severe fractures.Slow healing process is commonly observed in micro-crack healing [45].Furthermore, there is less cartilage formation over time in simple fractures [46].For a moderate fracture, the maximal production of the cartilage is observed around 10 days followed by a significant degradation, while the bone tissue production occurs after the first week.For a severe fracture, Figure 6 shows that there is a delay in the two tissues production compared with those given by moderate fractures, with the peak of the cartilage and bone productions observed at around day 16.

Immune-modulation therapeutic treatments to accelerate bone fracture healing
The administration of anti-inflammatory drugs and the injection and/or transplantation of MSCs and macrophages are two of the clinical trials that have been implemented in orthopedics to stimulate and accelerate bone fracture healing [5,15].In this section, Model (1)-( 10) is used to explore these possible therapeutic treatments to accelerate the healing of a broken bone under normal and pathological conditions such as severe fractures, advanced age, and senil osteoporosis [1].

Administration of anti-inflammatory drugs at the beginning of the healing process
A set of numerical simulations is presented to investigate the effect of the administrations of anti-inflammatory cytokines at the beginning of the healing process in healthy individuals and also in immune-compromised patients.In each case of the numerical simulations, c 2 (0) = 0, 10 and 100 ng/mL.Figure 7 shows that the administration of c 2 in the simple fracture slows down both the cartilage and bone productions.Figures 8 and 9 show that the administration of c 2 in the moderate fractures improves the tissues evolution but in a dose-dependent manner.On one hand, when D(0) = 2 × 10 7 the administration of c 2 has either a positive or negative effect on the two tissue productions.The administration of 10 ng/mL of c 2 enhances the early production of cartilage and increases the bone synthesis, while the administration of 100 ng/mL of c 2 results in the opposite effect.On the other hand, when D(0) = 5 × 10 7 the administration of c 2 enhances the earlier cartilage production and improves the synthesis of the bone for both concentrations, with 10 ng/mL being the optimal of the two doses.Next, the model is used to implement the administration of anti-inflammatory drugs under different pathological conditions.First, severe fractures in immune-compromised individuals are simulated by using the following parameter values: D(0) = 2 × 10 8 and k max = 0.0015, since in fractures of such individuals there is an increase in the accumulation of debris [46] and a decrease in the macrophages migration rate [47].Second, the following parameter values are used: and k 1 = 9 × 10 −6 to simulate bone fracture healing in aging individuals, since in this case, the macrophages phagocytic rate decreases and there is an increase of pro-inflammatory cytokine synthesis by M 1 [1,48].Finally, c 1 (0) = 100, k pm = 0.2, d m = 0.5, k pb = 0.16, and d b = 0.15 are used to simulate the healing process for an senil osteoporotic fracture, since in this case a high level of pro-inflammatory cytokines is observed and the MSCs and osteoblast functions decrease [1].Figures 10,11,and 12 show that the administration of anti-inflammatory cytokines under the above three under different pathological conditions always improve tissues productions; where the optimal dose of c 2 for both the advanced-age individuals and senil osteoporotic fractures is 10 ng/mL.

Cellular therapeutic interventions under immune-compromised conditions
Additions of MSCs to the injury site through injection and/or transplantation have been used in practice to stimulate and augment bone fracture healing [5].Another cellular intervention is the scaffold implants, where undifferentiated macrophages and MSCs are co-cultured together, and cytokines are slowly released to stimulate M 2 activation [1].The parameter values used in the numerical simulations that explore these possible therapeutic treatments are the same as in Subsection 6.4.1.
For severe fractures with immune-compromised conditions, the use of scaffold implants is simulated through a fast M 2 activation, i.e., k 02 = 0.3 and k 12 = 0.075, and also an increase in the  C m and M 0 densities, i.e., M 0 (0) =

Discussion and Conclusions
A new mathematical model was introduced to mathematically and numerically study the complexity of the molecular and cellular interactions during the bone fracture healing process.
The model examined the macrophages functions and their interactions with the tissue cells during the inflammation and repair phases of the healing process.Classically and alternatively activated macrophages were incorporated to mathematically represent their capabilities to modulate and resolve the inflammation.It also included the macrophages abilities to regulate the tissue cellular functions by the delivery of pro-and anti-inflammatory cytokines.In the new model, the resolution of the inflammation is initiated with the activation of the macrophages into their classical phenotype.The classically activated macrophages deliver the pro-inflammatory cytokine TNF-α as they engulf debris.
Then the alternatively activated macrophages and the MSCs modulate the inflammation by releasing the anti-inflammatory cytokine IL-10.Finally, the classically activated macrophages remove the remaining debris.The model also incorporated the different engulfing rates of activated macrophages, the saturation rates of phagocytes, and the maximal density of macrophages at the injury site thus allowing a better understanding of the interplay between macrophages and tissue cells during the bone fracture healing process.
The mathematical analysis revealed that there are three feasible fracture healing outcomes.Two of the outcomes represent a nonunion healing: one is the case when the cells deactivate or die out before the healing process finishes up and the other is the case when the tissue cells remain constant but the osteoclasts fail to completely remove the cartilage.The third outcome represents a successful healing, where the osteoblasts and osteoclasts are constantly producing and removing the woven bone.
The stability conditions of each outcome can be used to biologically explain why the fracture healing fails as well as to design therapeutic interventions to stimulate or accelerate the healing process.The new mathematical model allowed a variety of different types of numerical simulations to be performed quickly and cost effectively.It was used to monitor the progression of the healing of a broken bone as well as to predict the final outcome of the healing process.In particular, it was used to numerically simulate the administration of anti-inflammatory drugs to improve the bone fracture healing process.It was found that the administration of anti-inflammatory cytokines fails to accelerate the healing process in simple fractures, while it accelerates the healing process in moderate fractures depending on the cytokine concentrations, and always improves the healing process in severe fractures.Such results have been also clinically observed when corticosteroids and nonsteroidal anti-inflammatory drugs (NSAIDs) are administered in bone fractures [15].Therefore, based on the model findings, the concentration of debris must be carefully considered when administering anti-inflammatory drugs to enhance the fracture healing process [46].The model was also used to explore other potential cellular therapeutic approaches, such as the MSCs injection and transplantation.
It was found that such treatments can also improve the healing time of a broken bone, especially in immune-compromised patients.The model can also be easily adapted to other therapeutic approaches, such as the administration of different anti-inflammatory drugs, suggesting a variety of possible ways to guide clinical experiments and bone tissue engineering strategies. where Therefore the corresponding characteristic polynomial associated with J(E 1 ) is given by the product of the characteristic polynomials associated with each submatrix [50]: Therefore, the eigenvalues of J(E 0 ) are negative for the variables M 0 , M 1 , M 2 , c 1 , c 2 , C m , and C b and are equal to zero for D, m c , and m b .Since D (t) ≤ 0 for all the variables in the system (1)-( 10) and (D * , 0, 0, 0, 0, 0, 0, 0, m c , m b ) with D * = 0 is not an equilibrium point, then the solutions of the system (1)-( 8) are attracted to the set A = {(0, 0, 0, 0, 0, 0, 0, 0, m c , m b ) : m c ≥ 0, m b ≥ 0}.
Equations ( 9) and (10) imply that m c ≤ 0 and m b ≤ 0 for all m c > p cs /q cd 1 and m b > p bs /q bd .Therefore, the set B is a local attractor set of A [49].
Next, let us consider the case when k pm = d m and d b = k pb .Here, the eigenvalues of J(E 0 ) are the same as above except those associated with C m and C b , which are equal to zero.Therefore, in this case, by considering the second order approximations of the right hand sides of Equations ( 7) and ( 8), instead of just the first order approximations, and using similar arguments as above, proves that the set B is a local attractor set of A.
Proof of Theorem A2.The Jacobian matrix corresponding to the point E 1 is given by the following lower triangular block matrix where J 2 (E 1 ) has the same expression as J 1 (E 0 ) defined in Theorem A1 and  Proof of Theorem A3.The Jacobian matrix corresponding to the point E 2 is given by the following lower triangular block matrix where  are positive, then all the eigenvalues of J 1 (E 2 ), J 2 (E 2 ), J 3 (E 2 ) are negative except for the eigenvalue associated to D which is equal to zero.Therefore, since D ≤ 0 for all the variable system, then E 2 is locally stable.
, the cellular dynamics are represented by the circular shapes and solid arrows.The molecular concentrations and their production/decay are represented by the octagonal shapes and dashed arrows.The pro-and anti-inflammatory cytokines activation/inhibition effects on the cellular functions are represented by the dotted arrows.Removal of debris and the negative effect among the variables are represented by the dot-ending dotted arrows.

Figure 1
Figure 1 yields the resulting system of equations: and the new model (1)-(10).Next, numerical simulations are performed to support the theoretical stability results (successful and nonunion equilibria) and to numerically monitor the healing progression of a moderate fracture in normal conditions.Another set of numerical simulations is performed to analyze the effects of different debris densities on bone fracture healing.Finally, a set of numerical simulation results is presented to investigate the effects of different concentrations of anti-inflammatory cytokines and various cellular treatments on the fracture healing under numerous pathological conditions.All simulations are obtained by using the adaptive MATLAB solver ode23s and are initiated with densities of debris, macrophages, and MSCs set to D(0) = 5 × 10 7 , M 0 (0) = 4000, C m (0) = 1000, respectively, and the pro-inflammatory cytokines concentration set to c 1 (0) = 1.

6. 2 .
Different outcomes of the bone fracture healing process Next, a set of numerical simulations is presented to support the theoretical results.Accordingly to the qualitative analysis of the model there are three equilibria: E 0 , E 1 and E 2 where their stability conditions are determined by the tissue cells' proliferation and differentiation rates, k pm , k pb , d m and d b , respectively.The following parameter values are used: k pm = 0.5, d m = 1, k pb = 0.2202, and d b = 0.3, to demonstrate the stability of E 0 , since then k pm < d m and k pb < d b .The stability of E 1 is demonstrated using the following parameter values: d m = 1, k pm = 0.5, k pb = 0.2202, and d b = 0.15, since then k pm ≤ d m and k pb > d b .Finally, the following parameter values are used: k pm = 0.5 and d m = 0.1, to demonstrate the stability of E 2 , since then k pm > d m .Different time-periods are used in Figures 3-5 to better demonstrate the qualitative behavior of the system under different stability conditions.

Figure 3 .Figure 4 .
Figure 3. Cellular and molecular evolution of the resolution of the inflammation in normal conditions.

Figure 5 .
Figure 5. Cellular and molecular evolution of the repair process in a nonunion fracture healing.

Figure 6 .
Figure 6.Tissues evolution of a successful repair for different types of fractures.

Figure 13 .
Figure 13.Tissues evolution in a severe fracture without therapeutic innervation (solid line) and with M 0 (0) and C m (0) transplantation (dotted line).
Figures 13,14, and 15 show that the two cellular interventions increase both tissues productions.Furthermore, those interventions result in larger improvements in severe and senil osteoporotic fractures when compared to fractures in aging individuals.

J 2 (
d m − k pm ≥ 0 and k pb > d b and all the eigenvalues of J 1 (E 0 ) are non-positive values, then the eigenvalues of J(E 1 ) are negative except the eigenvalues associated with D and C m when k pm = d m , which are equal to zero.Therefore, E 1 is a locally stable node, since D ≤ 0 for all the variables of the system (1)-(10) and C m ≤ 0 when k pm = d m .

2 , H * 1 = a 12 a 12 +c * 2 , H * 2 = a 22 a 22 +c * 2 and J 11
is defined as in Theorem A1.Since all the eigenvalues of J 11 are negative (Theorem A1) and k pm > d m , and all equilibrium variables and parameter values

Table 1 .
Existence conditions for the equilibrium points and their biological meaning.

Table 2 .
Stability conditions for the equilibrium points.

Table 3 .
Parameter descriptions and units.