Modeling of an initial stage of bone fracture healing

In case of the secondary bone fracture healing, four characteristic steps are often distinguished. The first stage, hematoma and clot formation, which is an object of our study, is important because it prepares the environment for the following stages. In this work, a new mathematical model describing basic effects present short after the injury is proposed. The main idea is based on the assumption that blood leaking from the ruptured blood vessels propagates into a poroelastic saturated tissue close to the fracture and mixes with the interstitial liquid present in pores. After certain time period from the first contact with surrounding tissue, the solidification of blood in the fluid mixture starts. This results in clot formation. By assuming the time necessary to initiate solidification and critical saturation of blood in the mixture, the shape and the structure of blood clot could be determined. In numerical example, proposed mathematical formulas were used to study the size of the gap between fractured parts and its effect in blood clot formation.

RTG images of fractured long bone and formation of callus; from the left: 1 week after injury, 4 weeks after injury, 7 weeks after injury of cellular events. Therefore, without correct clot formation, growth of callus, working as a bridge between broken parts of bone, is difficult or sometimes impossible. It later on may cause nonunion or creation of false joint.
Two major schemes of fracture healing are commonly distinguished, namely primary and secondary healing, see e.g., [16,43]. In the primary (or direct) bone healing, the cortical remodeling is characteristic and predominating phenomenon. The main cell activities result in creation of so called cutting cones. It is a characteristic structure arranged initially by osteoclasts tunneling across the fracture and later on activities of osteoblasts. Haversian structures, the basic unit of bone, are formed by the activity of cutting cones. This scenario occurs in cases when a fracture is well stabilized and compressed when fractured fragments are in direct contact. In most cases, however, such situation does not happen, and then, the secondary bone healing is present that is an object of our study in this work. The great majority of fractures undergo secondary bone healing, which involves both endochondral and intramembranous ossification and requires some relative motions of fractured fragments. These motions represent an important aspect in healing process. They cannot be too small since they act as mechanical stimulus, and, on the other hand, too large motion results in damage of newly created structures.
Four conventional phases can be distinguished in secondary bone healing. (1) Hematoma and next clot formation. At this stage, blood from ruptured vessels propagates in porous matrix around the fracture, solidifies in contact with other tissues and forms a clot. (2) Inflammatory reaction and angiogenesis initiation. This phase is characteristic for removal of dead fragments of tissue and useless mineral elements by activities of macrophages and other specialized cells, and initiation of tissue differentiation. (3) Repair phase. The continuation of tissue differentiation results in callus formation, which functions as a bridge connecting fractured bone fragments. Some other effects including angiogenesis, cartilage formation and its calcification, and bone tissue formation are crucial in enabling mechanical stabilization of broken parts. Bone tissue, called woven bone, created at this stage is not fully functional and will be later replaced by high-quality lamellar bone. (4) Final remodeling phase. During this phase, callus is gradually being removed. New direct connection between broken bones takes over its mechanical functions. Finally, the woven bone is being replaced in remodeling process by fully functional, well-organized, strong, lamellar bone tissue [33]. In Fig. 1, RTG images of fractured bone at different healing stages are presented. Looking from left to right, the first picture illustrates the second stage when the inflammation process and initial formation of tissue are predominant; the second picture illustrates the next stage, when the creation of callus starts; the characteristic situation for the beginning of bone remodeling process is displayed in the last picture. Figure 2 illustrates the distribution and density of tissue at three cross sections of bone-far from the fracture, close to the fracture and at the site. These images were taken at the beginning of the last stage of healing. The difference of densities of old compact bone and newly created bone can be easily observed. This new tissue is still not fully mechanically functional and will be removed during remodeling process. The division of secondary fracture healing into four discrete phases is artificial. In reality, these phases overlap with each other. But it does make it easier to follow the processes involved.
The first phase, which is an object of this study, even though it is short compared to other phases, plays an important role in the whole process because it prepares an environment to the following three steps. First of all, the bleeding from broken blood vessels in periosteum and fractured bone leads to hematoma formation which serves as a source of cells and platelets involved in initiation of the inflammatory phase. Later on, solidification of blood contacting with other tissues results in creation of fibrous structure (clot). Clot among the others provides early mechanical stabilization of broken parts and plays a role as reservoir of different cells. These cells are active during inflammation and angiogenesis phase or are involved in later tissue regeneration by releasing signaling molecules and cytokines [15,17,38]. Clot domain defines a region of future activities present during next healing stages. Thus, the volume of hematoma and clot, their structure and correct formation are important issues affecting further healing steps.
It is well known and clinically observable that the distance between fractured parts of a broken bone and mechanical loading after injury are important structural and mechanical factors, which, in addition to other effects of biological nature, play a crucial role in healing and tissue regeneration processes [8,22]. The mechanical loading of a certain level is necessary to enable small micro-motions between fractured parts as to stimulate cells involved in tissue regeneration. On the other hand, newly formed structure connecting broken parts may be damaged by excessive loads. Such situation may lead to serious complications, for example, nonunion. Another important factor is the size between broken parts. Small and stabilized gap leads to direct healing, while too large gap size has negative effect on fracture healing which may lead to formation of false joint.
A considerable effort is also directed to build mathematical models and associated computational simulations of bone fracture healing and involved phenomena. Models to predict temporal and spatial evolution of ossification [3], and spatial pattern distribution of tissue [21,42] were built in order to unravel the complex relation between mechanical environment and tissue formation. Later on, a more advanced model which combines tissue differentiation rules with callus growth or resorption was developed by Lacroix [26] and Gómez [20,22]. Mathematical models are also capable to describe the cell behavior, such as migration and proliferation, and angiogenesis [4,5,8].
Besides development of mechano-biological rules and models, a great effort is made also in mathematical aspect. Numbers of advanced mathematical methods are involved to describe this healing procedure and associated effects. In [28], a general variational approach was proposed to derive bone remodeling rules in contrast to other works where such rules are postulated. Similar approach can be used to model the last stage of bone fracture healing. In [29], a model of interaction between bone and bone substitute material during regeneration process was proposed and extended for the second-gradient effects in [34]. A model of growth, mass transfer and remodeling in fiber-reinforced, multi-constituent materials was developed in [24]. These works together with [39] can serve as the basis for development of the complex model of bone fracture healing including various important effects.
Unlike all other models which focus on later healing stages, the aim of this work is to address mathematical description of the initial healing stage in which hematoma and clot are formed by the effect of bleeding into poroelastic tissue and later on blood solidification. The initial stage of bone fracture healing-hematoma and clot formation-is short (few hours) but important as it prepares the environment for following healing processes. We introduce a two-fluid mixture in porous matrix model to describe the bleeding into soft tissue. Soft tissue around the fracture is assumed to be a homogeneous poroelastic material, which is saturated by interstitial fluid. At this very initial period, a hematoma comes into being. After a certain time period from the first moment of blood contact with other tissues, the mixture of blood and other liquids starts to solidify and forms blood clot. In addition, a constitute equation, which describing matrix deformation, relates the increment in fluid content to volumetric strain, and incremental pore pressure is added. Numerical simulation of this model is performed to examine the effect of gap size in blood clot formation and its final shape in long bone fracture. This is a complex process where different effects are involved. These effects are important not only for healing but also in bone formation, regeneration, remodeling and some pathological changes in bones. Phenomena as for example tissue differentiation from mesenchymal stem cells and revascularization, which often play key roles in final effects of considered processes, are quite dependent on mechanical factors, see e.g., [6][7][8]23,25,26,37,40]. There are already many experimental and clinical investigations reported, see e.g., [2,6,8].
The further aim is to include the result of this work into the description of complete healing procedure composed of all four steps. That's because, as previously mentioned, the later healing steps undergo in the domain set by the shape and characteristics of the clot.

Methods
To build model of considered process and collect all of the required mathematical formulas, let us discuss the basic assumptions of our formulation.
Consider a long bone surrounded by an elastic porous tissue. Assume that bone is rigid in comparison with the surrounding poroelastic tissue. This poroelastic tissue is initially saturated with interstitial fluid. After injury, a fracture appears and the continuity of bone considered in our later numerical example as a simple tube is broken. The size of a gap between fractured parts of a bone plays an important role in fracture healing. To apply the condition of secondary healing, we assume that the broken parts of bone after injury do not contact with each other. In numerical simulations, the examination of the effect of gap size in healing process was studied and the results are shortly discussed in the next section.
Bone fracture and dislocation of its broken fragments is accompanied by a damage of blood vessels. Because of the pressure difference between blood in vessels and interstitial fluids in poroelastic tissue nearby bone fracture, blood leaks from the broken vessels and starts to penetrate into surrounding tissue. Blood mixes with interstitial fluid present in poroelastic surrounding tissue and propagates in pores. Concentration of blood in this fluid mixture varies in time and distance from the fracture. In addition to blood flow, we can also include in this theoretical consideration capillary diffusion effect. However, its contribution to the complete process should be examined experimentally in future investigations.
The fundamental effect that should be considered here is the solidification of blood. Blood outflow from vessels results in creation of hematoma close to the fractured site. After a certain time period from the first contact of blood with surrounding tissues, the solidification of blood starts and a clot is formed. This phenomenon plays an important role in organism and results in stopping blood loss. The solidification and formation of solid clot is possible in the domains where required concentration of blood in a mixture with interstitial fluid is assured. Thus, by assuming the time required to initiate solidification process and the value of minimal necessary concentration of blood, the domain occupied by a newly formed clot can be estimated.
In order to describe this mechanism, the flow of a mixture composed of two fluids, namely blood leaking from ruptured vessels and interstitial fluid present in poroelastic and deformable neighbor tissue should be considered. Making a more precise model, a diffusion of blood in porous media is also included. Additionally, blood solidification defined by two parameters-concentration of blood in flowing mixture and time period from the initial moment of blood leakage-is considered to determine the domain and structure of new clot. The structure, or more precisely-porosity of clot, can later be determined by a known value from the analysis of blood concentration in considered fluids mixture.
In Fig. 3, a schematic picture of the considered fractured region of a bone and surrounding soft tissue is presented. Boundaries with inflow condition locate at the two ends of fractured gap, representing ruptured blood vessels from periosteum and cortical bone. Boundaries with outflow condition located between inflow boundaries and at two vertical ends of considered soft tissue far enough form the fractured site to enable zero pressure condition.
Let us discuss in more details fundamental equations necessary to study described here effects. Consider a poroelastic material. The Navier's equation has a form (see, e.g., also [18] where it was used for modeling the confined compression of a cartilage sample in small displacement regime) where σ represents the total Cauchy stress tensor and is defined as and α B denotes Biot-Willis coefficient, p is the fluid pore pressure and C represents elasticity matrix defined for "drained" material. This is a standard formula for poroelastic material where a contribution of pore pressure is added in constitutive equation. Strain field is defined as where u denotes displacement field. Careful examination of Eq. 2 and splitting the stress into volumetric and deviatoric parts leads to the conclusion that the deviatoric part is pore pressure independent. Therefore, for the volumetric part, relation between total mean pressure p m and volumetric strain ε vol can be expressed as where and K d represents bulk modulus of the drained porous material. In the following paragraph, the equations for two-phase Darcy flow are collected. As we discussed earlier in this section, the flow of a mixture composed of fluids "1" and "2," namely blood and interstitial fluid, in a porous media is considered. Denote Darcy velocity by u, pressure of a fluid in pores by p, dynamic viscosity of mixture by μ and permeability of porous matrix by κ. According to Darcy's law, the velocity field is dependent on pressure gradient [1], Next important equation represents continuity condition. However, we deal with a mixture of two fluids. Therefore, some additional notation and formulas are necessary. Denote by κ 1 and κ 2 relative permeabilities for fluids "1" and "2," respectively, by μ 1 and μ 2 their dynamic viscosities and by s 1 and s 2 saturation (volume fractions) of both components. They satisfy the following condition, With this notation density of considered mixture can be expressed as, where ρ 1 and ρ 2 represent densities of component fluids, and dynamic viscosity of mixture follows from the formula The mass conservation condition can be written in the following form, where p denotes porosity of porous tissue and Q represents a source term. To obtain more convenient equation expressed in terms of pressure some, simple manipulations are necessary. Consider the first term in Eq. (10) and assume that both density and porosity are dependent of pressure. Then, simple differentiation leads to Applying in Eq. (11) definition of fluid compressibility and notation for storage S we obtain the storage model in a form Substitution of this formula and Darcy low (Eq. 6) leads to the equation expressed in terms of pressure, where the source term can be replaced by the rate of expansion of the pore space, In this theoretical consideration, the diffusion effect is also added. Let us denote capillary diffusion coefficient by D c and fluid "1" content by c 1 defined as With introduced here notation the transport equation can be written now as follows, The equations listed here describing penetration of blood after bone fracture in poroelastic tissue saturated initially with interstitial fluids. They have been used in simple numerical calculations to study the process of clot formation and gap size effect.

Results and discussion
The mathematical equations collected in the previous sections were used in a numerical simulation of a specific example. The goal of this calculation was to test proposed formulation and to examine the effect of fracture gap size in blood clot formation. In this two-dimensional example, a gap between fractured parts of bone and a poroelastic tissue close to the injury are considered, see Fig. 3. At the edges of broken bone, blood leaks from ruptured blood vessels present in periosteum and propagates in the domain occupied by poroelastic tissue saturated with interstitial fluid. The bone fragments attached to poroelastic tissue are considered rigid in comparison with this tissue. The composition of a mixture composed of blood and interstitial fluid and its distribution determine a domain where the clot is formed after blood solidification. Following are details of the examined case.  Figure 4 shows the evolution of blood propagation in soft tissue, and the dark red color represents its biggest saturation in pores. Assume initiation of solidification starts at 1,000 s after blood flowed out from ruptured vessels; the solidification is so quick that could be neglected comparing to whole procedure. In that way, we draw a clot shape base on saturation map after 1,000 s. When blood proportion in the fluid mixture is bigger than 0.7, clot may form. As in Fig. 5, the black lines show the borders of clot. Figure (a) shows small gap size case (gap = 0.5 cm), and clots from two fractured ends are fused well; (b) shows middle gap size case (gap = 1 cm), and clots from two fractured ends are barely contact with each other; (c) shows big gap size case (gap = 1.5 cm), and clots from two fractured ends are separated. As previously mentioned, clot is a source of different cells and signaling molecules which has the capacity to initiate the cascades of cellular events important for fracture healing. Therefore, domain near fracture without clot formed is hard for callus to grow and later on may cause nonunion. In addition to clot domain determination, the saturation of blood determines density and strength of a clot. This information is also important in considerations focused on the following stages of bone healing.

Conclusion
This simple model enables determination of shape, porosity and mechanical characteristics of blood clot and then initial shape of callus. It follows from our preliminary calculations that clot formation strongly depends on the gap size. This may explain the clinically observed effect that large gap size often results in unsuccessful formation of callus and failure of fracture healing.
Proposed model needs experimental validation and identification of the parameters introduced in mathematical description. In presented example, the porosity of surrounding tissue does not change in space and permeability decreases with distance-from maximal value at the interface with bone as to zero on the skin side. For better description of real situation, these assumptions can be improved by involving different important effects, as for example, randomly distributed cracks in soft tissue and pressure from the elasticity of skin. The biological effect associated with swelling also should be taken into consideration.
The results of numerical simulations are similar to clinical observations. Therefore, the calculated shape, structure and mechanical characteristics of clot can be used in analysis of next steps, in particular tissue differentiation which is strongly dependent on mechanical effects. It can be expected that among the others, the bifurcation of the solution might be one of the possible phenomenon determining the final effect. This is planned to apply the methods developed in [30][31][32] in the future study.