Finite memory model of bone healing in analysis of moving interface between mandible tissue and bone substitute material after tooth implant application

A previous bone remodelling model was presented elsewhere [30], and in the present paper, the same model was tested with new conditions; an interaction between bone tissue, bone substitute material and a dental implant was considered. The bone substitute material was assumed to be dead tissue, which does not synthesizes neither absorbs bone tissue, and it was considered, as well, resolvable. A moving border between the bone substitute material and the bone tissue was studied. The border moved as the newly synthesised bone tissue took over the bone substitute material. After the numerical calculations of time-steps, the whole bone substitute material was replaced by normal bone tissue and the implant was fixed in place only by bone tissue. Dynamical studies of the interaction of bone tissue and implant are used to improved implant design considering different factors, in this case, the presence of bone substitute material helping to fix the implant.

The model was tested against X-rays from dental practice elsewhere [30], and it was shown to be consistent with the actual biomechanical changes in bone's density in different scenarios. In the present study, the model is extended to examine the case where a bone substitute material is used to fill the space between the dental implant and the mandibular bone. That procedure is often used when the quality of the bone is not optimal and bone's density is increased artificially to improve the probability of good implant fixation. Often the actual patient's bone taken from another site in the mouth, or even body, is used as bone's substitute material.

Materials and Methods
Detail considerations about the mathematical simulation of bone remodelling are given elsewhere [30], and an alternative version was given in [22,23,27], but a summary formulation is reported in the following.
It is assumed that an osteocyte feels a density of strain stimulus and produce a signalling E x defined as: where the following notation was introduced: ξ describes the position of the osteocyte, i.e. the point where the signalling is produced; x describes the position of the osteoclast or the osteoblast, i.e. the point where the signalling is received; is the stimulus which depends on the density of strain energy and time; R is the distance from where the signalling is produced; X describes the bone's tissue mass density or bone's quality at the point ξ where the signalling is produced and describes how signalling decays with time, in our case we considered to be exponential. As well, signalling was considered to decay exponentially with the distance from the point where is produced. For a given point, at a given time, the signalling could be written as: where ρ denotes tissue mass density. Biomechanical processes were assumed quasistatic. Therefore, the linear static elastic equations used were: equilibrium equation: Linear strain tensor: General stress-strain tensor of material parameters, C i jkl , was considered to depends on tissue mass density: As a finite number of osteoblast and osteoclast work on the bone, either building or absorbing tissue, an upper and a lower limits for the bone's density were defined. An equilibrium state with a particular density of strain energy associated was assumed to exist in which the bone absorbed is equal to the bone produced. It was considered that the osteocyte emits a signalling ordering an action according to the strain felt, this signalling decayed exponentially to the distance from the osteocyte. Moreover, the signalling produced at any moment decays not only in space but also in time, which means that all past states are considered to determine the actual value of the signalling and stimulus. Osteoclast and osteoblast were considered to react to the total signalling emitted by all the osteocytes and either synthesise or absorb bone tissue or stay without action according to the difference between the biomechanical equilibrium density of strain energy and the actual density of strain energy at the point where the osteoblast or osteoclast is located.
In the bone substitute material region, no bone is neither absorbed not produced as the biomaterial is considered dead tissue. However, at the interface between bone tissue and biomaterial, the biomaterial is absorbed and new bone tissue is produced according to the remodelling rules discussed in the following paragraphs.
The ratio of change in bone's density was considered as an optimization of the difference between a function of the local density of strain energy and the equilibrium density of strain energy. The function of the local density used was the total stimulus E s . The ratio of change in bone's density was calculated using the relationship: Where d F a is the ratio of density change given as percentage and a is a constant that regulated density increment steps. A similar idea was used previously in the literature [31]. The density of strain energy that maintains the equilibrium between absorption and deposit of bone was represented by W 0 , and it was calculated with the following relationship: Where W s is the average of the density of strain energy. The integral is calculated over the whole mandible . W 0 was calculated only in the first cycle, and it was maintained constant in calculations throughout the 150 time-steps.
As it was already mentioned, the stimulus does not disappear instantly but decreases with time, i.e. the stimulus at any time will contribute to the stimulus of the future time-steps. From equation (2), the total stimulus E s is calculated by: With the mathematical justification given elsewhere in more detail [1,5,13,21,24,30] using the Prony series, γ to expand the time-dependent factor in 8 with finite memory of four time-steps, and coefficient of relaxation, γ , equal to 1, the total stimulus could be written as: where α i is a constant. We consider ε to be small enough to be neglected so we keep only the first two terms, which consider the current stimulus plus the stimulus of previous 4 cycles decaying exponentially. The contribution of the 5th element was calculated to affect less than 1% of the total value which is small enough to be neglected. A lazy zone was considered where stimulus is so low that no remodelling takes place. Considering the lazy zone limits, the signalling had the following relationship: very similar to the one used previously in the literature [31], where dzp is the lazy zone upper limit, and the dzn is the lazy zone lower limit. More details about this relationship are given elsewhere [30].
In numerical finite element simulation, the following relation was assumed between the tissue mass density and Young's modulus: where ρ denotes the tissue mass density and E the Young's modulus. The parameter b was chosen to give the values collected in Table 1 for a normal bone characteristics. The model's dimensions and more details concerning the position of the forces considered and other mechanical requirements were taken from the model described previously [30]. The considered 2D model is presented in Figure 2. The boundary between the biomaterial and the bone was defined with the used of 150 points with moving coordinates. For those points, the parameter given in equation 6 was calculated; the negative value indicated motion of the point away of the implant if it was positive towards the implant. To calculate the displacement, we used the following relation: where disp was the displacement of the point, F d a parameter proportional to the total stimulus strength at the border point at a particular time given by equation (9), and γ a parameter to regulate how fast the boundary moves. The actual displacement vector was normalized, and its direction was calculated as the bisector of the angle formed with the two points that were beside the point considered; one to the left and one to the right. Points at the edge of the mandible were considered to move along the mandible border. After the border displacement was calculated, the new density was calculated using equations (6), (9) and (10) for those new points located inside the bone's area. The correspondent Young's modulus was calculated using equation (11).
To perform the numerical testing of the model, the mathematical relations were implemented in a FEM software, COMSOL Multiphysics , to solve a set of differential equations. To control the iterations representing time, MATLAB ver. R2018b was used.
While the mathematical simulations were run at COMSOL Multiphysics , the Dirichlet boundary conditions were kept. The placement of the static forces over the model and other mechanical characteristic of the model are giving elsewhere [30].
The following material characteristics were considered for the mathematical simulation: The Comsol Multiphysics model was controlled by a MATLAB ver. R2018b script according to the flow diagram shown in Fig. 1.

Results
The initial configuration of the model is shown in Fig. 2. After running the model 50 time-steps, the boundary changed to a new position: The bone tissue in the mandible was very close to a biomechanical state of equilibrium at the beginning of the simulation, and no significant changes in bone density were noticed. The small picture at the right portraits a magnified region close to the implant to show in more detail the interface between bone tissue and bone substitute material. After 100 time-steps, the model changed to: At this stage, the border moved the middle position between the initial configuration and the implant. After 150 time-steps, the model configuration changed to: At this stage, the whole bone substitute material was absorbed and replaced by healthy normal bone.

Run the model with initial values
At the first cycle, calculate equilibrium parameters and feed them back to Comsol.
Save data of the new calculated tissue mass density and coordinates of the border between the bone substitute material and bone tissue.
Input the new tissue mass density and associated material parameters as well as the coordinates of the new border between bone substitute material and bone tissue as initial values.
Run the model in Comsol.
Chech convergence and stop calculations if changes of the model are negligible. Finite element method has been used to examine interactions between bone and dental implants [4,7,12,17,29,33,37]. Some authors concentrated on the prosthesis design [9,29]. Others focused on the stress distribution in the bone after the implant [4,11]. The distribution of stress on the implant was also tested [33], and the interaction between bone and implant was as well studied [7,37].
In the present paper, the relationship between bone tissue, bone substitute material and an implant was investigated. The border between the bone substitute material and bone tissue was considered dynamic in the sense that bone tissue, due to a stimulus to grow, takes over the bone substitute material. In our investigation, the bone effectively advances over the bone substitute material and fix the implant in place.
This examination of the bone tissue, bone substitute material and implant complements our previously presented model and proves finite element method and mathematical modelling as a strong tool to study complex and dynamic cases of bone remodelling in the presence of an implant and a bone substitute material.  As noted previously by some authors, the future design of implant should consider actual conditions of patient's bone to improve osteointegration [32,37]. More investigation is needed on the influence of several mechanical parameters of the implant design on osteointegration in the presence of bone tissue weakened by disease, i.e. osteoporosis, osteoarthritis, etc [15,[34][35][36].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.