Original article Biomath 1 (2012), 1209254, 1–5 B f Volume ░, Number ░, 20░░ BIOMATH ISSN 1314-684X Editor–in–Chief: Roumen Anguelov B f BIOMATH h t t p : / / w w w . b i o m a t h f o r u m . o r g / b i o m a t h / i n d e x . p h p / b i o m a t h / Biomath Forum An Upstream Finite Volume Scheme for a Bone Healing Model Yves Coudière∗, Mazen Saad† and Alexandre Uzureau ∗ ∗ LMJL UMR6629 CNRS-UN Université de Nantes, Nantes, France Emails: yves.coudiere@univ-nantes.fr, alexandre.uzureau@univ-nantes.fr † LMJL UMR6629 CNRS-UN École Centrale de Nantes, Nantes, France Email: mazen.saad@ec-nantes.fr Received: 12 July 2012, accepted: 25 September 2012, published: 17 October 2012 Abstract—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 nu- merical 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. Under- standing the bone healing is fundamental to create tri- dimensional 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 concen- trations 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 Citation: Y. Coudière, M. Saad, A. Uzureau, An Upstream Finite Volume Scheme for a Bone Healing Model, Biomath 1 (2012), 1209254, http://dx.doi.org/10.11145/j.biomath.2012.09.254 Page 1 of 5 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2012.09.254 Y. Coudière et al., An Upstream Finite Volume Scheme for a Bone Healing Model Der Meulen in [7] that reads ∂ts = f1(s, m, g) + div(Λ(m)∇s︸ ︷︷ ︸ diffusion − V (m)χ(s)∇m︸ ︷︷ ︸ haptotaxis ), (1) ∂tb = f2(s, b, m, g), (2) ∂tm = f3(b, m), (3) ∂tg = f4(b, g) + div (Λg∇g)︸ ︷︷ ︸ diffusion , (4) for t > 0 and x ∈ Ω, where Ω is an open bounded polyhedral and connected subset of Rd (d = 2, 3). The reaction terms fi (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 f1 describing the mitosis and the differentiation into osteoblasts of the stem cells f1(s, m, g) = α1 β21 + m 2 ms (1 − s)︸ ︷︷ ︸ mitosis − γ1 η1 + g gs︸ ︷︷ ︸ differentiation . The coefficient of diffusion Λ and the velocity V of haptotaxis of the stem cells are non-linear functions. The function χ is given by χ(s) = s(1 − s). 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 f2 because they are cling on the bone matrix. This term describes the mitosis, the decay of the osteoblasts and the osteoblastic differentia- tion f2(s, b, m, g) = α2 β22 + m 2 mb (1 − b)︸ ︷︷ ︸ mitosis + ρ γ1 η1 + g gs︸ ︷︷ ︸ differentiation − δ1b︸︷︷︸ removal . The term f3 describes the synthesis and the degradation of the bone matrix by the osteoblasts f3(b, m) = λ (1 − κm) b︸ ︷︷ ︸ synthesis and degradation . The rate of change of the growth factor is described by a diffusion term and a reaction term f4 describing the production by the osteoblasts and the decay of the growth factor f4(b, g) = γ2 (η2 + g) 2 gb︸ ︷︷ ︸ production − δ2g︸︷︷︸ decay . 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 con- ditions on s and g: (Λ(m)∇s − V (m)χ(s)∇m) · n = 0, (Λg∇g) · n = 0 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: s(0, x) = s0(x), b(0, x) = b0(x), m(0, x) = m0(x), g(0, x) = g0(x), for x ∈ Ω. III. THE NUMERICAL SCHEME In order to simulate bone healing, we need a nu- merical 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 A = [0, 1]×[0, ρ γ1 δ1 ]× [0, 1 κ ] × [0, ργ1γ2 4δ1η2δ2 ] ⊂ R4 for ρ γ1 δ1 ≥ 1. 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 Rd, called the edges, and a family P = {xK , 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 EK of E such that ∂K = ∪σ∈EK σ̄. For any edge σ ∈ E, either σ ⊂ ∂Ω or σ = K̄ ∩ L̄ for some K 6= 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 N (K) = {L ∈ T , σKL ∈ EK ∩ E?} ⊂ T the neighbors of K. • For any K ∈ T , the point xK belongs to K. For any σKL ∈ E, the line (xK , xL) is orthogonal to σKL. Biomath 1 (2012), 1209254, http://dx.doi.org/10.11145/j.biomath.2012.09.254 Page 2 of 5 http://dx.doi.org/10.11145/j.biomath.2012.09.254 Y. Coudière et al., An Upstream Finite Volume Scheme for a Bone Healing Model Additionally, for any σKL ∈ E?, we denote by nKL and dKL, respectively, the unit vector normal to σKL outward of K and the distance |xK −xL|. 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 tn = 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 (tn, tn+1). By using the divergence theorem, we obtain the following scheme |K|(sn+1K − s n K ) − ∆t ∑ L∈N (K) Λn+1KL sn+1L − s n K dKL |σKL| + ∆t ∑ L∈N (K) F (sn+1K , s n+1 L , V n+1 KL mn+1L − m n+1 K dKL )|σKL| = ∆t|K|f1(sn+1K , m n+1 K , g n+1 K ) |K|(bn+1K − b n K ) = |K|∆tf2(s n+1 K , b n+1 K , m n+1 K , g n+1 K ) |K|(mn+1K − m n K ) = |K|∆tf3(b n+1 K , m n+1 K ), |K|(gn+1K − g n K ) − ∆t ∑ L∈N (K) Λg gn+1L − g n+1 K dKL |σKL| = |K|∆tf4(bn+1K , g n+1 K ), where the unknowns sn+1K , b n+1 K , m n+1 K and g n+1 K approximate 1|K| ∫ K s(tn+1, x)dx, 1|K| ∫ K b(tn+1, x)dx, 1 |K| ∫ K m(tn+1, x)dx and 1|K| ∫ K g(tn+1, x)dx. An im- plicit 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 well- adapted, it does not ensure a maximum principle on the discrete solution. Consequently, the flux is defined such as follow: F (a, b, c) = c+ (χ↑(a) + χ↓(b)) − c− (χ↑(b) + χ↓(a)) where c+ = max(c, 0), c− = max(−c, 0), χ↑(a) =∫ a 0 χ′(s)+ds and χ↓(a) = − ∫ a 0 χ′(s)−ds. This flux verify the two classical properties of conservativity and consistency and an additional property of monotony: for any (a, c) ∈ R2, the mapping b ∈ R 7→ F (a, b, c) is non-increasing. This ensures the maximum principle. For this numerical scheme, we have proved the fol- lowing results. Theorem 3.1 (Existence of an admissible solution): If the initial data is physically admissible, specifically if (s0K , 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 unT = (s n K , b n K , m n K , g n K ) for all n ∈ N, which is physically admissible: ∀n ≥ 0, ∀K ∈ T , (snK , b n K , m n K , g n K ) ∈ A. Theorem 3.2 (Convergence to a weak solution): If the initial data is physically admissible and the discrete equivalent of the H1 semi-norm [10] of b0 and m0 are bounded, then there exists a subsequence (uh) 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 well- suited 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. Fig. 1: Voronoi diagram used for the simulation of a skull fracture (5884 control volumes) 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, Biomath 1 (2012), 1209254, http://dx.doi.org/10.11145/j.biomath.2012.09.254 Page 3 of 5 http://dx.doi.org/10.11145/j.biomath.2012.09.254 Y. Coudière et al., An Upstream Finite Volume Scheme for a Bone Healing Model like circular fractures for skull bones. We choose to use Voronoi diagrams, that verify the property of orthogo- nality 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 (m0(x, y) = 0.1 g.ml-1) and cellular cluster (s0(x, y) = 106 cells.ml-1 and g0(x, y) = 2 × 103 ng.ml-1). 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 (a) Bone matrix density 0 ≤ m ≤ 0.1 g.ml-1 (b) Concentration of stem cells 0 ≤ s ≤ 2.29 × 103 cells.ml-1 (c) Concentration of osteoblasts 0 ≤ b ≤ 9.9 × 105 cells.ml-1 (d) Concentration of the growth factor 0 ≤ g ≤ 180 ng.ml-1 Fig. 4: Bone matrix density, concentrations of stem cells, osteoblasts and the growth factor at t = 3 days. Biomath 1 (2012), 1209254, http://dx.doi.org/10.11145/j.biomath.2012.09.254 Page 4 of 5 http://dx.doi.org/10.11145/j.biomath.2012.09.254 Y. Coudière et al., An Upstream Finite Volume Scheme for a Bone Healing Model 0 0.173 0 0.05 0.1 x (cm) m (g/ml) 0 day 4 days 8 days 16 days 42 days Fig. 5: Bone matrix density along the line y = x as a function of distance to the origin (x, y) = (0, 0). 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. REFERENCES [1] C. Laurencin, A. Ambrosio, M. Borden, and J. Cooper Jr, “Tissue engineering: orthopedic applications”, Annual review of biomedical engineering, vol. 1 no. 1, pp. 19–46, 1999. http://dx.doi.org/10.1146/annurev.bioeng.1.1.19 [2] C. T. L. D. et Migonney V. B. Bujoli, “Ces matériaux qui réparent le corps”,CNRS Le journal, no. 252-253, p. 26, 2011. [3] A. Salgado, O. Coutinho, and R. Reis, “Bone tissue engineering: state of the art and future trends”, Macromolecular bioscience. vol.4 no. 8, pp. 743–765, 2004. http://dx.doi.org/10.1002/mabi.200400026 [4] L. Geris, A. Gerisch, J. Sloten, R. Weiner, and H. Oosterwyck, “Angiogenesis in bone fracture healing: a bioregulatory model”, Journal of theoretical biology, vol. 251 no. 1, pp. 137–158, 2008. http://dx.doi.org/10.1016/j.jtbi.2007.11.008 [5] S. Shefelbine, P. Augat, L. Claes, and U. Simon, “Trabecular bone fracture healing simulation with finite element analysis and fuzzy logic”, Journal of biomechanics, vol. 38 no. 12, pp. 2440– 2450, 2005. http://dx.doi.org/10.1016/j.jbiomech.2004.10.019 [6] M. Gómez-Benito, J. Garcı́a-Aznar, J. Kuiper, and M. Doblaré, “Influence of fracture gap size on the pattern of long bone healing: a computational study”, Journal of theoretical biology, vol. 235 no. 1, pp. 105–119, 2005. http://dx.doi.org/10.1016/j.jtbi.2004.12.023 [7] A. Bailon-Plaza and M. Van Der Meulen, “A mathematical framework to study the effects of growth factor influences on fracture healing”, Journal of Theoretical Biology, vol.212 no. 2, pp. 191–209, 2001. http://dx.doi.org/10.1006/jtbi.2001.2372 [8] T. A. Einhorn, “The cell and molecular biology of fracture healing”, Clinical Orthopaedics and Related Research, vol. 355, pp. S7–S21, 1998. http://dx.doi.org/10.1097/00003086-199810001-00003 [9] B. McKibbin, “The biology of fracture healing in long bones”, J. Bone Joint Surg. Br., vol.60-B no. 2, pp. 150–62, 1978. [10] R. Eymard, T. Gallouët, and R. Herbin, “Finite volume meth- ods”, Handbook of numerical analysis, vol. 7, pp. 713–1018, 2000. http://dx.doi.org/10.1016/S1570-8659(00)07005-8 [11] G. Cooper, M. Mooney, A. Gosain, P. Campbell, J. Losee, and J. Huard, “Testing the “critical-size”in calvarial bone defects: revisiting the concept of a critical-sized defect (csd)”, Plastic and reconstructive surgery, vol. 125 no. 6, 1685, 2010. http://dx.doi.org/10.1097/PRS.0b013e3181cb63a3 Biomath 1 (2012), 1209254, http://dx.doi.org/10.11145/j.biomath.2012.09.254 Page 5 of 5 http://dx.doi.org/10.1146/annurev.bioeng.1.1.19 http://dx.doi.org/10.1002/mabi.200400026 http://dx.doi.org/10.1016/j.jtbi.2007.11.008 http://dx.doi.org/10.1016/j.jbiomech.2004.10.019 http://dx.doi.org/10.1016/j.jtbi.2004.12.023 http://dx.doi.org/10.1006/jtbi.2001.2372 http://dx.doi.org/10.1097/00003086-199810001-00003 http://dx.doi.org/10.1016/S1570-8659(00)07005-8 http://dx.doi.org/10.1097/PRS.0b013e3181cb63a3 http://dx.doi.org/10.11145/j.biomath.2012.09.254 Introduction A model for Population Dynamics The Numerical Scheme Numerical Simulation: a Skull Fracture Conclusion References