On the effect of antiresorptive drugs on the bone remodeling of the mandible after dental implantation: a mathematical model

Bone remodeling identifies the process of permanent bone change with new bone formation and old bone resorption. Understanding this process is essential in many applications, such as optimizing the treatment of diseases like osteoporosis, maintaining bone density in long-term periods of disuse, or assessing the long-term evolution of the bone surrounding prostheses after implantation. A particular case of study is the bone remodeling process after dental implantation. Despite the overall success of this type of implants, the increasing life expectancy in developed countries has boosted the demand for dental implants in patients with osteoporosis. Although several studies demonstrate a high success rate of dental implants in osteoporotic patients, it is also known that the healing time and the failure rate increase, necessitating the adoption of pharmacological measures to improve bone quality in those patients. However, the general efficacy of these antiresorptive drugs for osteoporotic patients is still controversial, requiring more experimental and clinical studies. In this work, we investigate the effect of different doses of several drugs, used nowadays in osteoporotic patients, on the evolution of bone density after dental implantation. With this aim, we use a pharmacokinetic–pharmacodynamic (PK/PD) mathematical model that includes the effect of antiresorptive drugs on the RANK/RANK-L/OPG pathway, as well as the mechano-chemical coupling with external mechanical loads. This mechano-PK/PD model is then used to analyze the evolution of bone in normal and osteoporotic mandibles after dental implantation with different drug dosages. We show that using antiresorptive agents such as bisphosphonates or denosumab increases bone density and the associated mechanical properties, but at the same time, it also increases bone brittleness. We conclude that, despite the many limitations of these very complex models, the one presented here is capable of predicting qualitatively the evolution of some of the main biological and chemical variables associated with the process of bone remodeling in patients receiving drugs for osteoporosis, so it could be used to optimize dental implant design and coating for osteoporotic patients, as well as the drug dosage protocol for patient-specific treatments.

www.nature.com/scientificreports/ of denosumab and Ibandronate. As shown in 34 and in Fig. 3, the ash fraction decreases during the first stage of remodeling, because of the activity of osteoclasts tends to increase. When applying one of those two drugs, the ash fraction increases for all types of bones. In the cases of osteoporotic bone under σ = 1.0 MPa , and trabecular bone, the ash fraction first increases, and then decreases for doses of 1.0 and 3.0 mg/kg of denosumab. After this period, the ash fraction begins to increase again. The bone volume fraction increases for all types of bones, as shown in Fig. 3. This same figure also shows that the bone volume fraction tends to decrease with the stress level in trabecular and osteoporotic bones. The same trend can be observed when prescribing bisphosphonates. The bone volume fraction in osteoporotic bone increases in that latter case up to the maximum density allowed. Also, the increase rate for the bone volume fraction is higher when using denosumab than when applying bisphosphonates (Fig. 3). Finally, the damage level in cortical bone ( ρ = 2.05 g/cm 3 ) shows a greater increase than for the two other bone types. Using denosumab with doses of 0.1 and 0.3 mg/kg does not increase the damage level in osteoporotic and trabecular bones, while, on the contrary, for cortical bone, for those doses, the damage level increases. In the case of bisphosphonates, doses of 0.25 and 0.5 mg do not increase the damage level for osteoporotic bone. The change in ash fraction for ρ = 1.0 g/cm 3 subjected to σ = 7.0 MPa for for 0.25, 0.5, 1.0 and 2.0 mg of denosumab with respect to the control case at time t=3000 days was about -2%, 14%, -4% and 20% respectively, while these changes for 0.25, 0.5, 1.0 and 2.0 mg of Ibandronate were about -0.7%, -0.9%, -1.5% and 5% respectively. The evolutions of bone volume fraction, ash fraction, and damage for different doses of drugs and different dosage intervals for the osteoporotic bone, in particular, are depicted in Figs. 4 and 5. For a denosumab dose 0.1 mg/kg, increasing the dosage interval decreases the bone volume and ash fractions. In contrast, for a dose of 0.3 mg/kg and a time interval of 60 days, the damage reaches the maximum level allowed. After reaching that maximum damage, the ash fraction shows a higher reduction rate than for other doses. Simultaneously, the volume fraction for this dose is slightly reduced and then reaches a new steady value after the maximum damage is reached. For other time intervals and the same dose, the bone volume fraction increases with time, depending on the time interval, since damage does not reach its maximum level. In contrast, the ash fraction decreases when increasing the time interval. Increasing the drug dose increases damage, while shorter time intervals drive to a faster increase in damage up to its maximum level. Moving now to Ibandronate, a dose of 0.25 mg increases the damage up to its maximum value at time intervals of 60 and 90 days. The ash and bone volume fractions, unlike denosumab, do not stop increasing after reaching maximum damage. The same result is observed for other doses and dosage time intervals. It is possible to observe faster increases in the ash and bone volume fractions when decreasing the dosage time interval for all doses. Reductions in volume fraction for ρ = 0.5 g/cm 3 after 3000 days from the case of 60 days for 90, 120, 150 and 180 days of administration of 0.1 mg/kg of denosumab were about 13%, 31%, 37% and 37% respectively, while for Ibandronate these values were 45%, 45%, 56% and 57% respectively. Figure 6 shows the final density distribution in the whole mandible after application of the phenomenological bone remodeling model with additional views of several cross-sections to compare such results with corresponding CT images. Figure 6f shows a cut view of the density distribution of the reduced computational model, while Fig. 6g depicts the assumed density distribution after applying the density reduction due to osteoporosis.
Finally, the application of the coupled PK/PD and remodeling models here described to the mandible model after dental implantation, physiological mastication loads, and administration of different drug doses drives to the results shown in Figs. 7, 8, 9 and 10. The bone volume fraction (Fig. 7), damage level ( Fig. 8) and ash fraction ( Fig. 9) are compared with the base case with the application of no drug. For the two drugs analyzed here, any dose increases the bone volume and ash fractions and corresponding density. Doses of 1.0 and 3.0 mg/kg of denosumab and 1.0 and 2.0 mg of bisphosphonate produce higher damage. As consequence, the elastic modulus decreases in those regions (Fig. 10). The ash fraction ( Fig. 9), like damage, increases when increasing the drug dose, especially in the trabecular bone around the implant threads.

Discussion
One of the main objectives of this study was to investigate the effect of different doses of antiresorptive drugs on bone behavior. Antiresorptive drugs are used to treat diseases such as osteoporosis in which the balance between the activity of osteoclasts and osteoblasts is disturbed 37 . Denosumab and Ibandronate are antiresorptive osteoclast-targeting drugs used in the treatment of osteoporosis 25,38 . These drugs affect the bone remodeling process by different action mechanisms. Denosumab binds to RANK-L, reducing the binding between RANK and RANK-L, thereby the concentration of active osteoclasts on the bone surface. Bisphosphonates as Ibandronate, on the contrary, interfere with the osteoclast activity by binding to the bone mineral surface 25 . As a consequence, these drugs reduce the resorption stage in bone remodeling, increasing, therefore, the bone volume fraction and density, and, with that, improving the long-term bone quality in osteoporotic patients. Besides, the mineral content of bone increases, which causes bone to become more brittle and damaged. This may provoke local fractures despite the higher stiffness and strength of the treated bone.
Therefore, this ambivalent effect of antiresorptive drugs makes it difficult to predict their net effect on osteoporotic bone. This process is especially complex when treating with such drugs after implantation since in those case, it is not only the effect of drugs but also the critical change in the mechanical conditions of the surrounding bone which contributes to modifying the long-term bone internal microstructure. Implantation in patients with osteoporosis is, therefore, challenging 39 and its clinical treatment utilizing these drugs may negatively affect the success of the implant with increasing osteonecrosis 24,38 . This happens, for example, after dental implantation, a practice that has increased in recent years in the elderly, who have an increased risk of osteoporosis in the mandible bone, which justifies why this problem has attracted the interest of several authors 14 www.nature.com/scientificreports/ problems. In principle, an ideal mathematical model of bone remodeling should take into account the different bone cells involved and their main activities such as proliferation, differentiation, migration, death, biochemical signals production, changes in their expression due to biochemical or mechanical signals, and tissue resorption or www.nature.com/scientificreports/ production. Also, it should consider the diffusion, decay and production of growth factors, hormones, proteins, and other biochemical substances that control the cell behavior. Finally, the particular mechanical microenvironment and its interactions with the chemical reactions and cell protein expression should be taken into account. To implement all these processes requires very complex mechano-chemo-biological models with several coupled mechanisms not yet fully understood, a lot of parameters, many times unmeasured, and a difficult validation due to lack of experimental results in a sufficient number and variety of situations. Even with these limitations, mathematical modeling is a powerful tool for studying complex biological systems since they allow us to find out important trends, and to quantify, to a certain extent, the relationships between causes and effects, to test theoretical hypotheses, quantify the effects of the different parameters individually on the behavior of the biological system and to do virtual experiments in "what if " situations 42 .
In this paper, we present a combination of a PK/PD model and a fully-coupled chemo-mechano-biological bone remodeling approach that incorporates the stimulus effect on the signaling pathway between osteoclasts and osteoblasts. The effect of damage on the signaling pathway and the local material properties have also been considered. Finally, the mineralization level is monitored along the whole bone lifetime. We have proven that these types of models can be used for predicting the bone behavior in the mandible after dental implantation when using antiresorptive drugs for improving the long-term quality of osteoporotic bone. In order to study the applicability of this model, the long-term effects of different doses of denosumab (0.1, 0.3, 1.0 and 3.0 mg/kg) and Ibandronate (0.25, 0.5, 1.0 and 2.0 mg) on the mandibular bone surrounding a dental implant were studied.
We first examined the effect of different drugs and of the mechanical stimulus on the bone behavior. Denosumab does not bind to the bone mineral surface, unlike Ibandronate, so Ibandronate effects last longer after stopping its administration 24,43,44 . Considering this fact and comparing the increase in bone volume fraction after administration of these two drugs, we can speculate that the increase in volume fraction and in mineral content in bone induced by Ibandronate will be higher than that of denosumab. As commented, the drugs here analyzed inhibit the activity of osteoclasts, thereby reducing bone resorption and therefore, increasing the bone volume fraction and the mineralization. On the other hand, due to the effect of calcium, the bone becomes more brittle, and damage increases. By comparing Figs. 4 and 5, at low doses of 0.1 mg/kg for denosumab and 0.25 mg for Ibandronate, we found a more significant increase in the bone volume fraction when applying Ibandronate, while, contrarily, this drug produces higher increases in bone ash fraction and damage than denosumab, leading to a more brittle bone. Also, these drugs increase the ash fraction initially (Fig. 3), while the mineralized portion of the bone is also initially reduced. This stage is then followed by the filling of the resorbed bone by the  www.nature.com/scientificreports/ osteoblasts, forming new bone with the corresponding next mineralization. Increasing the drug dose in any of two cases considered also increases the bone volume fraction and damage for all types of bone. For cortical bone, for instance, which does not need drug treatment, even low drug doses cause bone to become highly brittle, reaching the maximum level of damage faster. In this bone, no decrease in the mineralized fraction was detected after drug treatment, so microdamage progresses in time, and a stress fracture may occur. The stress level also influences this behavior by promoting or delaying the damage rate. For trabecular or osteoporotic bone, low drug doses increase the bone volume fraction without substantial damage increase. This is more evident in osteoporotic bone. Consequently, drug application is beneficial in low-density bone, although this treatment always increases brittleness, which may compromise the success of dental implantation. This subtle control between these two opposite effects is essential when defining the treatment protocol  www.nature.com/scientificreports/ for a particular patient and the main reason for using and improving this type of mathematical models and "in silico" experiments. After this first test of the model, when we got similar results to other authors with comparable qualitative clinical conclusions 20-23 , we studied their effect on the bone surrounding a realistic dental implant using a 3-dimensional finite element model for an osteoporotic mandible after implantation and treatment with different drug doses. Qualitative comparison of the results with clinical studies in patients with osteoporosis who received bisphosphonate or denosumab agrees with the detected bone mineral content and bone mass increase after anti-resorptive treatment 45,46 .
Bone brittleness is not only linked to low bone mass but also to the increase in the mineral content of the tissue or to the accumulation of bone microdamage 47 . This latter could be due to bone resorption and the associated bone strength or to microcrack blockage of the bone remodeling activity, reducing microcracks repair 48 . Another  www.nature.com/scientificreports/ possibility is associated with the use of bisphosphonates or denosumab, since these drugs reduce the activity of osteoclasts, thereby reducing the bone remodeling rate, thus slowing down microdamage repair. As shown in our results, depending on the stress level, the bone type, and the drug dose, the damage level could highly increase when using these types of drugs. For example, in Fig. 3 for osteoporotic bone ρ = 0.5 g/cm 3 , after applying constant stress of σ = 1.0 MPa the amount of volume fraction increased along the whole process (43% for 0.3 mg/kg of denosumab and 30% for 0.5 mg of Ibandronate), but the increase in ash fraction in the early days leads to a maximum in the damage level. A similar effect is seen in Figs. 7 and 8. The first one shows an increase in the bone volume fraction in the patient with osteoporosis for higher doses of antiresorptive drugs. This increase is also seen in the vicinity of the implant threads. On the other hand, the second figure shows that microdamage also increases, especially in the threads neighborhood, which confirms the findings of other studies 47,49,50 .
Despite the remarkable potential of these models and the qualitatively accurate results produced, like any other mathematical models in systems biology, it still has important limitations that should be addressed in future studies. Among them, we can mention: 1. As mentioned in our previous study 34 , simplifications such as considering only the RANK/RANK-L/OPG biochemical pathway, or in the reaction rate assumed for the receptor-ligand reactions, or discarding the effect of osteocytes in the mechanosensing process, should be progressively overcome. 2. Due to the lack of knowledge about the effect of race, sex, age, etc, onto the biological parameters, those used in this study are taken from literature, so it is not possible yet to apply the model to patient-specific studies without additional experimental determinations. 3. The pharmacokinetics determines the dosage of the drug in the intended compartment. Therefore, developing a pharmaco-kinetic model with an adequate number of compartments and experimental parameters is essential for determining an accurate value for the concentration of drug reaching bone to better predict the bone behavior. Also, the parameters used in pharmacokinetic models should vary for different individuals, which has not been considered here. 4. Although the mineral part of the bone has a prevailing effect on the bone mechanical properties 51 , the effect of drugs on the organic component of bone was not considered, which may be important to assess bone fracture resistance, especially under tension 52 . 5. Different loading conditions have different effects on bone behavior. However, due to the lack of sufficient information, and considering that muscle forces vary for different individuals and are also modified after implantation, the loads applied in our finite element study were based on previous generic studies.
We have presented here a PK/PD model coupled with a mechano-chemo-biological bone remodeling model, that has been implemented in a standard commercial FE code to analyze the effect of two widely used drugs as denosumab and Ibandronate onto the evolution of the mandible bone after dental implantation. Using antiresorptive agents such as denosumab or bisphosphonates alter the evolution of the density and mechanical properties of bone by interfering in the bone remodeling mechanism, reducing the osteoclast activity, which may consequently increase bone brittleness by augmenting the mineral content and microdamage. This effect is corroborated by our results that show that using any of those drugs in osteoporotic patients increases the bone volume fraction, although, in parallel, it also increases bone brittleness, a well-known side effect of these treatments. Using patientspecific geometries and initial density values, this model may provide a good perspective to clinicians about these two contradictory effects of those drugs for the treatment of osteoporosis and to compare with longitudinal clinical results with different doses.
The model developed here is capable of capturing the main biological and chemical variables associated with the bone remodeling process in patients receiving drugs for osteoporosis treatment. Despite the limitations described, these models can be used to predict bone behavior, complementing costly, and time consuming clinical experiments, and for optimizing dental implant designs and coatings for osteoporotic patients, as well as patient-specific dosage protocols.

Methods
The bone remodeling model used here is a chemo-mechano-biological model that couples the effect of mechanical strains and microdamage with the biochemical RANK/RANK-L/OPG pathway, and, with the expression of the different cell phenotypes involved. It also models the tissue resorption-formation process accomplished by synchronized sets of osteoclasts and osteoblasts known as basic multicellular units (BMUs), followed by the tissue mineralization. This model closely follows the one presented in a previous paper of the authors 34 but adding now the PK/PD submodel that permits the analysis of the effect of drugs. As a result, the whole process and the effect on its main output variables, such as the local bone density and the mineral content, as well as their influence onto the bone mechanical properties, can be analyzed for different mechanical or chemical protocols. π TGF β and π RANK−L are the fraction of receptor sites occupied by receptor/ligand complex (i.e. π ∼ R • L/R ) of the transforming growth factor ( TGF β ) and of RANK-L, respectively. D R , D B and D C are differentiation rates of osteoblast progenitors to responsive osteoblasts, responsive osteoblasts to active osteoblasts, and osteoclasts precursors to osteoclasts, respectively. D A is the osteoclasts apoptosis rate caused by TGF β and, finally, k B is the death rate of active osteoblasts. The values of all these parameters appear in Table 1 with the corresponding references. After calculating the osteoblast and osteoclast populations, the rate of bone volume fraction can then be calculated as: where k form and k res are the rates of bone formation and resorption per unit cell, respectively, that depend on the bone type and location. The initial value of their ratio k form /k res is assumed to be equal to the initial ratio of the population of osteoclasts and osteoblasts C 0 /B a0 , that can be obtained by solving the steady-state expression of Eq. (1) (left-hand side equal to zero) and depend on the initial bone density.
Equations (1) and (2) provide the evolution along time of the distributions of the different cells involved in the process depending on the results of the chemical reactions between receptor and ligands involved in the RANK/RANK-L/OPG pathway that is briefly explained in the following section, and from them, the change in bone volume fraction.

Chemical reactions involved in bone remodeling. The receptor activator for nuclear kappa-B, RANK
(K) is a surface-bound molecule that binds to its ligand, RANK-L (L), serving as osteoclast activator 28 . Osteoprotegerin, OPG (O), a decoy receptor for RANK-L, is another protein expressed by osteoblasts and other tissues like the spleen, bone marrow, heart, liver, and kidney 26 . OPG inhibits RANK/RANK-L binding, so it plays a protective role against bone loss. Among the many systemic hormones that influence bone cell activity, PTH (P) is the most important calcium homeostasis regulator and bone remodeling hormone, so it is used as an anabolic agent in the treatment of osteoporosis 27 . The chemical reactions between the receptors and associated ligands are as follows 27 : www.nature.com/scientificreports/ with k i (i = 1, 3, 5) reaction binding rates, and k i (i = 2, 4, 6) reaction unbinding rates. p O , p L and p P are production fluxes and d O , d L and d P are destruction fluxes of OPG, RANK-L and PTH respectively. Since the reactions related to chemical bindings in biological systems occur faster than cell population changes, occupancy of the complexes ( π L = K•L K and π P = P r •P R P T = P P+P s ) as well as the OPG concentration (O) have been considered as pseudo-steady during the whole remodeling process. With all this, the OPG concentration and the fractions of RANK-L and PTH can be calculated as follows 27 : where P r is the number of free receptor of PTH, R P T stands for the number of PTH receptors per cell, P s = k 6 /k 5 , K P L is the maximum occupancy of RANK-L attached to each cell surface, P = S P k P is the PTH concentration, S P is the synthesis rate of systemic PTH, K P O represents the minimal rate production of OPG per cell and k O is the rate of OPG removal.
From Eqs. (1)-(4) we can calculate the evolution of B r , B a , C and v b without considering the effect of the mechanical environment on bone signaling. The next sections address this coupling between Mechanics and Biochemistry in bone remodeling.

Role of mechanical signal in chemical reactions.
Mechanical strain and microdamage are the most important signals affecting the receptor-ligand binding/unbinding energy fraction. In particular, mechanical forces change the energy barrier of the molecules and disrupt the binding/unbinding of receptors and ligands in chemical reactions, while damage interrupts the communication channels between cells, thus reducing the level of the mechanical signal. Therefore, we hypothesized in a previous paper 34 that the reaction rates of the receptor and ligand change as follows: where k r and k f are the unbinding and binding rates associated to the receptor and ligand interactions, k r 0 and k f 0 are unbinding and binding rate constants, γ is a constant relating the ratio between current and initial bone volume fractions ( γ = v b v b 0 ), S the normalized mechanical signal (see below) and Ŝ v is the normalized specific bone surface ( Ŝ v = S v S vmax , with S v max = 4.17mm 2 /mm 3 the maximum specific bone surface available) 53 . The mechanical signal (that is here considered as inhibitory of the action of the cells 54,55 ), depends on the mechanical stimulus ( ξ ) and on the damage level (d). We assume that the transmission of this signal in bone depends on its microstructure through the local bone volume fraction (i.e. v b ). In other words, that signal is transmitted easier in cortical than in trabecular bone and much easier than in osteoporotic bone. With all this, and following 55 , we can write: with c and a model parameters that, together with ξ , d and v b control the value of the inhibitory remodeling signal. The mechanical stimulus depends on an equivalent strain ε and the number of loading cycles N as in 56 : 55,57 , being u the strain energy density and E the elastic modulus, that correlates with the bone volume fraction and ash density as stated in the Eq. (21) below. Finally, m = 4 is an experimental constant 56 and i is the number of different load types. Drug effect on bone remodeling. Different therapies for bone diseases such as osteoporosis have been developed, influencing the RANK/RANKL/OPG pathway. For example, denosumab is one of such drugs that binds with RANK-L, reducing its binding potential with its receptor RANK 25 , and with that diminishing the differentiation to osteoclasts and their subsequent activation. In our model, this effect is modelled by adding the following additional reaction in Eq. (3) 33 : www.nature.com/scientificreports/ where C d is the drug concentration in the plasma compartment, and k on and k off are the association and dissociation rate constants of that reaction, respectively. After inclusion of this reaction, π RANK−L changes in (4) such as 33 : with K P L the maximum number of RANK-L attached to each cell surface 27 . Another drug usually used to treat osteoporosis is Ibandronate. Ibandronate, like other bisphosphonates, adheres to the bone mineral, preventing the activity of mature osteoclasts 25 . As a consequence, the osteoclast population (C) in Eq. (1) decreases 8 . In our model, we take this into account by adding the following equation: where D A and D C are the osteoclast apoptosis rate caused by TGF β and the differentiation rate of osteoclast precursors to osteoclasts, respectively, as defined in Eq. (1), I max is the maximal fractional extent of inhibition, IC 50 the concentration producing 50% of maximal inhibition and C b the concentration of bisphosphonate.
To calculate the drug concentrations ( C d and C b ), a pharmacokinetic-pharmacodynamic (PK/PD) model has been used 8,33,35,58 . Following 33 , the PK/PD model equations for subcutaneous administration of denosumab may be stated as: This model considers two compartments (subcutaneous and plasma) as shown in Fig. 11a. The first two equations determine the drug pharmacokinetics providing the drug concentration in the plasma compartment. In contrast, the third equation determines the pharmacodynamics in the action site, which reflects the maximal drug effect on the bone tissue. In these equations, Dose corresponds to the denosumab dose administrated subcutaneously during each time interval, C tot is the total drug concentration, sum of the free drug ( C d ) and the drug-ligand complex ( L • C d ), k a is the first-order absorption rate of the drug administrated subcutaneously into the plasma compartment, V c /F is the bioavailability-adjusted central compartment volume, k int is the drug-ligand complex internalization or degradation, k el the drug removal rate from the central compartment, R ss the concentration of free ligand in the steady-state situation, K D = k off /k on the equilibrium disassociation constant for the drug-ligand reaction (Eq. (10)), NTX the crosslinked N-telopeptide of collagen type I, a bone turnover biomarker measured in the serum, which determines the drug response in bone and, finally, k in and k out are the production and removal rates of NTX, respectively. The steady-state values of NTX ( NTX ss = k in /k out ) for multiple myeloma patients receiving 0.1, 0.3, 1.0 and 3.0 mg/kg of denosumab were chosen from 33,36 who gave values of 9.8, 9.1, 10.2 and 8.1 nM respectively.
Similarly, the PK/PD model for bisphosphonates as Ibandronate may be written as 35 : The first four equations establish the four-compartment pharmacokinetic model for bisphosphonates 35 , providing the Ibandronate concentration in plasma, bone and two peripheral compartments (see Fig. 11a). In these equations, the subscripts pl, p1, p2 and b stand for plasma, peripheral1, peripheral2 and bone compartments, respectively. C i (i = pl, p1, p2, b) is the drug concentration in the compartment i, Q i (i = p1, p2, b) is the intercompartmental clearance, V i (i = pl, p1, p2, b) the volume of each compartment, and CL the renal clearance. www.nature.com/scientificreports/ Dose corresponds to the Ibandronate dose administrated intravenously during each time interval and C pl is the sum of current drug concentration and the amount administrated at the beginning of the dosage interval. The fifth equation shows the pharmacodynamics of Ibandronate in bone, with uCTX the Urinary C-Terminal Telopeptide of Collagen Type I, a bone turnover biomarker measured in urine, KS and KD the uCTX formation and degradation rates, respectively, R tar the limit value for the uCTX formation rate defined by the rate constant of k qq , IC 50 the Ibandronate concentration in the bone compartment producing 50% of maximum response of uCTX and n the Hill coefficient 35 . The PK/PD parameters used in this study are shown in Table 2, with the corresponding references.
Mineralization. Newly created bone tissue is mainly composed of collagen osteoid that is, then, progressively mineralized. Mineralization is generally divided in two phases: the first one is fast (it lasts only several days) reaching about 60% of the maximum mineral content in mature bone, while the second is much slower, lasting several years to achieve the final mineral content in fully developed bone 55,57 . Here, as in our previous model 34 , and in 55 , the first phase was assumed as instantaneous, following bone mineralization an exponential evolution during the second stage, that is: with α max and α 0 the maximum and initial ash fractions (bone residual after calcination that essentially corresponds to the mineral content), respectively, and κ the exponential mineralization constant. Taking this into account, and for a representative volumen element, the mean ash fraction at a certain time after periods of new bone formation and resorption may be written as 55 : with v b,0 the initial bone volume fraction, while, from Eq. (1), v F = dB a dt and v R = dC dt are the rates of bone volume fraction of newly created and resorbed bone, respectively and h 0 is the initial microcrack density defined as the number of microcracks per unit volume. Bone, as a living tissue, is able to repair those micro-cracks, so damage increases when having high stresses/ strains, ḋ acc (damage accumulation rate), while, at the same time, microcracks are removed in regions where bone is resorbed, ḋ rep (damage repair rate) 2 . We can then write: As stated in 34 the damage accumulation for a certain number of cycles is a function of the load amplitude and the type of stress state (tension, d acc,t , or compression, d acc,c ) can be written as: which δ 1,2 , γ 1,2 and C 1,2,3 are parameters. N is the number of load cycles and ε = √ 2u/E is the equivalent strain in each of those cycles, which described before, with u the strain energy density and E the elastic modulus. E * is the the reference elastic modulus which for undamaged ( d = 0 ) cortical bone the ratio E/E * is equal to one.
Fatigue life ( N f ) in compression and tension calculated as:  where δ i and β are constants. ε u is ultimate strain which has relation with calcium content: Finally damage repair evolution is calculated as follows: where v R = k res C is the rate of bone volume fraction due to osteoclasts activity.
Mechanical properties. Finally, as a first approach, and despite the well-known local orthotropy of bone tissue 61 , we assumed bone tissue as heterogeneous and isotropic with its mechanical properties defined by the following correlation between the volume fraction ( v b ), ash fraction ( α ) and damage (d), with the bone elastic modulus 55,62 : A summary of the mechanical parameters used in these study are presented in Table 3.
Finally, a scheme of the chemo-mechano-biological bone remodeling model, coupled with the PK/PD models for the two drug contents, is illustrated in Fig. 11.

Finite element simulation.
Computed tomography images of a healthy adult woman were used to construct the 3D geometric model of the mandible. After segmenting the mandible and teeth in MIMICS 10 (Materialise, Leuven, Belgium) and obtaining the STL files, CATIA (CATIA V5, Dassault Systèmes, Vèlizy-Villacoublay, France) was used to create the final three-dimensional geometry of the mandible, teeth, and PDLs (see Fig. 12). The gaps between the mandible and the teeth were used to obtain the geometry of the PDLs. The implant selected for this analysis is based on the INTRI design without internal resilient parts 68 . The height of the implant was 11 mm and its diameters at the top and bottom were 5.1 mm and 4.5 mm with two threaded steps of 1 and 0.5 mm respectively. Finally, the crown was designed for the first right molar, taking into account the implant neck.
The model here developed is based on physiological mechanisms and properties, so, contrary to other phenomenological bone remodeling models, it leads to wrong results when the initial density distribution is not physiological and related to the initial values of the cell concentrations. Therefore, the geometry in Fig. 12 was used first to obtain the initial density distribution for the next simulations. A phenomenological bone remodeling model 56,69 was used for this purpose, considering the mastication muscles' reaction forces, while the boundary conditions were applied to each teeth involved in the mastication process based on previous studies 63, 70    www.nature.com/scientificreports/ simulation of the complete mandible (540 simulation steps), we correlated such initial density distribution with the elastic modulus point-wise, using the following correlations: E = 1736ρ 3.2 and E = 2014ρ 2.5 for cortical and trabecular bone, respectively 63 . Finally, those values were modified to take into account osteoporosis. Reductions of 33% 66% have been reported for the modulus of elasticity for cortical ( ρ > 1.2g/cc ) and trabecular bones ( ρ < 1.2g/cc ), respectively, in osteoporotic patients 17 . Therefore, we modified the initial distribution of the elastic moduli in each of the bone types considering such values.
To reduce the computer time, a cut of bone, which includes the premolar tooth and its PDL, the second molar and its PDL and the implant, was isolated to perform the next simulations. The whole geometric model, together with the density distribution (and corresponding material properties), was then exported to ABAQUS (ABAQUS 6.11, Dassault Systèmes, Vèlizy-Villacoublay, France) to perform the finite element simulations. A four-noded solid tetrahedral mesh was built in ABAQUS-CAE. The final number of elements was obtained after ensuring sufficient accuracy in a previous convergence analysis. The final resulting number of elements for the whole mandible, teeth, PDLs, implant, and crown in the final model was 860064, 233083, 30369, 46939 and 5325, respectively, while the number of elements in the model of the section cut for bone, teeth, and PDLs was 448265, 12896 and 4729 respectively.
Complete osseointegration was assumed for the bone-implant, bone-PDL and PDL-tooth interfaces. Displacement at the nodes of the mesial and distal surfaces of the cut model was imposed with values derived from the results of the complete mandible model (Fig. 12). The Titanium implant and crown materials were assumed as linearly elastic with E = 118 GPa, ν = 0.35 and E = 82.8 GPa, ν = 0.33 respectively 71,72 . Finally, the bone mechanical properties change during the bone remodeling process, so a user material (UMAT) subroutine of ABAQUS was implemented to compute such properties along the loading process according to the model described above. To simulate the effect of drug treatment in the bone surrounding the dental implant, different drug doses were used and the corresponding results compared with the control model without the drug.