The Inﬂuence of Loading Program on the Stimulated Callus Mineralization

: The aim of this work is to describe the change of physical properties of the callus material with the use of the proposed mathematical model for callus remodeling. Callus tissue can be considered as a biomaterial where it’s properties change over time due to the stimulated healing process. The proposed model is based on the mechanical stimulus theory. It is used to estimate the stress-stimulated change in the callus, Young’s modulus, and the density in the case of a mandible fracture. Three healing loading programs are discussed and compared: optimal, intermittent, and intermittent with residual load. Here, the optimal loading program is understood as the in-time change of stimulating loads, which results in the shortest necessary healing time and, simultaneously, in the most uniform distribution of material density in the analyzed domain. The necessary healing time is a period after which the callus density (and hence the Young’s modulus) reaches the desired value. The results of the study suggest a signiﬁcant di ﬀ erence in the value of the callus maximal density between all three analyzed loading programs for a given healing time interval. The highest values of the density are obtained using the optimal loading program, however, all three programs provide satisfactory density distributions. The analytical results are compared with the callus density estimation based on the computer tomography (CT) medical data.


Introduction
The mandible (lower jaw) is the only load-bearing, moveable bone in the skull. Fractures in the mandible may happen as a result of infection, car accident, falling, sudden hit, or after a resection for pathology, e.g., cancer, treatment. As the only moveable bone of the cranium (skull), the mandible is subject to stresses and because of its shape, location, and function, the loads are unique to it [1]. This work focuses on the mandible fractures and healing process of the callus tissue after trauma. Such processes are known and thoroughly described for long bones, but in literature, there are no such works about mandible where complex states of stress and strain occur. The callus is a tissue formed between two bone parts after bone trauma. This process consists of three phases: reactive, reparative, and the remodeling phase during which different cells take part. In the first hours after the trauma, fibroblasts replicate and form a loose aggregate of cells, known as granulation tissue. Then, the fibroblasts within the granulation tissue develop into chondroblasts and this process culminates in a new mass of heterogeneous tissue, which is known as a fracture callus. The mineralized matrix is penetrated by osteoblasts, which form new lamellar bone in the form of the trabecular bone [2]. During the healing process, the properties of the callus tissue change make it mechanically stronger. The Young's (elasticity) modulus of the callus, which increases in this process of mineralization, is used as an indicator of material mechanical properties. The soft callus tissue transforms into bone tissue [2]. There is a need for the formulation of mechanical models for biomaterial transformation, To sum up, only a few models proposed in the literature take into account bone resorption due to overload [3,10,12] and there are no such studies in relation to the callus material, which in the process of mineralization transforms into a compact bone tissue. Due to these facts, the aim of this work is to provide a new approach based on the mechanical stimulus theory to describe the change of the callus density in the process of healing. With the use of a mathematical model for callus remodeling [21,22], the analyses of various loading programs are performed. In this work, we mainly focus on the elastic modulus, E, of the callus tissue, i.e., we aim at obtaining its sufficiently high value and homogeneous distribution.
For the initial model validation, the medical data of the fractured mandible is analyzed. Firstly, the callus region based on the computer tomography data using Mimics software is marked, then the density measurements are performed. The aim of this analysis is to determine whether it is possible to determine the modulus of elasticity of the callus and compare its distribution qualitatively with the numerical analysis results.

Geometry of the Callus Tissue
The proposed callus remodeling model [21,22] describes density change based on the mechanical stimulus theory using strain energy density (SED) in the stimulus definition. It should be clarified that the proposed model of tissue remodeling is used to analyze density changes only in the narrow section of the broken mandible in the process of stimulated healing and not in the mandible bone tissue. This section is a thin layer of the callus tissue of a shape similar to a rectangle and of 0.1-1.0 mm thickness (Figure 1), possibly stimulated during healing time. The section is subjected to loading as an effect of biting and muscular forces action. The result of load reduction to its centroid are mainly bending moment and shear force. In order to investigate the complex (not uniaxial) state of stresses and strains, a rectangular cross-section (b × h) loaded with bending moment and transverse force is considered. The initial value of the Young's modulus (E 0 ) is set to 2 MPa [23] and its initially uniform distribution is assumed all over the section. Due to non-uniform distribution of the stimulus, non-uniform healing is observed. The Young's modulus which indicates that a bone may be loaded physiologically is set to 200 MPa after Knets [24]. Various loading programs may be considered for in-time change of the bending moment and shear force. The goal is to propose such a program that results in the fastest mineralization during mandible healing. In the literature, a bone is considered to be an orthotropic or anisotropic material. However, in most finite element simulations, bone material (cortical tissue) is assumed to exhibit isotropic behavior [15,17,[25][26][27][28][29][30][31]. There are mainly three reasons for that. Firstly, the relations between Young's modulus and bone mass density are established experimentally mainly based on assumptions of material isotropy. Secondly, due to complex bone internal architecture, it is difficult to determine the orientation of the material properties in a bone. Thirdly, Milewski [32] conducted a finite element comparative analysis for the mandible with two different material properties (i.e., isotropic and orthotropic) and only a small difference was found between the results. From the biological point of view, in the process of mineralization, the callus tissue differentiates into mature bone tissue that exhibits the features of the cortical bone [2,33,34]. For these reasons, the isotropic material behavior for callus tissue is assumed in our analyses. It is worth mentioning that there are known various relations between tissue density and its Young's modulus [13,[35][36][37][38]. Some mathematical relations between the bone mass density and the mechanical properties (Young's modulus) of the bone were experimentally established by other authors [25,31,37]. In this work, the elasticity modulus E (MPa) of the callus tissue is related to tissue mass density ρ (g/cm 3 ): where C and γ are the constants to be experimentally determined [39][40][41].
In this study C = 3790 MPa·cm 9 /g 3 and γ = 3 are accepted following [13,14]. The callus tissue is modeled as a linear, elastic, and isotropic material for which constitutive relation is described with the use of two constants (the Young's modulus, E, and the Poisson ratio, ν).
According to Frost's theory [9], a bone can adapt its internal architecture of the trabeculae together with its mechanical properties with regard to various mechanical stimuli. The strain energy density (SED) is widely used by other authors to predict and evaluate the remodeling process [3,10,11,13,20]. Moreover, the SED is very convenient to use because of its scalar character. In the proposed model, the local (point-level) callus density change is a function of the chosen mechanical stimulus, following the equation of remodeling rate [21]: where ρ is the bone mass density (g/cm 3 ), ψ denotes SED/ρ (J/g), B and D are remodeling constants, K min , K max , and K over are the remodeling threshold values (J/g) and t (time unit) is the time parameter. K w is the maximal rate of density change (g/(cm 3 ·time unit)) calculated as the top of parabola (Equation (2) row (3)) for ψ = K over . The ordinary differential equation (Equation (2)) is solved numerically using the forward Euler method: Table 1 sums up the parameters used in the numerical analyses of bone remodeling in the mandible. It is worth pointing out that such parameters as K = K min +K max 2 , the effect of lazy zone, and the time constant regulating the remodeling rate are quite different in the literature.  [12] 0.004 (J/g) -1 ((g/cm 3 ) 2 (MPa·time unit) −1 ) Mellar et al. 2004 [13] 0.004 (J/g) 10% 1 ((g/cm 3 ) 2 (MPa·time unit) −1 ) Lian et al. 2010 [15] 0.016 (J/g) 20% 1 ((g/cm 3 ) 2 (MPa·time unit) −1 ) δ-percentage denoting the region of the lazy zone.
The following values are used in this work for the model parameters: stands for the end of interval with positive remodeling rate (ψ ∈ (K max , K d )), K min = 0.0036 (J/g), K max = 0.0044 (J/g), K over = 0.0358 (J/g), δ = 10%, K w = 0.0314 (g/(cm 3 ·time unit)). In our simulation, at each time step, the change in callus density is calculated using Equation (2) at every material point of the rectangular section. Then, the corresponding Young's modulus of the bone tissue is locally updated according to the Equation (1) and the analysis is carried out using the new, modified material properties. The iterative process continues until an acceptable distribution of the callus elasticity modulus is reached.

Analytical Formulation
As stated above, the rectangle is used to model the callus section of a broken mandible ( Figure 1). The rectangle is loaded with the bending moment (couple) M together with the shear (transversal) force T. This results in non-homogeneous stress distribution over the section surface. Namely, they vary along the section height (coordinate z ∈ [−h/2, h/2]). Normal stress is expressed as: where I E (mm 4 ·MPa) denotes the weighted (generalized) moment of inertia, and the shear stress is: where S E (mm 3 ·MPa) denotes the weighted static moment about a neutral axis of section cut-off lying on one side of level z. Hence, the stimulus becomes a function of z, as do the Young's modulus, E, and density, ρ. A complex stress state implies that the value of the strain energy density, SED, varies along the height of the section:

Optimal Loading Program
The previous studies [21,22], showed that for a constant-in-time load value, it is impossible to obtain a sufficiently high Young's modulus even after an infinitely long time of healing. The question arises of how to choose a loading program (understood as functions in time for M(t) and T(t)) so as to obtain the desired distribution of E(z). In order to compare various distributions of the callus Young's modulus over the healing domain (rectangle), as well as to assess sufficiency and homogeneity of a particular distribution, three parameters are defined and calculated: maximal value (E max ), mean value (E m ), and deviation parameter (s) The less is s, the less is difference E max − E m , the better (the more homogeneous) is the obtained distribution of E(z). For the most desired uniform distribution, E(z) = const z , the deviation parameter s = 0. The closer E max is to 200 (MPa), the more sufficient the distribution is. Hence, the optimal loading program is defined as either functions of M(t) and T(t), which minimize s ( s → min ) and maximize E max ( E max → max ) in a given time interval, or equivalently, such functions of M(t) and T(t), which minimize s ( s → min ) and minimize the time necessary to obtain a desired value of E max (e.g., E max = 200 (MPa)).

Results
The results are presented in Figures 2-4 for three analyzed loading programs. On the right side, each of ten colored lines represent the Young's modulus distribution along the height of the analyzed rectangular section ( Figure 1C) for the given time step (t = 1, 2, 3, . . . , 10). The first line on the left represents the starting point (t = 0) with the initial value of the callus elasticity modulus equal to 2 MPa. Moving to the right, along the "time increase" arrow, we observe an increase of the Young's modulus maximal value and the mean value in the callus. The three parameters of homogeneity are calculated for the final time step (t = 10 (time unit)).    Figure 2 presents the outcome of a progressive remodeling process for the optimal loading program. The values of the bending moment, M, and shear force, T, in time, are calculated using the control algorithm described in [21,22]. This control algorithm is used for the load parameters so as to improve homogeneity of the elasticity modulus with possibly maximal remodeling rate. At each time step, we require that the value of the stimulus parameter ψ is kept in a narrow vicinity of the optimal value, K over , for all section points, which denotes that the remodeling rate remains close to the maximal possible value, K w . Satisfying min z ψ(z) = (1 − α)K over and max z ψ(z) = (1 + α)K over (α determines the vicinity), we come to a set of equations to determine the optimal values of the bending moment, M, and the shear force, T. They are related to the distribution of E(z) obtained at the preceding time steps, starting from E(z) = const z max and z min denote coordinates of section points where the stimulus is maximal or minimal, respectively. With this optimal loading program it is possible to reach unrestricted value of the elasticity modulus with its distribution kept close to homogeneous. In Figure 2 It is unlikely that in practical treatment of a patient, the load parameters (M and T) can be changed in continuous and optimal manner. The loading program might be interrupted and for a designated time interval, the loads may diminish to zero. The aim of investigating this intermittent loading program is to assess its influence on the final distribution of the callus elasticity modulus. The optimal load values (M and T) are applied for a given time interval (∆t load = 1 (time unit)), hence the callus density (and elasticity modulus) increases due to Equation (2). Then the loads drop to zero for equal time intervals and the modulus decreases due to the resorption effect (ψ < K min ). The process is repeated that way for 10 time units. Figure 3 presents the outcome of the remodeling process for the intermittent loading program in which one may observe that, in general, a satisfactory homogeneous distribution of the Young's modulus in a section is achieved (s = 0.019). However, the maximal value of the elasticity modulus is almost five times lower than the value obtained for the optimal loading program after the same healing time (Table 2). This loading program can still provide a sufficient value for the elasticity modulus of the callus tissue but a longer time period is required ( Table 3). The homogeneity of the "long time" distribution is worse than the one obtained for the continuous optimal loading program. Also, remodeling effectiveness depends on the ratio between loading (∆t load ) and unloading (∆t unload ) time intervals, and the bigger the ratio is, the higher the value of E max can be remodeled in the given time. This loading program was analyzed as physiological-like when during the healing process, the callus tissue is repeatedly stimulated and not stimulated for given time intervals.  Figure 4 presents the outcome of a progressive remodeling process for the intermittent loading program with the 10% residual load. That means the load values are assumed to reduce down to 10% of ones acting before unloading, instead of total vanishing. This program was analyzed only as a possibility to happen and the 10%, 20%, or 30% values for the residual loads are chosen to test how the model of the callus remodeling would respond to this type of loading scheme. Comparison of the load paths in Figures 2-4 indicates significantly different values. It is explained that optimal load values, calculated with Equations (9) and (10), are not established once for good. In contrary, they strongly depend on healing process history, exhibited by the modulus distribution E(z), obtained after all the preceding time steps. If there is no unloading, as in the optimal program, the modulus increases continuously in time at all section points, and hence the required values of optimal M and T increase too. For the sub-optimal intermittent programs, the modulus increases during loading but decreases during each unloading period. In general, if a following loading is to start after an unloading time, the optimal values of loads are less because of smaller values of E(z).
One may observe for the intermittent programs that for a given time interval (Table 2) the higher is the residual load the better are values of E max , but the distribution homogeneity is worse. Similarly, Table 3 indicates that for the intermittent programs the higher is the residual load the shorter is healing time, necessary for satisfactory E max value. However, the distribution homogeneity is worse. The explanation is, that during total unloading the resorption occurs with the same rate (−BK min (g/(cm 3 ·time unit)) at each section point, since the stimulator ψ = 0 for each coordinate z. That means, the modulus is reduced for the same value everywhere and the distribution homogeneity is not changed after unloading. For partial unloading (intermittent with residuum), it happens that the modulus decreases in some section regions (where ψ < K min ), or remains unchanged in some regions (where K min < ψ < K max ), or increases elsewhere (where ψ > K max ). The last possibility is more likely to happen for higher residual loads. This results in the deterioration of the distribution homogeneity which is improved after the loading period.

Callus Density Measurements
The clinical studies of computer tomography data are strictly limited and mostly difficult to obtain. The patient's medical scans are performed during clinical appointments, which are rarely regular as mostly they take place with unevenly spaced time intervals. Therefore, it is difficult to collect the supportive clinical data for model validation, especially for estimating its time-scale parameters. The idea is to compare a series of scans and hence, knowing the true healing time, to scale the model and establish the true time unit. Thus, in this research, the data set includes only two CT scans of a fractured mandible. The first CT scan was obtained a day after bone trauma, the second one 18 days after. In this study, we focus on determining the callus Young's modulus based on the CT scans using Mimics software (V20.0, Materialise Co. Ltd., Leuven, Belgium). For this reason, we have excluded the first set of CT scans as there are no visible signs of the callus and the broken bone parts are evidently separated ( Figure 5). For the second set of CT scans, eight planes perpendicular to the forming callus tissue are chosen for analyses. They are shown in Figure 6 where the callus area for the density measurements is marked as a region between bone parts (red arrows) and shown as planes from the molar area to the area beside the fixing plate. Then, the measurements of radiological density units HU (Hounsfield scale) are performed around the formed callus tissue at every plane in the region between the bone parts. One may observe a sudden drop of the radiological density indicating the formed callus tissue, which is significantly lower than the one for the mandible bone parts. The tissue mass density is calculated by the software's built-in equations. Then, using Equation (1), Young's modulus of the callus tissue is calculated at every investigated point. The Young's modulus distribution of the callus tissue along profile lines is presented in Figure 7. One may observe a sudden drop of the Young's modulus along the profile lines indicating the formed callus tissue, where the modulus of elasticity is significantly lower than the one for the mandible bone parts. The minimum value for every profile line is chosen for further analysis. It is assumed that the minimum value indicates the region of forming of the callus tissue and represents its properties. The results of the Young's modulus calculations for every plane are shown in Figure 7. The highest values of the elasticity modulus are observed for the planes at the upper part of the mandible, near the molar. The callus tissue Young's modulus decreases in the direction from the molar to the fixing plate (Figure 5b blue part). The observed distribution of the modulus remains in good agreement with the calculated results in the qualitative sense. It indicates that the highest value is met at the top of the callus layer which is stimulated the most because of the highest value of normal stress.   It is not possible to compare directly or verify the results of the numerical or analytical analyses mainly due to the different loading program during patient's rehabilitation. The exact values of the loads and forces are quite impossible to establish post-healing. Therefore the goal of this part of the work is to present a way to mark the callus region based on the CT data using Mimics software. As presented, it is possible to determine the modulus of elasticity based on the radiological density and compare its distribution qualitatively with the numerical analyses results.
For more adequate results, it seems necessary to obtain medical data from a patient at subsequent time intervals during the healing process. Unfortunately, CT scans are rarely performed in a short time space because of the fear for the patient's health. Due to the X-rays used in that study, there is a risk of osteoradionecrosis (cell death due to radiation). The presented correlation between clinical data and numerical results are satisfactory at this moment of analyses but need further investigation with the full 3D model of the mandible.

Discussion
In our study the computational callus remodeling model is used to estimate how the loading program effects the callus density distribution during the stimulated healing process. Also, an attempt to establish the optimal loading program was made, leading to the most homogeneous distribution of material properties with nominal values kept on a desired level. Two sub-optimal approaches are studied in order to get closer to possible clinical treatment. Based on the numerical results, the callus remodeling response can be evaluated. The results of this study suggest a significant difference in the callus maximal density value between all three analyzed loading programs performed in the same time interval. The highest values of the density is obtained using the optimally-controlled model. All three programs provide satisfactory density distribution, although the intermittent programs require more time to achieve the designated 200 MPa Young's modulus value.
There are three assumptions that have to be considered as limitations to accurately predict the outcomes of the callus distribution. First, the callus remodeling algorithm used in this work does not include the individual differences of the tissue density due to sex, age, diseases, osteoporosis, pharmacology treatment, and other patient-specific factors. Therefore, the influence of these factors on the callus remodeling process should be considered in the future study. That could mainly be achieved by variation of the model constant parameters. Due to that fact, the results presented and discussed here must be taken qualitatively and quantitative data will be different for each particular patient case. Second, the starting point of the iterative analysis is a homogeneous distribution of tissue density with Young's modulus equal to 2 MPa. In fact, the callus density in the first day of the healing process was not fully investigated or thoroughly described [23,[42][43][44]. Third, isotropic material properties were considered instead of anisotropic, which are closer to the true material. However, as mentioned before, in the process of mineralization, the callus tissue acquires the features of a cortical bone, which in the literature is mainly modeled as an isotropic material. Moreover, the very small thickness of the callus layer justifies this simplification.
One may observe from the above results that the proposed loading programs have a big influence on the bone density distribution. The predicted outcomes are not far different to those experienced in reality. It is known that too big a load can cause substantial damage or irreversible fracture in the tissues which are not strong enough [45,46]. Therefore, to strengthen the callus tissue after trauma, the load must be applied at an adequate rate and of an adequate value dictated by the bone density at the particular point of healing time. Qualitative and quantitative comparison between results shows a sufficient degree of resemblance between the computational results and the clinical CT data. Also, numerical studies prove that efficient healing of the tissue may be obtained in the realistic load-controlled treatment (intermittent program).
Simulations of callus remodeling, validated and compared with clinical data, are proven as a helpful way to elucidate the mechanisms behind bone healing. The effect of callus remodeling and the influence of the loading program on the healing process are critical for the development of planning the rehabilitation process. It is expected that this analyses of callus remodeling can contribute to the development of new loading programs used during patient rehabilitation. In the future, a detailed 3D analysis will be performed so that more a complex model of the callus and the mandible can be considered and more thorough solutions can be reached. The outcome of this research demonstrates the importance of the computational biomechanics due to the limitations in obtaining clinical data and experimental results. It substantiates the understanding of callus remodeling and provide new insights into bone healing. Additionally, it extends the potential of computer simulations in the field of bone biomechanics and hopefully, eventually reduce animal studies and clinical trials. Thoroughly examined numerical model of the callus stimulated transformation can be implemented into a finite element analysis of realistic (in shape and loading scheme) bones with thin layers of callus after trauma.