Novel description of bone remodelling including finite memory effect, stimulation and signalling mechanisms

A new mathematical model is presented for bone remodelling that includes a finite memory effect. In this new model, stimulus resulting from mechanical loading is separated from the signalling to grow or absorb bone. Also, a signal decaying exponentially to the distance from the point produced as well as an effect of decaying signal in time are considered. In addition, the model presented correctly predicts the synthesis and resorption of the bone tissue in a normal healthy mandible and in cases when teeth are missing. The model presented could be implemented to study cases of bone diseases where the signalling between the cells is disrupted and to predict changes in bone caused by several anomalies, such as missing a tooth or the presence of a dental prosthesis.

bone's density were considered. Another factor considered that influence stimulus is the amount of strain felt by the osteocyte [29,30,41,49,65]. The osteocytes emit signals that activate either osteoblasts or osteoclasts, which means bone synthesis, if the signal is positive, and bone resorption if the signal is negative. A "lazy zone" proposed in the article by Huiskes previously cited will be considered, where signal is not enough strong for cells to start working on depositing or absorbing bone. The "lazy zone" is considered to be related to a biomechanical equilibrium state that was assumed to exist with a particular density of strain energy associated called reference energy. The idea is not new as it is used either implicitly nor explicitly in the literature see for example [36,63,64]. Acting cells were considered to react to the stimulus and either synthesize or reabsorb bone to stay without action according to the reference energy and the density of strain energy at the point where the cell is located.
The signal that emits the osteocyte will be assumed that decays exponentially with distance from it, meaning that the signal's radius of effect is finite. Also, the signalling will be considered to decay exponentially with time; the signalling is not a zero-one response but creates a "memory" effect of the past strain states for a finite period of time. The mathematical model of remodelling proposed is very dynamic and tends to emulate the true working of the biology of the process.
A new mathematical model is proposed that takes into consideration previous suggested ideas, incorporates new ones and meets new challenges. But most importantly, the model proposed has been tested and verified with medical observations and examined to predict normal bone density distribution of mandibular, and the distribution in some edentulism cases. In the literature generally considered stimuli is of the type where the stimulus St x depends on some function ϕ of the density of strain energy W s . A new nomenclature needs to be introduced before getting into more detail. Osteocytes receive a mechanical stimulus and they produce signalling that activates the functioning of the osteoclast or the osteoblast. The stimulus that the osteocytes receive could depend on many factors, in our case, a strain energy dependence is considered. The signalling that the osteocytes produce is considered to depend on the type of stimulus (activating either osteoblast or osteoclast), how strong the stimulus is, the relative distance to the osteoblast or osteoclast and time. So, a signalling E x is introduced where: − ξ a point where the signal is produced (position of the osteocyte), − x a point where the signal is received (position of the osteoclast or the osteoblast), − Ψ the stimulus which depends on the density of strain energy and time, − R the distance from where the signalling is produced, − X a function which depends on the bone's density or bone's quality at the point ξ where the signalling is produced, − Φ a function that determines how signalling decays with time, in our case we considered to be exponential.
Having included this broad approach to the study particular cases could be considered, for example, signalling is faulty, the bone's quality is not perfect, or the signalling does not decay exponentially and many others.
Taking into account the considerations previously mentioned about the signalling, E x could be written as where ρ denotes tissue mass density.
In this formulation, we assume that biomechanical processes can be considered quasistatic. The effects of infinite memory associated with integral over time span < 0, t > are included, but it is assumed that only finite past time interval is considered since the older effects decay exponentially, their effect could be neglected. Therefore the standard, linear static elastic equations are used. Equilibrium equation Linear strain tensor

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 mass density.
Input the new mass density and assosiated material parameters as initial values.
Run the model in Comsol again.

Fig. 1 Flow chart of the script
General stress-strain relation where tensor of material parameters depends on microstructure of tissue (in our formulation a simplified case is assumed where C i jkl depends on tissue mass density)

Materials and methods
To perform a numerical validation of the proposed model, the discussed earlier mathematical relations were implemented in COMSOL Multiphysics , which is widely used to solve partial differential equations using finite element method, combined with MATLAB ver. R2018b to control iterations representing time steps.
The following initial material characteristics were considered for calculations using in COMSOL Multiphysics : The COMSOL calculations were controlled by a MATLAB ver. R2018b script according to the flow diagram shown in Fig. 1.
As discussed in the previous section, the signal 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. In reality, after a certain time, the signal diminishes and its effect can be neglected. Therefore, the finite memory is assumed and considers a finite period of time in the past affecting the present state of bone.
To show details of the memory effect, a 2 cells model is considered. In the next graphic, the interaction of 2 cells emitting a signal according the strain energy density is shown. A memory effect was added where the signalling at any given moment in past decays exponentially on the following time steps. So the signalling at any moment is the integral of signalling over the vicinity of the acting cell plus a time effect. The model was run 100 time steps and an example of distribution of signal from two neighbour cells at some intermediate moment is presented in Fig. 2.
The density of strain energy W 0 that maintains the equilibrium between resorption and synthesis of bone was calculated with the following relationship where W s is the density of strain energy over the bone's domain, and V denotes the total volume considered (bone/sample/mandible). The integral is calculated over the whole mandible. W 0 was calculated only at the beginning of the process, and it was maintained constant throughout considered process of 150 time steps during numerical calculations. Following formula (3), outside the lazy zone, the mathematical relationship between the stimuli, density of strain energy, the density of the bone and the distance of the signalling is where E n is the signalling at the step n (time t n ) for a particular osteocyte, W s the density of strain energy, ρ the bone's density and R the distance travelled by the signal. The parameter d was calculated considering that at the average distance between 5 cells, the stimulus dropped 90% of its peek value. The stimulus does not disappear instantly but decreases with time, adding it contribution to the next stimulus felt by the actor cells in the future coming time steps. The total stimulus E s is calculated by Time is considered mainly as a parameter controlling slow evolution of the actual values of tissue's biomechanical characteristics. This is worth mentioning that geometrical state of bone is considered dependent on history of changes what makes this formulation more reasonable. In numerical calculation, the time domain is divided into time steps and integral over time domain is replaced by a sum of a finite number of steps. A common approach to treat time dependence is the use of Prony where γ is the relaxation coefficient and α is a constant. Relaxation and creep studies use normally these type of series, and they are very useful to study timedependent strain. Similarly, stimulus was considered not to relax instantly but with exponential time decreases making Prony series ideal for its study.
Using the Prony series with four time steps as finite memory, and coefficient of relaxation equal to 1, the total stimulus could be written as We consider ε to be small enough to be neglected so we keep only the first two terms, which considers the current stimulus plus the stimulus of previous 4 steps decaying exponentially. The contribution of the 5th time  [6,8,13,25] step, using a coefficient of relaxation equal to 1, is < 1 per cent of the total value which is small enough to be neglected.
The ratio of density change was calculated by where F a is the ratio of density change given in percentage. It was possible to modulate how big was the maximum ratio of change for every cycle with the parameter a. It was also considered the maximum ratio of change 1 and the minimum −1. So, the parameter a was chosen accordingly. The signalling has the following relationship, very similar to the one used previously in the literature where dzp is the lazy zone upper limit, and the dzn is the lazy zone lower limit. In the relationship considered, it is not needed a third case, as previously done in the literature [36], because the F a changes sign accordingly to the stimulus; if stimulus is above the value of W 0 , F a has a positive sing signalling bone's synthesis; and if the stimulus is below this value, F a is negative signalling bone's resorption.
In numerical simulation the following relation between the density ρ and the Young's modulus, E, was assumed during the process of bone remodelling E = bρ 2 (13) parameter b was chosen to obtain the values considered in Table 1 for a normal bone. The dimensions to use in our geometrical model were taken from tomographic examination of a real mandible. The place for the muscles working on the mandible was localized on our model, and the model was run over COMSOL Multiphysics using MATLAB ver. R2018b for 150 time steps. The considered 2D geometrical model is presented in Fig. 3, and it shows the initial density distribution for all the situations analysed to enable comparison between cases with the same initial distribution of material and loads but different configurations of teeth. It is worth mentioning that one side of the model in considered examples was always left unchanged as a reference while the other one was modified by removing selected teeth.
During bone's strength analysis, a crucial factor is to choose the way bone's is loaded similar problems could be seen in the literature of structures [61]. The actual place of the muscles was considered as locus of the applied forces.
Simplified 2D model of mandible was loaded in three different places, at the localization of the temporalis muscle, at the localization of the masseter muscle and on the teeth. Forces at the teeth were considered according to values reported on the literature [18,26,27,39] and the forces applied by the muscles where calculated accordingly to make total momentum zero over the mandible [38,50,54] and: When F F F m is the muscle's vector forces, F F F c is the reaction vector forces at the Condyle, and the F F F b is the bite's vector force, and r stands for the action radius of each force. As a boundary condition, COMSOL Multiphysics roller feature over the external boundary was added at the Condyle simulating its actual biomechanical function (slip between head and cup of the joint).   Fig. 12.
In the first of these three examples, a case was considered with missing the 8th and the 9th teeth. In the following figures, a sequence of time steps stages is presented to show how the remodelling process results in redistribution of the tissue. The model after 50 steps changed to: After 100 steps: The model after 150 steps changed to: After 130-150 time steps the mass density reached a stable state with stimulus falling mainly into the "lazy zone" in most of the bone's area, and no significant further changes were noted. At this stage, the calculations were ended.
For a comparison the picture from X-ray examination is presented in the following graphic, Fig. 7:  Areas with high bone density are marked in Fig. 7 with an H and areas with low bone density are marked with an L. In Fig. 6 the model is shown after 150 steps and a comparison between areas of high/low bone density could be done between both images showing that our model predicts the areas of growing bone (high density) and the areas of resorbed bone (low density).
For further validation, two teeth from considered model were removed, the model was tested 150 time steps for similar reason than the previous case, mass density distribution had reached equilibrium, and images  For a comparison of the results of numerical calculation the X-rays image is presented in Fig. 11.
Comparing Figs. 10 and 11, typical areas of high and low bone density can be observed in both pictures (marked with H and L, respectively). There is also agreement in both pictures regarding the location of the areas where the tissue has been significantly weakened due to a lack of teeth (marked with an R).
Another case considered for validation with different configuration of teeth missing is displayed in Fig. 12.
In the following three figures, Figs. 13, 14, and 15, the results of computer simulation after 50, 100 and 150 steps are shown.
Comparing Figs. 12 and 15, a match between the areas of high and low bone's density in both pictures can be observed. As well it can be identified an area of bone resorption of the bone due to lack of teeth. This case is more remarkable because in both figures it could be seen that the teeth on the side of the missing teeth area and alone is loosing the bone's grip.
In the area studied, where the teeth are missing in the X-rays and in the model, the resorption behaviour of the bone was predicted correctly. The effects of the missing teeth are also predicted to be over the bone localized in the vicinity. The loosing of bone material around the teeth that is left alone at the left end of the mandible is correctly predicted as well in proposed model.

Conclusion
A new model with memory effects was presented which describes the synthesis and bone resorption considering a strain-based stimulus that provokes a decaying exponentially in time and space signal activating the cells involved in the natural processes of the bone. The model takes into account the bone's density influence on the stimulus on the signalling, as stated before, more density means more active osteocytes signalling. In the areas where density is lower the stimulus to grow bone is lower and vice versa. This effect is mostly omitted in the literature.
The fact that stimulus decays exponentially throughout time is the new feature proposed in the model. Models of interaction between bone and teeth have been proposed in the literature but in the model proposed here formulation of stimulus and signalling are treated separately. The consideration of the signalling and its decays in time and distance makes the model very "dynamic" and useful to study the lack of teeth, and possible diseases of the bone related to faulty signalling between cells.  The presented model correctly predicts the synthesis and bone resorption in areas where the teeth are missing and it was validated by the medical observations using pantographic X-rays pictures of dental patients. However, considerable changes of bone density may lead to significant tissue gradient deterioration and fatigue and damage analysis could be helpful in future more advanced studies [17,53]. In the dental implant vicinity, in the case of loosening implant, another material model should be applied for instance granular material, see Misra and Poorsolhjouy [47].
The current model will be tested in cases of analysis of the interaction between bone substitute material, living tissue, and implant. However, in order to study the behaviour of bone substitute material and implant, different constitute models should be included in our formulation, for example, models of truss and pantographic-like microstructures, see, e.g., [1,11,16,21,31,56]; models including micropolar theory, see, e.g., [2,[22][23][24]; and models including higher gradients, see, e.g., Dell'Isola et al. [19,20].