Numerical Analysis of an Osseointegration Model

In this work, we study a bone remodeling model used to reproduce the phenomenon of osseointegration around endosseous implants. The biological problem is written in terms of the densities of platelets, osteogenic cells, and osteoblasts and the concentrations of two growth factors. Its variational formulation leads to a strongly coupled nonlinear system of parabolic variational equations. An existence and uniqueness result of this variational form is stated. Then, a fully discrete approximation of the problem is introduced by using the finite element method and a semi-implicit Euler scheme. A priori error estimates are obtained, and the linear convergence of the algorithm is derived under some suitable regularity conditions and tested with a numerical example. Finally, oneand two-dimensional numerical results are presented to demonstrate the accuracy of the algorithm and the behavior of the solution.


Introduction
The study of the phenomenon of osseointegration in endosseous implants is a very interesting issue because of the increasing use of many types of implants in clinical practice. Dental implantation stands out as the area that has benefited the most from the innovation and continuous development of bone implants in the last few years. This is probably due to the outstanding clinical results (see, e.g., [1]). These dental implants are artificial systems that consist of an endosteal component, which is completely inside the mandible, and an abutment that connects the endosseous component with the oral cavity, in order to replace a missing tooth [2].
The behavior of bone tissue depends on the chemical properties of the materials used on the implant and on the ability of the implant to introduce a mechanical state of stress, which induces the osteogenesis processes. Therefore, there are two possible ways to study this process. The first one involves clinical trials, where known failure mechanisms of implants are evaluated and their incidence assessed. The second one consists of preclinical studies, which allow performing a preliminary analysis of the system using computer tools.
Currently, the number of papers dealing with the numerical simulation of the behavior and stability of dental implants is huge (see, for instance, ). This work also follows the aforementioned research line about preclinical studies, continuing the studies already published in [27,28]. There, the authors presented a new biological model to predict the osseointegration of the implant in the bone. In the first part [27], they considered the full model, performing two-dimensional numerical simulations (implemented in ABAQUS) of the bone ingrowth process around the dental implant, assuming two different surface properties. Then, in the second part [28], two simplified versions of this problem were considered, providing a stability analysis and showing the existence of traveling wave-type solutions. In this work, we also continue the research started in [29,30] (where both simplified models were numerically analyzed) considering the full osseointegration model.
The paper is structured in the following way. The biological problem and its variational formulation are presented in Section 2. Then, fully discrete approximations are introduced in Section 3. We use the finite element method to approximate the spatial variable and an Euler scheme to discretize the time derivatives. A priori error estimates are then proven, from which the linear convergence is derived. Finally, one-and two-dimensional numerical simulations are presented in Section 4, where we show the accuracy of the numerical resolution and the behavior of the solution.

Biological Problem and Its Variational Formulation
In this section, a brief description of the model is presented. To obtain further details, we refer to the paper where the model was presented [27]. Denote by "·" and · the inner product and the Euclidean norm on R d . Let Ω ⊂ R d , d = 1, 2, 3, be a domain with a Lipschitz boundary Γ = ∂Ω, which can be decomposed into two disjoint parts Γ 1 and Γ 2 such that measure(Γ 1 ), measure(Γ 2 ) > 0. Let [0, T], T > 0, and ν = (ν i ) d i=1 represent the time interval of interest and the outward unit normal vector to Γ, respectively. Finally, let x ∈ Ω be the spatial variable and t ∈ [0, T] the time variable. We do not indicate the dependence of the functions on them in order to simplify the writing.
According to [27], we describe the biological model in what follows. Once the surgical procedure is finished and the implant is placed, blood fills the cavity, and proteins and tissue fluids are absorbed into the implant surface. Then, platelets are activated due to the contact with the surface, and they release a number of growth factors. These growth factors stimulate the migration and proliferation of osteoblasts. Once these processes begin, the osteogenic cells coming from the old bone migrate to the surface of the implant, and they differentiate into osteoblasts, which later synthesize bone matrix.
Here, following [27], we assume that there are two generic growth factors. The first factor corresponds to the release of activated platelets (PGDF, TGF-β), and the second one accounts for the molecules secreted by osteoblasts and osteogenic cells (BMPs, TGF-β superfamily).
First, let us denote by c ∈ R the density of platelets, whose evolution is determined as follows, The evolution of c was modeled as a linear diffusion (multiplied by D c ), the cell adhesion to the implant surface, and a kinetic term. Cell adhesion was modeled with a linear "taxis" term depending on the gradient of adsorbed proteins ∇p, with a coefficient H c . We note that the value of p is assumed to be known, and it depends on the microtopography of the implant surface. Its maximum value appears at the surface of the implant and then decreases as we move away from this surface, reaching soon a value of zero that is maintained in the remaining domain [27]. We assume a high platelet concentration at the beginning, so the only kinetic term appears from cell removal due to inflammation, with a linear rate A c .
Secondly, the density of the osteogenic cells, m, is modeled by the equation: Here, besides the linear diffusion and kinetics, there is a chemotaxis term along the gradients of the growth factors s 1 and s 2 and logistic growths with these factors. These growths are bounded by a limiting cell density N. Further explanations on the modeling can be found in the original reference.
Next, the density of osteoblasts b is modeled by the following equation: In this case, we obtain an ordinary differential equation instead of a partial differential one because osteoblasts remain on the surface of the bone matrix, so there is no flux of these cells. The source term comes from the differentiation from the osteogenic phenotype, and the decay accounts for differentiation into osteocytes A b .
Finally, following [27], we describe the evolution of both generic growth factors s 1 and s 2 . The first generic factor s 1 is assumed to satisfy the following constitutive partial differential equation: Again, random dispersion of this factor is modeled with a diffusive term with a linear coefficient D s1 . The kinetic term accounts for the secretion of s 1 by platelets (c) together with the concentration of the absorbed proteins p and the growth factor s 1 itself, both as logistic terms. As in the previous equations, a natural decay of the factor with linear rate A s1 appears.
The second generic factor is assumed to satisfy a similar equation: Now, there are two source terms depending on logistic functions of s 2 and concentrations of m and b. The natural decay is multiplied by the factor A s2 .
We note that these constitutive equations are nonlinear, and so, for mathematical reasons, some of the terms must be truncated. In order to do that, for some constants M 1 > 0 and M 2 > 0, we introduce the truncation operators R 1 : R → [0, M 1 ] and R 2 : R → [0, M 2 ] given by: Remark 1. We point out that these truncation operators are introduced to assure some mathematical properties. Anyway, we note that, although we have not been able to prove these boundedness properties theoretically, in the bone remodeling processes, as well as in the performed numerical simulations, they were observed.
In order to simplify the writing, we define the following functions f : R 3 → R, g : R 2 → R, and h : R 3 → R: Concerning the boundary conditions for the density of the osteogenic cells m, we assume a known value on Γ 1 since there is a contribution from the bone. There is no flux through Γ 2 . Moreover, we assume that there is no flux through the whole boundary Γ for the remaining variables (density of platelets and generic growth factors).
Finally, we denote by c 0 , m 0 , s 10 , s 20 , and b 0 the initial conditions for c, m, s 1 , s 2 , and b, respectively. Collecting all these equations, the full bone integration problem is then written as follows.
∂b ∂t ∂m ∂ν Remark 2. We note that in [27], the authors also derived the equations for the volume fractions of the fibrin network, the woven bone, and the lamellar bone, but they were obtained from the above densities and growth factors. Therefore, in this work, we did not consider them. However, it could be interesting to analyze the evolution of the mass density [31][32][33].
To obtain the variational formulation of Problem P, let Y = L 2 (Ω), H = [L 2 (Ω)] d and E = H 1 (Ω), and denote by (·, ·) Y , (·, ·) H , and (·, ·) E their respective scalar products, with corresponding norms · Y , · H , and · E . Moreover, let us define the variational space V as follows, with scalar product (·, ·) V and norm · V . Integrating in time Equation (5) and using the initial condition (9), we obtain: Plugging Equation (10) into Equation (4), applying Green's formula, and then using the boundary conditions (6)-(8), we derive the following variational formulation of Problem P, written in terms of densities c and m and concentrations s 1 and s 2 .
In the following result, we state that the above problem has a unique weak solution.
Theorem 1. If we assume that the constitutive coefficients D c , A c , D m , A m , D s1 , A s1 , D s2 , and A s2 are strictly positive, then Problem VP admits a unique solution with the regularity: The proof of the above theorem is rather technical, and it is based on well known results on evolution equations with monotone operators (see, for instance, [34][35][36]) and an application of the fixed-point theorem. We omit the details for the sake of reading.

Fully Discrete Approximations and an a Priori Error Analysis
In this section, we introduce a finite element algorithm to approximate solutions to Problem VP. Then, we also provide an a priori error analysis.
The discretization of Problem VP is performed as follows. First, we assume Ω to be a polyhedral domain, and we consider two finite dimensional spaces E h ⊂ E and V h ⊂ V, which approximate the variational spaces E and V, given by: where P 1 (K) represents the space of polynomials with a global degree less than or equal to one in K.
Following [37], we denote by (T h ) h>0 a regular family of triangulations of Ω, which is compatible with the partition of the boundary Γ = ∂Ω into Γ 1 and Γ 2 . Therefore, the finite element spaces E h and V h are composed of functions that are continuous and piecewise affine. Let h K be the diameter of an element K ∈ T h so h = max K∈T h h K denotes the spatial discretization parameter. Moreover, we assume that the discrete initial conditions, denoted by c h 0 , m h 0 , s h 1 0 , and s h 2 0 , are given by: where P h is the L 2 (Ω)-projection operator over E h .
To obtain a discretization of the time derivatives, let 0 = t 0 < t 1 < . . . < t N = T be a uniform partition of the time interval [0, T], and we define k as the time step size given by k = T/N. If r(t) is a continuous function, then we define r n = r(t n ), and for a sequence {w n } N n=0 , we let δw n = (w n − w n−1 )/k be its usual divided differences.
Therefore, using a Euler scheme, we obtain the following fully discrete approximation of Problem VP.
Problem VP hk . Find the discrete density of platelets c hk = {c hk n } N n=0 ⊂ E h , the discrete density of osteoblasts m hk = {m hk n } N n=0 ⊂ V h , the discrete concentration of the first growth factor s hk and the discrete concentration of the second growth factor s hk , and s hk 2 0 = s h 2 0 , and, for n = 1, . . . , N and for all u h ∈ V h and v h , where the discrete density of osteoblasts b hk = {b hk n } N n=0 is given by: We note that we have used the explicit expression in the nonlinear terms in order to simplify the numerical resolution of this discrete problem.
Since the resulting equations of Problem VP hk are now uncoupled because they can be solved separately, using Lax-Milgram lemma, it is easy to obtain the existence of a unique discrete solution we derive some a priori estimates for the numerical errors c n − c hk n , m n − m hk n , s 1 n − s hk 1 n s 2 n − s hk 2 n , and b n − b hk n . Thus, we assume that Problem VP has a unique solution with the following regularity: For the sake of simplicity, we assume that all the constants involved in Problems VP and VP hk are equal to one; i.e., We note that it is immediate to extend the results presented below to the general case.
First, we obtain some estimates on the density of osteoblasts. Then, we subtract Equation (10) at time t = t n and Equation (22) to obtain: where, here and in what follows, C > 0 is a positive constant that depends on the data of the problem and the continuous solution; however, it is independent of the discretization parameters h and k, and I n is the integration error given by: Then, we obtain some error estimates for the density of the platelets c. Finally, we subtract the variational Equation (11) at time t = t n and the discrete variational Equation (18) for a test function v = v h ∈ E h , obtaining: Thus, Using the following property of the divided differences: where the notation δc n = (c n − c n−1 )/k represents the corresponding divided differences, and applying several times Young's inequality: we have, for all v h ∈ E h , Thus, by induction, it follows that, for all where we used the regularity c ∈ C 1 ([0, T]; Y). Now, we obtain the estimates on the density of osteoblasts m. Therefore, subtracting the variational Equation (12) at time t = t n for a test function v = v h ∈ V h ⊂ V and the discrete variational Equation (19), we find that: Taking into account that: where we have used the notation δm n = (m n − m n−1 )/k, the properties of the truncation operators R 1 and R 2 , and the expression of function f , and applying several times inequality (26), it follows that, for all u h ∈ V h , Summing up to n, keeping in mind that m ∈ C 1 ([0, T]; Y), we find that, for all {w h j } n j=0 ⊂ V h , Finally, we estimate the numerical errors on the concentration of the first and second growth factor, s 1 and s 2 , respectively. Then, subtracting the variational Equation (13) at time t = t n for a test function w = w h ∈ E h ⊂ E and the discrete variational Equation (20), we have, for all w h ∈ E h , Therefore, we find that: Taking into account that: where we have used the notation δs 1 n = (s 1 n − s 1 n−1 )/k, the properties of the truncation operator R 1 , and the expression of function g, and applying several times inequality (26), it follows that, for all Therefore, summing up to n, since Proceeding in a similar way, we obtain the following estimates for concentration s 2 , for all Combining now the estimates (24)- (30) and keeping in mind that (see [38]): using the discrete version of Gronwall's inequality (see [39]), we reach the following a priori error estimate result.

Theorem 2.
Let the assumptions of Theorem 1 hold. Assume that Problem VP has the additional regularity (23), and let us denote by (c, m, s 1 , s 2 , b) and (c hk , m hk , s hk 1 , s hk 2 , b hk ) the solutions to problems VP and VP hk , respectively. Therefore, we obtain that there exists a positive constant C > 0, assumed to be independent of the discretization parameters h and k, but depending on the solution (c, m, s 1 , s 2 , b) and the problem data, such that, for all {u h n } N n=0 where integration error I j was defined in (25).
We notice that the error estimates (31) can be used to analyze the convergence rate of the algorithm. Hence, under some additional regularity conditions, we obtain the linear convergence of the algorithm that is stated in the next corollary.
Once these densities and growth factors are calculated, the discrete density of osteoblasts at time t n , b hk n , is obtained from (22). This equation can be expressed in terms of the solution in the previous step, which is more convenient since the summation over all previous steps is avoided, resulting in: This discrete problem leads to an uncoupled linear system of equations. It is solved separately by using the classical Cholesky method. The corresponding numerical algorithm was implemented on a 3.2 GHz PC using FEniCS [40,41] for the 1D simulations, where a typical 1D run (h = k = 0.01) took about 0.84 seconds of CPU time. The two-dimensional problem was solved using FreeFem [42].

Numerical Convergence
In order to show the convergence of the algorithm proven in Corollary 1, a one-dimensional example was solved with all parameters set to one: The function p was considered to take the form p(x) = 10 x − 1, and the initial conditions c 0 = m 0 = s 1 0 = s 2 0 = b 0 = 1. Here, to simplify the description of the example, homogeneous Neumann conditions were employed for all the variables.
Since the exact solution for the coupled problem is unknown, the solution with a refined mesh (h = k = 3 × 10 −5 ) was considered as the reference. The norm of the numerical errors defined in Corollary 1 is shown in Table 1, and the main diagonal of this table is plotted in Figure 1. As we can see, the convergence of the numerical algorithm was found, and the linear convergence, stated in the previous section, seemed to be achieved.

One-Dimensional Examples
Since, to the best of the authors' knowledge, there are no one-dimensional examples regarding this model in the literature and they can provide some insights with an easy spatial domain, some examples are shown. The domain represents a horizontal section of the computational domain presented in [27] (the same domain used for the two-dimensional examples in this section); thus, Ω = [0, 0.6]. The simulation ran until a final time of one day was reached.
First, we address the definition of function p(x). As explained in [27], this function must reach its maximum value (which is a model parameter defined as p m from now on) on the implant surface. Furthermore, it decreases fast to the value of zero (in [27], the function decreases to zero at a distance of 0.1 mm from the implant surface), remaining zero in the rest of the domain. If the support of the function is considered linear, a non-differentiable point appears when the zero value is reached. This could lead to problems since the gradient of this function appears in Equation (1). To avoid problems when the spatial discretization is performed, it is better to create the mesh in such a way that this singular point is not inside an element. This precaution is easy to fulfil in one dimension, but not as easy in a two-or three-dimensional domain.
The results for the concentration of platelets (c) and both growth factors (s 1 and s 2 ) are shown in Figures 3 and 4, respectively. For the case of platelets, it can be seen that both solutions were qualitatively (and to a great extent, also quantitatively) similar. However, for the case of the piecewise definition of p, a non-differentiable point appeared at x = 0.5, while in the other case, the solution was smoother. Regarding growth factors, they showed again a similar behavior, slightly differing for values close to x = 0.6.   To complete the one-dimensional analysis, we show the difference for the two values of p m , as discussed in [27]. We solved the same problem as before, but now comparing the value for p m = 0.5 and p m = 0.1. We chose the sigmoid definition for p in this case.
The results for these simulations are shown in Figure 5 (concentration of platelets) and Figure 6 (concentrations of osteogenic cells and osteoblasts). The results showed similar behavior as that found in [27], the shown variables reaching similar values in the left part of the domain (where the values of p were similar) and greatly differing in the right hand side, where the gradient of p was much lower.
Again, we compared the cases of p m = 0.1 and p m = 0.5. The obtained results are shown in Figure 8 (platelets density), Figure 9 (osteogenic cell density), and Figure 10 (osteoblasts density). These results agreed qualitatively with published results of the literature and were consistent with the 1D computations.   As in the one-dimensional case, when p m = 0.1, the concentration of cells near the implant was much smaller in general. It is also worth noting that for the osteogenic cells (m), the concentrations were higher in the upper and lower implant boundaries than in the right one. This effect also appeared with the osteoblasts, but more accentuated.

Conclusions
In this paper, a numerical analysis of a bone remodeling model for endosseous implants was performed. The variational form was obtained, and from that, an algorithm to solve the discrete problem was proposed. The full discretization was done using the well known finite element method and a combination of Euler schemes. This allowed obtaining mathematical results on the convergence of the algorithm. This convergence was exemplified with a one-dimensional simple test example. Some other one-dimensional cases were solved and completed with a two-dimensional solution.
We point out that this work was a first step toward dealing with the study of these osseointegration models. A possible extension could be to consider, for instance, the so-called joint congruence, that is to analyze how the articular surfaces mate each other (see, for instance, [43,44]).