ANALYSIS OF A BONE REMODELING MODEL

. In this work we study, from the numerical point of view, a bone remodeling model. The variational formulation of this problem is written as an elliptic variational equation for the displacement ﬁeld, coupled with a ﬁrst- order ordinary diﬀerential equation, with respect to the time, to describe the physiological process of bone remodeling. Fully discrete approximations are in- troduced, based on the ﬁnite element method to approximate the spatial variable, and on an Euler scheme to discretize the time derivatives. Error estimates are obtained on the approximate solutions, from which the linear convergence of the algorithm is derived under suitable regularity conditions. Finally, some numerical results, involving examples in one, two and three dimensions, are presented to show the accuracy and the performance of the algorithm.

1. Introduction. In this paper we introduce an analyze a numerical method to solve a mathematical model for the bone remodeling process. This model, derived by Cowin and Hegedus ([3,4,13]), is a generalization of the nonlinear elasticity, and it is based on the fact that the living bone is continuously adapting itself to external stimuli. Since this process has an enormous effect on the overall behavior and health of the entire body, the ability of these models to predict the bone remodeling is of great importance.
During the last ten years, some papers dealt with mathematical issues of these models as the existence and uniqueness of weak solutions under some quite strong assumptions (see, e.g., [14,15,16,17,19]), the analysis of an asymptotic rod model ( [6,7,8,9]) or the numerical stability of finite element models ( [11]). Here, our aim is to continue [16], providing the numerical analysis of a fully discrete algorithm and performing some numerical simulations which demonstrate its accuracy and its behavior.
The paper is outlined as follows. In Section 2 we describe briefly the mechanical problem and we derive its variational formulation. An existence and uniqueness result for this problem is also recalled. Then, a fully discrete scheme, based on the finite element method to approximate the spatial variable and on an Euler scheme to discretize the time derivatives, is introduced in Section 3. A main error estimates result is proved from which, under suitable regularity conditions, the linear convergence of the algorithm is deduced. In Section 4, a numerical algorithm to solve the fully discrete problem, is described and some numerical simulations, 2. Mathematical formulation of the Cowin-Hegedus model. Let us denote by · the euclidean inner product in R d and its corresponding norm by | · |. Let S d be the space of second order symmetric tensors on R d , or equivalently, the space of symmetric matrices of order d, and let : be its inner product and | · | its norm.
Denote by Ω ⊂ R d , d = 1, 2, 3, an open bounded domain and denote by Γ = ∂Ω its outer surface which is assumed to be Lipschitz continuous and it is divided into two disjoint parts Γ D and Γ N . We denote by x = (x i ) d i=1 a generic point of Ω = Ω ∪ Γ, and for x ∈ Γ, ν(x) = (ν i (x)) d i=1 denotes the outward unit vector, normal to Γ at point x (see FIGURE 1).
The body occupying the volume Ω (in this work a bone) is being acted upon by a volume force of density f , it is clamped on Γ D and surface tractions with density g act on Γ N . Finally, let [0, T ], T > 0, be the time interval of interest. Let u : (x, t) ∈ Ω × [0, T ] → u(x) = (u i (x, t)) ∈ R d be the displacement field, σ : (x, t) ∈ Ω × [0, T ] → σ(x) = (σ ij (x, t)) ∈ S d the stress field and ε(u) : (x, t) ∈ Ω × [0, T ] → ε(x, t) = (ε ij (u)) d i,j=1 ∈ S d the strain tensor given by ε ij (u(x, t)) = 1 2 To measure the change of volume fraction from a reference configuration we introduce a function e : (x, t) ∈ Ω × [0, T ] → e(x, t) ∈ R, which we will name as bone remodeling function. According to [4,13], the bone is assumed elastic and the constitutive law is written as follows, where ξ 0 represents the reference volume fraction and B(e) = (B ijkl (e)) d i,j,k,l=1 is a constitutive tensor depending on the bone properties that will be described below. We notice that if ξ 0 = 1 and e = 0, the constitutive law corresponds to the classical Hooke's law for linear elastic materials.
The evolution of the bone remodeling function is obtained from the following first-order ordinary differential equation (see [4,13]), where a(e) is a constitutive function and A(e) = (A ij (e)) d i,j=1 denote the bone remodeling rate coefficients. A dot above a function denotes the time derivative of this function.
We notice that the functions B(e), a(e) and A(e) characterize the material properties and there are few experimental data concerning their form. We adopt the polynomial approximations employed in several papers (see, for instance, [13]): ij are constants representing the bone properties.
For a given constant L > 0, let us define the truncation operator Finally, the process is assumed quasistatic and therefore, the inertia effects are neglected. Moreover, let e 0 denote the initial bone remodeling function.
The mechanical problem, derived from the continuum mechanics laws in the framework of small displacements theory, is the following (see [13]).
Problem P. Find the displacement field u : Ω × [0, T ] → R d , the stress field σ : Ω × [0, T ] → S d and the bone remodeling function e : Ω × [0, T ] → R such that e(0) = e 0 and, Here, γ > 0 is the density of the full elastic material which is assumed constant for the sake of simplicity.
We turn now to obtain a variational formulation of problem P. We use the classical spaces L 2 (Ω) and H 1 (Ω) equipped with the classical norms (see [1]): the variable x i . Let us denote by Y = L 2 (Ω) and H = [L 2 (Ω)] d , and define the following spaces equipped with the product norms: The following hypotheses are assumed on the problem data in order to obtain an existence and uniqueness result (see [16]).
The elasticity coefficients B ijkl are assumed to satisfy the following properties: (a) There exists L B > 0 such that The constitutive function a and the bone remodeling rate coefficients A ij are Lipschitz and bounded functions in R; that is, there exist L a , L A , M a and M A such that: The reference volume fraction ξ 0 satisfies the following conditions for some 0 < ξ m 0 < 1, The density forces have the regularity, and the initial value of the bone remodeling function e 0 verifies that e 0 ∈ C(Ω).
For every e ∈ L ∞ (Ω), let us define the following bilinear form c(e; ·, ·) : and the linear form L(e; ·) : V → R given by Throughout this work we systematically use the following identification for every function s : Ω × [0, T ] → R. Applying Green's formula, we derive the following variational formulation for the mechanical problem P.
Problem VP.
e(t) = a(e(t)) + A(e(t)) : ε(u(t)) in D ′ (0, T ; L 2 (Ω)), where D ′ (0, T ; L 2 (Ω)) is the space of distributions valued on L 2 (Ω) (see [1]). The following result, which states the existence of a unique weak solution to Problem VP, is obtained in [16] by using a priori estimates, the Cauchy-Lipschitz-Picard and Schauder theorems (full details can be seen in [18]). Theorem 2.1. Let the assumptions (6)-(10) hold. Assume that, for any given function e ∈ C 1 ([0, T ]; C(Ω)), the unique solution to the following problem: 3. Numerical analysis of a fully discrete scheme. We now introduce a finite differences scheme and a finite element algorithm for approximating solutions to Problem VP and derive an error estimate. First, we consider two finite dimensional spaces V h ⊂ V and B h ⊂ L ∞ (Ω) ⊂ Y , approximating the spaces V and L ∞ (Ω), respectively. Here, h > 0 denotes the spatial discretization parameter.

Remark 1.
In the numerical simulations presented in the next section, V h and B h consist of continuous and piecewise affine functions and piecewise constant functions, respectively; that is, where Ω is assumed to be a polygonal domain, T h denotes a finite element triangulation of Ω composed by d-symplex denoted by T r, and P q (T r), q = 0, 1, represents the space of polynomials of total degree less or equal to q in T r.
Secondly, the time derivatives are discretized by using a uniform partition of the time interval [0, T ], denoted by 0 = t 0 < t 1 < . . . < t N = T where t j = jk, k being the time step size defined as k = T /N and N ≥ 1. Moreover, for a continuous function S(t) we let S n = S(t n ).
In this section, C denotes a positive constant which depends on the problem data and the continuous solution, but it is independent of the discretization parameters h and k.
Using the forward Euler scheme, we propose the following fully discrete approximation of Problem VP: where e hk 0 ∈ B h is an appropriate approximation of the initial condition e 0 . From properties (6) it is straightforward to obtain the existence and uniqueness of the discrete solution which we state in the following.
Below we focus our attention on deriving estimates for the numerical errors u n − u hk n V and e n − e hk n Y . The following theorem is the main result of this work.
Proof. First, let us obtain an error estimates on the bone remodeling function e. By substraction of (12) at time t = t n and (17) we find thaṫ Integrating in Ω and using the norm in Y , we have, Now, using (7) we get and a(e n ) − a(e hk n−1 ) Y ≤ C( e n − e n−1 Y + e n−1 − e hk n−1 Y ). Taking into account the regularity e ∈ C 1 ([0, T ]; C(Ω)) ⊂ C 1 ([0, T ]; Y ) and applying the mean value theorem we have e n − e n−1 Y ≤ k ė C[0,T ;Y ] . Now, by an induction argument we conclude the following estimate on the bone remodeling function e, Next, let us estimate the numerical errors on the displacement field. Thus, we write equation (11) for all v = v h ∈ V h ⊂ V and we substract it to equation (16) to obtain, c(e n ; u n , u n − u hk n ) − L(e n ; u n − u hk n ) − c(e hk n ; u hk n , u n − u hk n ) + L(e hk n ; u n − u hk n ) = c(e n ; u n , Keeping in mind that , we can write the previous equation as, V . Then, we obtain the next inequality, From properties (6) and (9) it follows that c(e hk n ; u n , u n − u hk n ) − c(e n ; u n , u n − u hk n ) n Y u n − v h V , and therefore we can deduce that, Using several times the inequality, taking ǫ small enough we find, for all v h ∈ V h , Combining (19) and (21), we have, for all v h n ∈ V h . Using now a discrete version of Gronwall's inequality (see [12] for details), it leads to (18).
Error estimates (18) are the basis for the analysis of the convergence rate of the algorithm, presented below using the finite element method.
Let Ω be a polyhedral domain and denote by T h a finite element triangulation of Ω compatible with the partition of the boundary Γ = ∂Ω into Γ D and Γ N . Let V h and B h be defined by (14) and (15), respectively, and assume that the discrete initial condition e hk 0 is given by e hk 0 = π h e 0 , where π h : C(Ω) → B h is the standard finite element interpolation operator (see, e.g., [2]).
Assume the following additional regularity conditions on the continuous solution: The next result follows from estimates (18).
Proof. First, we recall the well-known approximation property of the finite element space V h (see [2]), Using the mean value theorem, it is easy to check that From the definition of the interpolation operator π h we have (see again [2]), Using again the approximation property of the finite element space V h and taking into account that Then, from (21) with n = 0 and (26) we have: The estimation (23) is now concluded from (18) (14) and (15), respectively. First, u hk n is obtained solving problem (16): u hk n ∈ V h , c(e hk n ; u hk n , v h ) = L(e hk n ; v h ) ∀v h ∈ V h . This is a linear problem equivalent to a linear system which is solved by using Cholesky's method.
Next, let n ∈ {1, . . . , N } and suppose that u hk n−1 and e hk n−1 are known. The discrete bone remodeling function e hk n is calculated from (17) as: e hk n = e hk n−1 + ka(e hk n−1 ) + kA(e hk n−1 ) : ε(u hk n−1 ). The numerical scheme was implemented on a 3.2Ghz PC using MATLAB, and a typical 1D run (h = k = 0.001) took about 3.5 seconds of CPU time. A run for the 2D example spent about 8.5 seconds for each time iteration and for the 3D example took about 2 minutes for each time iteration.

4.2.
Numerical results for one-dimensional problems.

4.2.1.
A first one-dimensional example: numerical convergence. As a one-dimensional example, the following problem is considered.
We observe that functions B(e) and a(e) do no satisfy the boundedness assumptions presented in (6) and (7). However, it is easily solved by using the truncation function Φ L . Anyway, we use value L = 10 6 , so it is large enough and we can assume that this truncation does not modify the results. Using discretization parameters h = k = 0.001, the pointwise errors of the displacement and bone remodeling fields at times t = 0.2, 0.4, 0.6, 0.8, 1 s. are plotted in FIGURE 2. As it can be seen, the highest pointwise errors are concentrated near the right end. In FIGURE 3, the evolution in time of the errors of the displacement and bone remodeling fields at point x = 1 is shown. As it was expected, the errors increase with respect to the time.   In FIGURE 6 the bone remodeling and displacement fields are plotted at different times (t = 0, 2, 5, 12, 20, 100 days). We can see that, the bone remodeling function increases with time and, consequently, converges to the homogeneous steady case (e ≈ 0.021677), the stiffness is larger and the displacements decrease. u(x,t) x Displacements t=0 t=2 t=5 t=12 t=20 t=100 Figure 6. Example 1D-2: Bone remodeling function and displacement field at several times.
Assuming now that the initial bone remodeling function vanishes (e 0 = 0) and that the lower corner is fixed, in FIGURE 7 the bone remodeling function and the displacements are shown at different times. Now, we can observe that the bone remodeling function is always homogeneous and it converges to a steady state. The displacements also decrease because the stiffness increases.  where the fourth-order tensors B 0 = (B 0 ijkl ) 2 i,j,k,l=1 and B 1 = (B 1 ijkl ) 2 i,j,k,l=1 and the second-order tensors A 0 = (A 0 ij ) 2 i,j=1 and A 1 = (A 1 ij ) 2 i,j=1 have the following components (see [10]): Using the time discretization parameter k = 0.01, in FIGURES 9 and 10 the displacements (multiplied by 20) and the bone remodeling function are shown at initial time (left) and after twenty days (right). As can be observed, the displacements decrease since the bone remodeling function is positive and so the stiffness increases. Finally, we notice that the bone remodeling function tends in time to the steady and homogeneous case (e ≈ 0.03156).

4.3.2.
Second two-dimensional example: two linearly increasing compression forces. As a second two-dimensional example we consider a similar setting than in the previous case (see FIGURE 11). The unique differences are: there is not initial bone remodeling (e 0 = 0), the final time is now T = 108 days and the compression forces are supposed to be linearly increasing with respect to x and they have the following form: Taking k = 0.01 as the time discretization parameter, the displacements fields (multiplied by 20), at initial time and at final time, and the bone remodeling function at final time are plotted in FIGURES 12 and 13, respectively. As we can see, the body bends and there is a recuperation at long time in the displacements because of the bone remodeling. Moreover, the bone remodeling function is positive on the right part since the body is compressed there and it is negative on the left one since there is an extension there. 4.3.3. Third two-dimensional example: a linearly increasing compression force. As a final two-dimensional example, we consider a similar physical setting than in the previous case. Here, we assume that the final time is now T = 60 days and that only a linearly increasing compression force acts on the upper horizontal part (with Ω Γ N Γ N g g Figure 11. Example 2D-2: Physical setting. the same expression that in the above example). The lower horizontal boundary remains clamped (see FIGURE 14).
Using the time discretization parameter k = 0.01, in FIGURES 15 and 16 the displacements (multiplied by 20), at initial time and after 60 days, and the bone remodeling function at final time are shown. Again, we can observe that there is a decrease in the displacements because of the bone remodeling. Moreover, the bone remodeling function takes the highest values where the body compresses.     Taking k = 0.01 as the time discretization parameter, the initial configuration and the displacement field (multiplied by 50), at initial time and at final time, and the bone remodeling function, at final time, are shown in FIGURES 19 and 20. As we expected, the body bends along the direction of the applied forces. Figure 19. Example 3D-2: Initial configuration and displacements (×50) at initial time (left) and after 85 days (right).
As it happened in the previous case and also in the two-dimensional examples, we notice that the deformation decreases through the time due to the bone remodeling, and the bone remodeling function takes positive values where the compression is produced and negative ones where there is extension.

5.
Conclusions. The paper dealt with the numerical analysis, including numerical simulations in one, two and three dimensions, of a bone remodeling model introduced in [4]. A numerical algorithm for the variational problem, based on the finite element method to approximate the spatial variable and an Euler scheme to discretize the time derivatives, was proposed, an error estimate on its solutions was obtained and its linear convergence was established under suitable regularity assumptions.
The algorithm was implemented using MATLAB code and some examples were computed. First, a one-dimensional setting was choosen in such a way as to show the numerical convergence of the algorithm. As can be seen in FIGURE 4, the linear convergence of the algorithm was achieved. Then, two similar one-dimensional examples were performed under compression forces. In both cases it was observed that there was a decrease in the displacements since the stiffness increases along with the bone remodeling function. Secondly, three two-dimensional examples were presented. In the first example, constant compression forces were applied and a constant bone remodeling function was obtained. Then, two examples including compression forces were provided from which a decrease in the displacements was observed due to the bone remodeling. Finally, two three-dimensional examples were described, assuming that one or two linearly increasing compression forces acted on the body. The same phenomena were reproduced. All of them allowed us to observe the good behavior of the model and the numerical scheme.