An Upstream Finite Volume Scheme for a Bone Healing Model

This paper is devoted to the introduction of a numerical scheme for a bone healing model and simulation of skull fractures. The mathematical model describes the evolution of mesenchymal stem cells, osteoblasts, bone matrix and osteogenic growth factor. We propose a numerical scheme based on an implicit finite volume method constructed on an orthogonal mesh. The efficiency and robustness of the scheme are shown in simulating a skull fracture in rats. Keywords-bone healing; finite volume method; fracture fimulation


I. INTRODUCTION
Bone is a tissue with a remarkable ability to regenerate itself.But for large gap sizes bone fails to heal itself in a clinically reasonable period of time, this problem is called non-union or delayed union.Approximately 5 − 10% of the 5.6 million fractures occurring annually in the United States develop into non-unions or delayed unions [1].There exists different methods to help bone healing, such as autograft or synthesis materials.In Europe, 1.5 million bone grafts are carried out to treat these fractures [2].An exterior help often allows the complete bone healing but there is still clinical cases which don't heal.For these complex cases, biology and medicine researchers would create ex vivo tissues that would then be reimplanted.The most common process therefore consists in the growth, in a bioreactor, of mesenchymal stem cells seeded on a synthesis material [3].Currently, this process yields good results only for two-dimensional bone growth in Petri dishes.Understanding the bone healing is fundamental to create tridimensional ex vivo bones and it is an important field of tissue engineering researches.There exists several mathematical models that simulate the bone healing [4], [5], [6].In this paper, we propose one such mathematical model (based on the one developed in [7]), describe a numerical scheme to approximate its solution and show some numerical illustrations that simulate a healing process in skull bones.Stability and convergence results for the numerical scheme are stated.

II. A MODEL FOR POPULATION DYNAMICS
The bone healing is a complex phenomenon involving a cascade of cellular and tissue events [8], [9].In this paper, we study a simplistic model but well-adapted to the bone growth in bioreactor describing the rates of change, with respect to time and space, of the concentrations in the mesenchymal stem cells s, the osteoblasts b, the extracellular bone matrix m and the osteogenic growth factor g.This dimensionless model is a simplified version of the model proposed by Bailón-Plaza and Van Der Meulen in [7] that reads for t > 0 and x ∈ Ω, where Ω is an open bounded polyhedral and connected subset of R d (d = 2, 3).The reaction terms f i (i = 1, . . .4) describe the exchange, production, decay, etc of the four populations of interest.The evolution of the stem cells is described by a diffusion term, a haptotaxis term (directional motility of the cells up a gradient of cellular adhesion sites, here the bone matrix) and a reaction term f 1 describing the mitosis and the differentiation into osteoblasts of the stem cells The coefficient of diffusion Λ and the velocity V of haptotaxis of the stem cells are non-linear functions.The function χ is given by It allows to avoid an infinite accumulation of the stem cells due to the haptotaxis term.The osteoblasts are only regulated by a reaction term f 2 because they are cling on the bone matrix.This term describes the mitosis, the decay of the osteoblasts and the osteoblastic differentiation The term f 3 describes the synthesis and the degradation of the bone matrix by the osteoblasts synthesis and degradation .
The rate of change of the growth factor is described by a diffusion term and a reaction term f 4 describing the production by the osteoblasts and the decay of the growth factor The parameters α i , β i , γ i , η i , δ i (i = 1, 2), ρ, λ and κ are real positive numbers.These four equations are completed by the homogeneous Neumann boundary conditions on s and g: for t > 0 and x ∈ ∂Ω, where n is the outward unit normal of ∂Ω; and by the data of initial conditions on s, b, m and g:

III. THE NUMERICAL SCHEME
In order to simulate bone healing, we need a numerical scheme that approximate the solutions to the equations ( 1) to (4).The approximate solution must remain bounded in the region of physically bounded solutions s, b, m and g, defined by Hence the approximate solution converges towards a weak solution of the equations.
The space discretization is based on an admissible mesh as defined in [10].It is a finite family T of polygonal open convex subsets K of Ω, called the control volumes such that Ω = ∪ K∈T K, together with a finite family E of disjoint subsets of Ω consisting in non-empty open convex subsets σ of affine hyperplanes of R d , called the edges, and a family P = {x K , K ∈ T } of points in Ω, called the centers verifying the following properties.
• Each σ ∈ E is contained in ∂K for some K ∈ T and for any K ∈ T , there exists a subset E K of E such that ∂K = ∪ σ∈EK σ.For any edge σ ∈ E, either σ ⊂ ∂Ω or σ = K ∩ L for some K = L in T .In the latter case, we denote σ = σ KL , called the interfaces.We denote by E ⊂ E the subset of all the interfaces and, for any K ∈ T , by Additionally, for any σ KL ∈ E , we denote by n KL and d KL , respectively, the unit vector normal to σ KL outward of K and the distance |x K −x L |.The measure of K ∈ T is denoted by |K| and the (d − 1)-dimensional measure of σ ∈ E is denoted by |σ|.
The time discretization is the sequence of discrete times t n = n∆t for n ∈ N where ∆t > 0 is a given time-step.
The numerical scheme is obtained by using the finite volume method: equations (1) to (4) of the model are integrated on each control volume K and interval of time (t n , t n+1 ).By using the divergence theorem, we obtain the following scheme , where the unknowns |K| K m(t n+1 , x)dx and 1 |K| K g(t n+1 , x)dx.An implicit time stepping strategy is used for all the terms.The approximations of Λ and V at an interface are calculated with an arithmetic mean between the two neighboring control volumes.The haptotaxis term is approximated by a flux F .Although an upstream flux would be welladapted, it does not ensure a maximum principle on the discrete solution.Consequently, the flux is defined such as follow: For this numerical scheme, we have proved the following results.

Theorem 3.1 (Existence of an admissible solution):
If the initial data is physically admissible, specifically if (s 0 K , b 0 K , m 0 K , g 0 K ) belongs to A for all control volumes K in T , then the discrete system of equations has a solution for all n ∈ N, which is physically admissible:

(Convergence to a weak solution):
If the initial data is physically admissible and the discrete equivalent of the H 1 semi-norm [10] of b 0 and m 0 are bounded, then there exists a subsequence (u h ) of discrete solutions that converges to a function u = (s, b, m, g) almost everywhere in [0, T ] × Ω (for any T > 0).This function u is an admissible (u(t) ∈ A for all t > 0) weak solution of the model.

IV. NUMERICAL SIMULATION: A SKULL FRACTURE
In order to validate the interest of the model, we simulate the healing of a skull fracture in rats.It is wellsuited to our model because there is no cartilage in this kind of fractures although it is essential in many cases.The simulation corresponds to an experience presented in [11].In this article, defects were created using a 2.3 mm outer diameter trephine in the parietal bones of 6 adult sprague-Dawley rats and the rats were allowed to heal for 42 days.To compute numerical solutions, we implement a simplified semi-implicit time-stepping technique with the Newton's method coupled with a biconjugate gradient method to solve the nonlinear system arising from the discretization.A difficulty in the implementation is to construct admissible meshes satisfying the orthogonality condition.Structured rectangular meshes are admissible meshes but they can not be used for complex geometries, like circular fractures for skull bones.We choose to use Voronoi diagrams, that verify the property of orthogonality between the interfaces and their respective centers.But most mesh generators give Voronoi diagrams with very small interfaces.To avoid this, we give a set of points well distributed in the domain, next we construct the Voronoi diagram associated to these points (figure 1), which represents the centers of the control volumes.Consequently, the number of interfaces is different for each control volumes.Fig. 2: Full geometry of the fracture.Fig. 3: One quarter of the domain where the black area corresponds to the bone (m 0 (x, y) = 0.1 g.ml -1 ) and cellular cluster (s 0 (x, y) = 10 6 cells.ml - and g 0 (x, y) = 2 × 10 3 ng.ml - ).Elsewhere, there is nothing initially.
The geometry of the circular skull fracture is reported on figure 2. The symmetries about fracture line and bone axis implies that only one-quarter of the domain needs to be considered (figure 3).Initially, the domain contains only the bone and two cell clusters along the broken bone made up of stem cells and growth factor (figure 3).The mesh used is a Voronoi mesh made up of 5884 control volumes (figure 1) and the time-step is fixed at ∆t = 14 minutes and 24 seconds.After 3 days, we  observe the formation of osteoblasts (figure 4c) where the stem cells are initially concentrated, these osteoblasts synthesized a new bone matrix (figure 4a).The stem cells moved towards the center of the fracture (figure 4b).The osteoblasts trapped in the new bone become osteocytes.
This model successfully simulates the evolution of the mineralization front (figure 5).At t = 42 days, since all stem cells disappear, we can consider that the healing is terminated.The defect heals approximatively 42% (close to the results of the article [11]), it is a nonunion fracture because the initial defect is too large.

V. CONCLUSION
In this article, we proposes a model and a finite volume numerical method to simulate bone regeneration that is well-suited to the growth in bioreactor.The approximate solutions remain in a physically admissible region, and the convergence to a weak solution of the equation is guaranteed.The simulation of a skull fracture allows to validate this model for bone healing without cartilage.Now, we plan to develop another model in order to simulate the growth in a bioreactor.It ought to include the flow environment and its interaction with the previous model.Modeling this interaction is an important challenge in order to understand three-dimensional ex vivo tissue culture.

Fig. 1 :
Fig. 1: Voronoi diagram used for the simulation of a skull fracture (5884 control volumes)

1 Fig. 4 :
Fig. 4: Bone matrix density, concentrations of stem cells, osteoblasts and the growth factor at t = 3 days.

Fig. 5 :
Fig. 5: Bone matrix density along the line y = x as a function of distance to the origin (x, y) = (0, 0).